[R-sig-ME] SIR in MCMCglmm - syntax question

[Jarrod Hadfield]

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>