# [R-sig-ME] Another MCMCglmm question - fixing correlations

Wed Dec 1 11:28:43 CET 2010

Hi,

There is a way, but its quite involved and you may find it easier to
fit this type of model in something like ASReml if you can get a

I have implemented simultaneous-recursive (SIR) mixed models  in
MCMCglmm, although currently they are not well tested or documented
(See Chapter 9 of the CourseNotes) and cannot be used with non-
Gaussian data or certain patterns of missing data. Nevertheless, in
this example the existing code should work fine.  A simple SIR model
has the form:

y_{i} = mu + \lambda*y_{j} + u_{i} + e_{i}

where mu is the intercept, u is a random effect (e.g. breeding value
in your case) and e is a residual. The new part of the model is the
structural parameter \lambda multiplied by y_{j}, the jth element of
the response vector, which can be useful for modelling the effects of
behavioural interactions or time series etc.  However, placing y's on
the LHS and RHS complicates the likelihood, but MCMCglmm deals with
this.

If we set i=j then the above equation can be rearranged:

y_{i}  - \lambda*y_{i} = mu + u_{i} + e_{i}
y_{i}(1  - \lambda) = mu + u_{i} + e_{i}
y_{i} = (mu + u_{i} + e_{i})/(1  - \lambda)

from this it can be seen that mu, u, e, and \lambda cannot be uniquely
estimated without restrictions.  To illustrate we'll simulate some data:

fac<-gl(25,4)
y<-rnorm(25)[fac]+rnorm(100, -1, sqrt(2))

# mu = -1, VAR(u)=1, VAR(e) = 2, \lambda = 0

dat<-data.frame(y=y, fac=fac)

prior1 = list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=0.002)))

m1<-MCMCglmm(y~1, random=~fac,  data=dat, prior=prior1)

# fitting a non-SIR model gives estimates consistent with the true
values, which you can verify through summary(m1). Now for the SIR
model. Because we know all parameters cannot be uniquely estimated
I've set the residual variance to one arbitrarily.

prior2 = list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=0.002)))

m2<-MCMCglmm(y~1+sir(~units, ~units), random=~fac,  data=dat,
prior=prior2)   # sir(~units, ~units) sets i=j.

# The estimates do not look consistent with the real values.  However,
we can rescale them by 1-\lambda and the estimates  are close to being
identical up to Monte Carlo error (the prior will have slightly
different effects)

HPDinterval(m1$VCV[,1]) HPDinterval(m2$VCV[,1]/(1-m2$Lambda)^2) HPDinterval(m1$VCV[,2])
HPDinterval(m2$VCV[,2]/(1-m2$Lambda)^2)

HPDinterval(m1$Sol) HPDinterval(m2$Sol/(1-m2$Lambda)) # And what's the point? Lets Imagine the same data are associated with 2 genders: dat$sex<-gl(2,1,100, labels=c("Female", "Male"))

# and that Females are twice as big as males, and on top of this have

dat$y[which(dat$sex=="Female")]<-2*dat$y[which(dat$sex=="Female")]
dat$y[which(dat$sex=="Female")]<-dat$y[which(dat$sex=="Female")]
+rnorm(50)

For Males   mu = -1, VAR(u)=1, VAR(e) = 2 as before
For Females   mu = -2, VAR(u)=4, VAR(e) = 8+1=9.

# Importantly the correlation between u's in the different sexes is 1,
despite them having the same variance. By only fitting a sir model to
females

prior3=list(R=list(V=diag(2), nu=0.002), G=list(G1=list(V=1, nu=0.002)))

m3<-MCMCglmm(y~1+sir(~units:at.level(sex,"Female"),
~units:at.level(sex,"Female")), random=~fac, rcov=~idh(sex):units,
data=dat, prior=prior3)

# we can extract the male estimates and they are consistent with what
we expect.

mean(m3$VCV[,1]) # var(u) males mean(m3$VCV[,3])  # var(e) males
mean(m3$Sol) # mu males # we can then rescale the female estimates as before: mean(m3$VCV[,1]/(1-m3$Lambda)^2) # var(u) females mean(m3$VCV[,2]/(1-m3$Lambda)^2) # var(e) females mean(m3$Sol/(1-m3\$Lambda))  # mu  females

# and the estimates should be good.

Note in this last model we did not need to place restrictions on the
variance components through the prior. Restrictions are required, but
we achieved this by having ~fac rather than something  like
us(sex):fac or idh(sex):fac.  More complicated models could be
generated which place other types of restriction on the covariance
matrices.

Cheers,

Jarrod

On 30 Nov 2010, at 22:19, Szymek Drobniak wrote:

> Hi,
>
> is it possible in MCMglmm to fix the (co)variance matrix so that the
> resulting variances could vary, but covariances had always the value
> ensuring that r=1? It would be useful in testing cross-sex rG where
> the H0 should be rG=1 rather than rG=0 (as defined in idh()). I looked
> through available (co)variance structures but none seems suitable.
>
> best,
> sz.
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.