Hi,
I need to fit a model that tests nongenetic predictors of interest while accounting for the relatedness between individuals. For example, in the simplified simulated example below, I have a response variable, Y, and a fixed effect, X, measured on 20 subjects:
N < 20
Subj <factor(1:N)
X < rnorm(N, 0.2, 4)
Y < rnorm(N, X*0.4, 6)
Data < data.frame(Subj, X, Y)
I also have a hypothetical relatedness matrix (Rmat):
Rmat < diag(N)
corrs < runif(N*(N1)/2, 0, 1)
Rmat[lower.tri(Rmat)] <corrs
Rmat < t(Rmat)
Rmat[lower.tri(Rmat)] <corrs
rownames(Rmat)<Subj; colnames(Rmat)<Subj
Can I use GenABEL to fit a mixed model that tests for the association between X and Y while accounting for the relatedness/ kinship matrix? I seem to be getting close using the polygenic function, but this doesn't provide standard errors/pvalues for the fixed effect(s), while if I try to run grammar it seems to want genetic/SNP data as my predictors of interest.
Any help/ pointers greatly appreciated!
Many thanks,
Nick
Testing nongenetic predictors in polygenic mixed model
Forum rules
Please remember not to post any sensitive data on this public forum.
The first few posts of newly registered users will be moderated in order to filter out any spammers.
When get a solution to the problem you posted, please change the topic name (e.g. from "how to ..." to "[SOLVED] how to ..."). This will make it easier for the community to follow the posts yet to be attended.
Please remember not to post any sensitive data on this public forum.
The first few posts of newly registered users will be moderated in order to filter out any spammers.
When get a solution to the problem you posted, please change the topic name (e.g. from "how to ..." to "[SOLVED] how to ..."). This will make it easier for the community to follow the posts yet to be attended.

 GenABEL senior expert
 Posts: 151
 Joined: Wed Feb 09, 2011 3:24 pm
Re: Testing nongenetic predictors in polygenic mixed model
HI,
you have basically three choiches:
1) use polygenic in combiantion with the likelihood ratio test, but you will not get standard errors and if you add more than one variable at a time you will not get individual variables pvalues. In this case you should first fit a null model like:
Then you fit the full model:
then you can use the likelyhood ratio test.
Beta estimates are in yourmod$h2an$estimate. They do not present names but they are ordered in the order you put them. The vector is 2 elements longer since heritability and residual variance are there as well. Remember that if you want this to make sense you have to use exactly the same people in the 2 models (model.frame() is great to achieve this). Also if you add more than 1 variable in the test model, remember to add degrees of freedom when you estimate the pvalue.
2) use polygenichglm which will give you the effects and the standard errors and everything else, but is a hierarchical mixed model which is not exactly the same.
3) use GRAMMAR+, in this case you first regress out kinship, you get the enviromental corrected residuals ($grresidualY) and then use this to run a normal regression with glm.
Hope this solves your problems.
Best.
Nicola
you have basically three choiches:
1) use polygenic in combiantion with the likelihood ratio test, but you will not get standard errors and if you add more than one variable at a time you will not get individual variables pvalues. In this case you should first fit a null model like:
Code: Select all
nullmod=polygenic(y~sex+age,data=yourdata, kinship=yourkinship)
Then you fit the full model:
Code: Select all
yourmod=polygenic(y~sex+age+ yourtest.variable,data=yourdata, kinship=yourkinship)
then you can use the likelyhood ratio test.
Code: Select all
p=pchisq(nullmod$h2an$minimum  yourmod$h2an$minimum, df=1, lower=F)
Beta estimates are in yourmod$h2an$estimate. They do not present names but they are ordered in the order you put them. The vector is 2 elements longer since heritability and residual variance are there as well. Remember that if you want this to make sense you have to use exactly the same people in the 2 models (model.frame() is great to achieve this). Also if you add more than 1 variable in the test model, remember to add degrees of freedom when you estimate the pvalue.
2) use polygenichglm which will give you the effects and the standard errors and everything else, but is a hierarchical mixed model which is not exactly the same.
3) use GRAMMAR+, in this case you first regress out kinship, you get the enviromental corrected residuals ($grresidualY) and then use this to run a normal regression with glm.
Hope this solves your problems.
Best.
Nicola
Re: Testing nongenetic predictors in polygenic mixed model
Hi Nicola,
Thanks very much for your response. The first two solutions you offered seem to address my problem perfectly.
For your third proposal, if I were to fit the "enviromental corrected residuals" in a regular regression, wouldn't this assume that each individual is unrelated (i.e. independent)? I can see how regressing out the kinship would account for any confounding (such as that due to population structure), but am struggling to understand how any relatedness/ nonindependence of individuals would be addressed using this approach.
Best regards,
Nick
Thanks very much for your response. The first two solutions you offered seem to address my problem perfectly.
For your third proposal, if I were to fit the "enviromental corrected residuals" in a regular regression, wouldn't this assume that each individual is unrelated (i.e. independent)? I can see how regressing out the kinship would account for any confounding (such as that due to population structure), but am struggling to understand how any relatedness/ nonindependence of individuals would be addressed using this approach.
Best regards,
Nick
Who is online
Users browsing this forum: No registered users and 1 guest