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

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
L[2,1]<--lambda

# 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, beta2, beta1[2:3], beta2[2:4],lambda)

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

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>
```