1
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ This SCREAMS Bayes to me, especially since s() is available in brms. Maybe someone can come along and demonstrate this, but if now I will return after some meetings to do it. $\endgroup$ – Demetri Pananos Jun 16 '20 at 12:48
  • $\begingroup$ It would be much appreciated. I can see you are a biostatistician and this is actually a statistical genetics problem - locating mutational hotspots. The data are empirical case-control data and the priors are previously published data. $\endgroup$ – Adam Waring Jun 16 '20 at 15:14
  • 1
    $\begingroup$ Adam, what does pos mean? Does it mean "position"? Are you trying to determine with your model how "position" affects the chance of a study subject being a case or a control? The reason I ask this is because your current gam model, the way it is formulated, does not acknowledge that your outcome variable, cohort, is a binary outcome variable taking only two possible values: 0 (control) and 1 (case). Your gam model would need to use the option family = binomial(link="logit") the way it is currently stated to exploit the binary nature of the outcome variable. $\endgroup$ – Isabella Ghement Jun 16 '20 at 15:26
  • 1
    $\begingroup$ Adam, are you sure you don't want to use cohort status as a predictor of pos (e.g., position)? It's hard to comment more on this question, because there is too much information missing - what is the purpose of the study? What does pos mean? $\endgroup$ – Isabella Ghement Jun 16 '20 at 15:32
  • 1
    $\begingroup$ @AdamWaring I agree with Isabella, I think we need some more information. Let's start here: What are you trying to study and what variables are you using? $\endgroup$ – Demetri Pananos Jun 16 '20 at 15:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.