Using prior information for a GAM smooth function to reduce standard errors

Question

I have some data that I want to model in a GAM. However, there are few observations, generally leading to high standard errors.

cases = c(36, 94, 101, 128, 143, 143, 143, 143, 143, 149, 149, 151, 154, 
          155, 156, 156, 156, 163, 173, 173, 174)
controls = c(3, 3, 3, 3, 6, 12, 19, 19, 19, 19, 20, 21, 21, 21, 22, 24, 
             24, 24, 25, 27, 27, 27, 28, 29, 31, 31, 31, 31, 31, 32, 32, 36, 
             36, 36, 46, 47, 47, 47, 47, 47, 47, 48, 51, 51, 51, 51, 51, 51, 
             51, 54, 54, 54, 59, 59, 62, 63, 63, 63, 63, 63, 65, 68, 68, 75, 
             76, 79, 79, 79, 79, 79, 79, 79, 81, 81, 81, 81, 81, 81, 81, 82, 
             82, 82, 82, 82, 82, 82, 82, 92, 94, 94, 101, 104, 104, 104, 106, 
             116, 116, 125, 130, 131, 134, 138, 138, 141, 143, 143, 143, 146, 
             146, 146, 147, 149, 152, 154, 154, 154, 154, 154, 154, 154, 156, 
             156, 156, 156, 156, 156, 156, 156, 156, 162, 162, 163, 165, 168, 
             173, 174, 174, 174, 178, 186, 186, 191, 192, 192, 194)
x = data.frame(pos=c(cases, controls), cohort=c(rep(1, length(cases)), rep(0, length(controls))))

library(mgcv)
plot(gam(cohort~s(pos), family="binomial", data=x))

I have additional data but I combine this directly as it is from a different source and I do not want it to directly modify the estimates. However, I do want it to reduce the standard errors where it does conform to my data. Is there a way to use this information as a prior to achieve this?

prior_cases = c(3, 6, 9, 22, 25, 46, 49, 51, 56, 57, 63, 64, 65, 69, 77, 79, 
                82, 92, 94, 95, 98, 101, 116, 127, 128, 128, 133, 134, 138, 138, 
                140, 142, 145, 149, 149, 149, 150, 151, 152, 155, 161, 161, 170, 
                173, 173, 174, 180, 191)
prior_controls = c(2, 16, 18, 31, 57, 57, 63, 65, 74, 94, 98, 123, 154, 159, 177, 
                   178)
plot(density(prior_cases), col="red")
lines(density(prior_controls), col="blue")

The prior data is generally supporting the same relationship but I cannot figure out how to use this information. This is one example but I have many more and most with more complex relationships. Any guidance on how to proceed will be much appreciated.

Thanks for the comments.

More Background:

pos = position of a rare-missense variant in a linear protein sequence cohort 1 = affected (HCM) cohort 0 = unaffected (controls)

I am trying to model mutational hotspots in Mendelian disease-genes, that is, regions in the linear protein sequence more likely to harbour a case variant. From this data it suggests the right-tail of the protein has a higher variant density but the standard errors do not give great confidence in this. I have some previous data that also kind of aligns with this from public resources: ClinVAR, common variants in population cohorts. Ideally, I would like to also use this information to give more support to the trend observed in the case-control data. GAMs may not even be the best approach here, they have been successful previously for some genes, however I am open to suggestions of alternative approaches.

I have tried using "rstanarm" i.e. stan_gamm4, but cannot find an appropriate way to introduce this prior data into the prior_smooth argument.