Dear R users, I have recently started to use the nlme and lme4 packages for fitting mixed effects models. Previously I used Genstat for doing this. Current I am struggling with analysing an experiment while modelling the correlation in the residual term in two directions, along rows and down columns. The data is below. It is an agricultural experiment with variety as a fixed effect and yield as the response. The experiment was laid out as a lattice square design but I wish to ignore this blocking structure and model the spatial correlation in the residuals. The layout is 10 rows by 15 columns of plots, each planted with a different variety. I wish to model the correlation with an AR(1) model. I have below example code based on the nlme package. I can model correlation in 1 direction and get the same results as I do with Genstat. slate.ar1 <- gls(Yield~as.factor(Variety),corr = corAR1(form = ~ 1|FieldRow),data=Slate) summary(slate.ar1) #Generalized least squares fit by REML # Model: Yield ~ as.factor(Variety) # Data: Slate # AIC BIC logLik # 552.6045 628.969 -249.3022 # #Correlation Structure: AR(1) # Formula: ~1 | FieldRow # Parameter estimate(s): # Phi #0.7325946 But to model spatial correlation in 2 directions along rows (FieldRow) and down columns (Field Column) I get an error message. slate.ar1*ar1 <- gls(Yield~as.factor(Variety),corr = corAR1(form = ~ 1|(FieldRow+FieldColumn),data=Slate)) #Error in corAR1(form = ~1 | (FieldRow + FieldColumn), data = Slate) : # unused argument(s) (data = Slate) Thank you for any advice you can give. Tom Data: Variety Yield FieldRow FieldColumn 1 10.03 1 1 2 13.56 1 2 4 14.12 1 3 3 12.39 1 4 5 15.08 1 5 19 19.67 1 6 23 15.72 1 7 2 19.69 1 8 6 17.47 1 9 15 15.98 1 10 18 16.3 1 11 25 16.33 1 12 9 12.55 1 13 11 12.77 1 14 2 15.72 1 15 6 15.31 2 1 7 15.4 2 2 9 12.5 2 3 8 16.58 2 4 10 11.85 2 5 8 16.05 2 6 12 15.5 2 7 16 15 2 8 25 16.42 2 9 4 15.04 2 10 5 16.8 2 11 7 15.26 2 12 16 14.52 2 13 23 14.8 2 14 14 14.82 2 15 21 11.26 3 1 22 14 3 2 24 13.29 3 3 23 12.87 3 4 25 15.55 3 5 11 13.95 3 6 20 16.96 3 7 24 15.7 3 8 3 14.04 3 9 7 12.85 3 10 6 14.73 3 11 13 17.61 3 12 22 16.95 3 13 4 13.64 3 14 20 17.9 3 15 11 12.61 4 1 12 14.23 4 2 14 11.1 4 3 13 17.35 4 4 15 16.17 4 5 22 18.2 4 6 1 13.51 4 7 10 12.97 4 8 14 14.12 4 9 18 15.06 4 10 24 15.12 4 11 1 13.55 4 12 15 15.24 4 13 17 14.78 4 14 8 13.71 4 15 16 14.58 5 1 17 20.36 5 2 19 21.19 5 3 18 19.12 5 4 20 18.93 5 5 5 17.48 5 6 9 14.5 5 7 13 17.4 5 8 17 14.5 5 9 21 15.23 5 10 12 13.64 5 11 19 16.9 5 12 3 13.34 5 13 10 12.39 5 14 21 15.57 5 15 3 16.23 6 1 18 18.62 6 2 8 16.45 6 3 13 18.88 6 4 23 15.27 6 5 16 16.06 6 6 24 18.42 6 7 10 11.86 6 8 13 14.62 6 9 2 12.42 6 10 10 10.82 6 11 4 13.04 6 12 17 12.67 6 13 11 12.66 6 14 23 12 6 15 1 13.31 7 1 16 14.17 7 2 6 16.11 7 3 11 14.54 7 4 21 17.9 7 5 12 17.67 7 6 20 19.17 7 7 1 12.64 7 8 9 10.6 7 9 23 9.51 7 10 12 11.3 7 11 6 12.66 7 12 24 12.89 7 13 18 12.6 7 14 5 11.74 7 15 5 12.11 8 1 20 14.11 8 2 10 11.83 8 3 15 15.5 8 4 25 16.6 8 5 4 15.26 8 6 7 16.81 8 7 18 15.45 8 8 21 12.9 8 9 15 9.76 8 10 19 12.4 8 11 13 11.81 8 12 1 9.17 8 13 25 12.87 8 14 7 9.75 8 15 2 13.88 9 1 17 14.53 9 2 7 13.84 9 3 12 16.69 9 4 22 17.38 9 5 25 18.45 9 6 3 17 9 7 14 15.28 9 8 17 13.73 9 9 6 12.4 9 10 21 12.52 9 11 20 15.91 9 12 8 14.28 9 13 2 15.09 9 14 14 12.73 9 15 4 14.43 10 1 19 16.67 10 2 9 15.49 10 3 14 14.59 10 4 24 17.22 10 5 8 15.83 10 6 11 14.9 10 7 22 16.07 10 8 5 13.15 10 9 19 11.74 10 10 3 14.43 10 11 22 16.49 10 12 15 14.07 10 13 9 13.15 10 14 16 13.18 10 15
Fitting a fixed effect model while modelling correlation in the residuals in 2 directions
3 messages · Thomas Bishop, Adam D. I. Kramer, Ben Bolker
On Tue, 15 Sep 2009, Thomas Bishop wrote:
slate.ar1*ar1 <- gls(Yield~as.factor(Variety),corr = corAR1(form = ~ 1|(FieldRow+FieldColumn),data=Slate))
...data=Slate is one paren too deep. Should be ...+FieldColumn)),data=Slate) :) --Adam
From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] On Behalf Of Adam D. I. Kramer [adik at ilovebacon.org]
Sent: Tuesday, September 15, 2009 12:22 AM
To: Thomas Bishop
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Fitting a fixed effect model while modelling correlation in the residuals in 2 directions
Sent: Tuesday, September 15, 2009 12:22 AM
To: Thomas Bishop
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Fitting a fixed effect model while modelling correlation in the residuals in 2 directions
On Tue, 15 Sep 2009, Thomas Bishop wrote: > slate.ar1*ar1 <- gls(Yield~as.factor(Variety),corr = corAR1(form = ~ > 1|(FieldRow+FieldColumn),data=Slate)) ...data=Slate is one paren too deep. Should be ...+FieldColumn)),data=Slate) Does corAR1 really work for spatial problems? I would have thought the definition would have been slightly different ...