SIR in MCMCglmm - syntax question

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)

SIR in MCMCglmm - syntax question

Thread (2 messages)