[R-sig-ME] SIR in MCMCglmm - syntax question
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Apr 6 20:21:16 CEST 2011
Hi,
Code is below. Be VERY careful with SIR models generally (working out
whether the structural parameters are idenitifiable is not an easy
task) and in MCMCglmm in particular (they are not yet thoroughly
tested or documented). With respect to the first issue, I have set
the residual covariance to zero in the example - my gut feeling is
that if you estimate the residual covariance too then things start to
become non-identifiable. However, I'm far from sure.
Cheers,
Jarrod
library(MASS)
library(MCMCglmm)
G<-matrix(c(1,0.5,0.5,1),2,2)
E<-diag(2)
# generate covariance matrices for random effects (G) and residuals (E)
u<-mvrnorm(100, c(0,0), G)
e<-mvrnorm(400, c(0,0), E)
# sample random effects and residuals
beta1<-rnorm(3)
# draw the three fixed effects for trait 1: intercept, age and smoking
beta2<-rnorm(4)
# draw the four fixed effects for trait 2: intercept, age, drinking
and exercise
lambda<-runif(2, -1,1)
# draw the two structural parameters
dat<-data.frame(smoking=rbinom(400,1,0.5), age=runif(400),
exercise=rnorm(400), drinking=rnorm(400), group=gl(100,4))
# generate covariates and group (random effect) identifiers
X1<-model.matrix(~age+smoking, dat)
# design matrix for fixed effects on trait 1
X2<-model.matrix(~age+drinking+exercise, dat)
# design matrix for fixed effects on trait 2
Ly<-cbind(X1%*%beta1, X2%*%beta2)+u[dat$group,]+e
# generate data with out feed back
L<-diag(2)
L[1,2]<--lambda[1]
L[2,1]<--lambda[2]
# generate mini-Lambda (i.e. 2x2 rather than 800x800)
y<-t(apply(Ly, 1, function(x){solve(L,x)}))
dat$y1=y[,1]
dat$y2=y[,2]
# solve for y for each set of observations
m1<-MCMCglmm(cbind(y1, y2)~trait-1+at.level(trait,1):age
+at.level(trait,1):smoking+at.level(trait,2):age+at.level(trait,
2):drinking+at.level(trait,2):exercise+sir(~at.level(trait,1):units,
~at.level(trait,2):units)+sir(~at.level(trait,2):units,
~at.level(trait,1):units), random=~us(trait):group,
rcov=~idh(trait):units, data=dat, family=rep("gaussian", 2))
# fit model
tpar<-c(beta1[1], beta2[1], beta1[2:3], beta2[2:4],lambda)
devAskNewPage.orig <- devAskNewPage()
devAskNewPage(TRUE)
for(i in 1:9){
hist(cbind(m1$Sol, m1$Lambda)[,i], main=colnames(cbind(m1$Sol,
m1$Lambda))[i], breaks=30)
abline(v=tpar[i], col="red")
dev.interactive()
}
# compare posterior with true values
devAskNewPage(devAskNewPage.orig)
On 6 Apr 2011, at 14:25, Szymek Drobniak wrote:
> Dear all,
>
> I know that SIR models are not yet well documented or widely used but
> I've started thinking about using them in my research - and ran into
> some interpretation problems. As far as I understand from
> Gianola&Sorensen 2004 SIR models should allow for estimation of Lambda
> elements. However I find it difficult to translated this Lambda matrix
> into the S matrix as generated by sir() in MCMCglmm. Have any of you
> used these models? Or otherwise - how should be the correct way of
> specyfying a SIR model for the example in Gianola&Sorensen about
> smoking and drinking (pp. 1418-1419). That would probably enlighten me
> ;)
>
> Cheers
> sz.
>
> --
> Szymon Drobniak || Population Ecology Group
> Institute of Environmental Sciences, Jagiellonian University
> ul. Gronostajowa 7, 30-387 Kraków, POLAND
> tel.: +48 12 664 52 19 fax: +48 12 664 69 12
>
> www.eko.uj.edu.pl/drobniak
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110406/7c86066c/attachment-0004.pl>
More information about the R-sig-mixed-models
mailing list