0
$\begingroup$

snp1, snp2, and snp3 are "versions" of geneA. I want to determine if geneA significantly effects disease. I have sex as a potentially confounding factor. If geneA affects disease I want to do pairwise comparisons to see which snp affects disease most significantly. I have tens of thousands of these genes to look at, so I am trying to limit type 1 error by only looking at pairwise comparisons in genes that have significant overall effect.

This is my dummy data for geneA

data <- data.frame("snp1"=c(runif(n=150, min=0,max=2),
                          c(runif(n=50, min=0,max=2))),
                  "snp2"=c(runif(n=50, min=0, max=.2),
                           runif(n=50, min=0, max=.2),
                           runif(n=50, min=1.5, max=2),
                           runif(n=50, min=1.5, max=2)),
                  "snp3"=c(runif(n=50, min=0, max=.2),
                          runif(n=50, min=0, max=.2),
                          runif(n=50, min=1.5, max=2),
                          runif(n=50, min=1.5, max=2)),
                  "sex"=runif(n=50, min=0, max=1),
                   "disease"=c(rbinom(150, 1, 0.1),
                               rbinom(50, 1, 0.9)))

This is my likelihood ratio test approach for the whole gene, with snp1 left out as reference:

multi_snp_full <- glm(disease ~ snp2 + snp3 + sex, data=data, family="binomial")
null <- glm(disease ~ sex, data=data, family="binomial")
anova( null, multi_snp_full , test="Chisq")

If I wanted to go back and see how snp2 affect disease would it just be this? (eg no LR test):

single_snp_test <- glm(disease ~ snp2 + sex, data=data, family="binomial")
summary(single_snp_test )
Improve this question
$\endgroup$
4
  • $\begingroup$ Since you have snp3 in your full model as well, you should look at the effect of snp2 in multi_snp_full $\endgroup$
    – Demetri Pananos
    Commented Nov 20, 2020 at 16:29
  • $\begingroup$ >you should look at the effect of snp2 in multi_snp_full I get that snp3 is in the full model, but I don't really get what you mean by the second half of you comment. Could you please clarify what you mean? $\endgroup$
    – curious
    Commented Nov 20, 2020 at 16:37
  • $\begingroup$ Wait are you saying that if I was intersted in determining the effect of snp2 or snp3 by themselves I would just go: summary(multi_snp_full) and look at the effect and p value for snp2 or snp3 in the readout? So no need to make another model like I did with single_snp_test $\endgroup$
    – curious
    Commented Nov 20, 2020 at 16:47
  • $\begingroup$ Yes, that is correct. $\endgroup$
    – Demetri Pananos
    Commented Nov 20, 2020 at 17:43

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.