Hi, I am trying to parameterize the following mixed model (following Piepho and Ogutu 2002), to test for a trend over time, using multiple sites: y[ij]=mu+b[j]+a[i]+w[j]*(beta +t[i])+c[ij] where: y[ij]= a response variable at site i and year j mu = fixed intercept Beta=fixed slope w[j]=constant representing the jth year (covariate) b[j]=random effect of jth year, iid N(0,sigma2[b]) a[i]=random effect of the ith site, iid N(0, sigma2[a]) t[i]=random effect of ith site, iid N(0, sigma2[t]) c[ij]=random error associated with ith site and jth year I would like to assume that an unstructured relationship applies to a[i] and t[i] (i.e., I would like to assume that the random effects a[i] and t[i] are drawn from a multivariate normal distribution with non-zero covariance parameter). These random effects are assumed to be independent from the b[j]'s and from the c[ij]'s. I have tried several approaches, but cannot seem to duplicate the results presented in Piepho and Ogutu using R's lme function (but I can reproduce the results using SAS proc mixed). In SAS, the model is fit using: proc mixed method=REML nobound; class year site; model y=w site/ddfm=satterth s; random int/sub=year; random int w/sub=site type=un; run; Any help would be greatly appreciated! Reference: Piepho, H-P. and J.O.Ogutu. 2002. A simple mixed model for trend analysis in wildlife populations. Journal of Agricultural, Biological, and Environmental Statistics, 7(3):350-360. Thanks, John
Hi, try something like: lme(y~w,random=list(~1|year,~1+w|site)) Christian ----- Original Message ----- From: "John Fieberg" <John.Fieberg at dnr.state.mn.us> To: <R-help at stat.math.ethz.ch> Sent: Wednesday, April 02, 2003 10:36 PM Subject: [R] lme parameterization question> Hi, > > I am trying to parameterize the following mixed model (following Piepho > and Ogutu 2002), to test for a trend over time, using multiple sites: > > y[ij]=mu+b[j]+a[i]+w[j]*(beta +t[i])+c[ij] > > where: > y[ij]= a response variable at site i and year j > mu = fixed intercept > Beta=fixed slope > w[j]=constant representing the jth year (covariate) > b[j]=random effect of jth year, iid N(0,sigma2[b]) > a[i]=random effect of the ith site, iid N(0, sigma2[a]) > t[i]=random effect of ith site, iid N(0, sigma2[t]) > c[ij]=random error associated with ith site and jth year
Have you looked at Pinheiro and Bates (2000) Mixed Effects Models in S and s-Plus? "lme" is great, but I couldn't make it work until I spent some time with that book. Hope this helps. Spencer Graves John Fieberg wrote:> Hi, > > I am trying to parameterize the following mixed model (following Piepho > and Ogutu 2002), to test for a trend over time, using multiple sites: > > y[ij]=mu+b[j]+a[i]+w[j]*(beta +t[i])+c[ij] > > where: > y[ij]= a response variable at site i and year j > mu = fixed intercept > Beta=fixed slope > w[j]=constant representing the jth year (covariate) > b[j]=random effect of jth year, iid N(0,sigma2[b]) > a[i]=random effect of the ith site, iid N(0, sigma2[a]) > t[i]=random effect of ith site, iid N(0, sigma2[t]) > c[ij]=random error associated with ith site and jth year > > I would like to assume that an unstructured relationship applies to > a[i] and t[i] (i.e., I would like to assume that the random effects a[i] > and t[i] are drawn from a multivariate normal distribution with non-zero > covariance parameter). These random effects are assumed to be > independent from the b[j]'s and from the c[ij]'s. I have tried several > approaches, but cannot seem to duplicate the results presented in Piepho > and Ogutu using R's lme function (but I can reproduce the results using > SAS proc mixed). > > In SAS, the model is fit using: > > proc mixed method=REML nobound; > class year site; > model y=w site/ddfm=satterth s; > random int/sub=year; > random int w/sub=site type=un; > run; > > Any help would be greatly appreciated! > > Reference: > Piepho, H-P. and J.O.Ogutu. 2002. A simple mixed model for trend > analysis in wildlife populations. Journal of Agricultural, Biological, > and Environmental Statistics, 7(3):350-360. > > > Thanks, > > John > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-help