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

[Jarrod Hadfield]

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  
license. If not ......

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  
some additional residual variation:

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.