Dear Jarrod,

I'm very interested in this SIR feature of MCMCglmm, which I hadn't been aware of until your e-mail earlier today. Do I understand correctly that this should allow for feedback effects between units, reflecting the fact that y(i) and y(not-i) affect each other? If so, this could be a useful way of fitting, for example, spatial autoregressive models, where some units may affect each other by virtue of being neighbours, and we would like to know the magnitude of these effects. (You mentioned time series, but some web searches haven't turned up much on MCMCglmm and time series. And time only flows in one direction, whereas units in space can have two-way effects.)

However, in looking at the documentation and code, I've been struggling a bit with the "sir(~XX, ~XX)" part. Is there a way to specify a "connectivity" matrix directly, rather than using the "sir" function? Looking at the MCMCglmm code, as I understand it, "sir" is both a way of generating a connectivity matrix AND a way of telling MCMCglmm to treat the sir(~XX, ~XX) differently than other covariates.

In the case of spatial autoregressive models, each element on the main diagonal of the connectivity matrix is 0, and rows are often standardised to sum to 1. Such a connectivity matrix might, for example, look like:

N <- 100
W <- matrix(rbinom(N^2, 1, 0.1), N, N)*upper.tri(matrix(1, N, N)) # chance of being neighbours is 0.1
W <- W + t(W) # matrix is symmetrical, with zeros on the main diagonal

But I'm having a hard time figuring out how to use "sir" to specify such a matrix in the middle of a call to MCMCglmm.

If you see what I mean, can you offer any suggestions? I might take a stab at modifying MCMCglmm to make this possible, but before I do that (and I'm not certain I'll have the know-how) I wanted to check whether you had any ideas.

Many thanks for any assistance or clarification you can provide.
- Malcolm



> Message: 5
> Date: Wed, 1 Dec 2010 10:28:43 +0000
> From: Jarrod Hadfield <j.hadfield@ed.ac.uk>
> To: Szymek Drobniak <geralttee@gmail.com>
> Cc: r-sig-mixed-models@r-project.org
> Subject: Re: [R-sig-ME] Another MCMCglmm question - fixing
> 	correlations
> Message-ID: <0AF131D6-557C-4324-A65A-A175E41DBA8F@ed.ac.uk>
> Content-Type: text/plain; charset=US-ASCII; format=flowed; delsp=yes
> 
> 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))


	[[alternative HTML version deleted]]


