# lme4 or other open source R package code equivalent to asreml-R

I want to fit mixed model using lme4, nlme, baysian regression package or any available.

Mixed model in Asreml- R coding conventions

before going into specifics, we might want to have details on asreml-R conventions, for those who are unfamiliar with ASREML codes.

y = Xτ + Zu + e ........................(1) ;


the usual mixed model with, y denotes the n × 1 vector of observations,where τ is the p×1 vector of ﬁxed eﬀects, X is an n×p design matrix of full column rank which associates observations with the appropriate combination of ﬁxed eﬀects, u is the q × 1 vector of random eﬀects, Z is the n × q design matrix which associates observations with the appropriate combination of random eﬀects, and e is the n × 1 vector of residual errors.The model (1) is called a linear mixed model or linear mixed eﬀects model. It is assumed

where the matrices G and R are functions of parameters γ and φ, respectively.

The parameter θ is a variance parameter which we will refer to as the scale parameter.

In mixed eﬀects models with more than one residual variance, arising for example in the analysis of data with more than one section or variate, the parameter θ is ﬁxed to one. In mixed eﬀects models with a single residual variance then θ is equal to the residual variance (σ2). In this case R must be correlation matrix. Further details on the models are provided in the Asreml manual (link).

Variance structures for the errors: R structure and Variance structures for the random eﬀects: G structures can be specified.

variance modelling in asreml() it is important to understand the formation of variance structures via direct products. The usual least squares assumption (and the default in asreml()) is that these are independently and identically distributed (IID). However, if the data was from a field experiment laid out in a rectangular array of r rows by c columns, say, we could arrange the residuals e as a matrix and potentially consider that they were autocorrelated within rows and columns.Writing the residuals as a vector in field order, that is, by sorting the residuals rows within columns (plots within blocks) the variance of the residuals might then be

are correlation matrices for the row model (order r, autocorrelation parameter ½r) and column model (order c, autocorrelation parameter ½c) respectively. More specifically, a two-dimensional separable autoregressive spatial structure (AR1 x ­ AR1) is sometimes assumed for the common errors in a field trial analysis.

The example data:

nin89 is from asreml-R library, where different varities were grown in replications / blocks in rectangular field. To control additional variability in row or column direction each plot is referenced as Row and Column variables (row column design). Thus this row column design with blocking. Yield is measured variable.

Example models

I need something equivalent to the asreml-R codes:

The simple model syntax will look like the follows:

 rcb.asr <- asreml(yield ∼ Variety, random = ∼ Replicate, data = nin89)
.....model 0


The linear model is specified in the fixed (required), random (optional) and rcov (error component) arguments as formula objects.The default is a simple error term and does not need to be formally specified for error term as in the model 0.

here the variety is fixed effect and random is replicates (blocks). Beside random and fixed terms we can specify error term. Which is default in this model 0. The residual or error component of the model is specified in a formula object through the rcov argument, see the following models 1:4.

The following model1 is more complex in which both G (random) and R (error) structure are specified.

Model 1:

data(nin89)

# Model 1: RCB analysis with G and R structure
rcb.asr <- asreml(yield ~ Variety, random = ~ idv(Replicate),
rcov = ~ idv(units), data = nin89)


This model is equivalent to above model 0, and introduces the use of G and R variance model. Here the option random and rcov specifies random and rcov formulae to explicitly specify the G and R structures. where idv() is the special model function in asreml() that identifies the variance model. The expression idv(units) explicitly sets the variance matrix for e to a scaled identity.

# Model 2: two-dimensional spatial model with correlation in one direction

  sp.asr <- asreml(yield ~ Variety, rcov = ~ Column:ar1(Row), data = nin89)


experimental units of nin89 are indexed by Column and Row. So we expect random variation in two direction - row and column direction in this case. where ar1() is a special function specifying a first order autoregressive variance model for Row. This call specifies a two-dimensional spatial structure for error but with spatial correlation in the row direction only.The variance model for Column is identity (id()) but does not need to be formally specified as this is the default.

# model 3: two-dimensional spatial model, error structure in both direction

 sp.asr <- asreml(yield ~ Variety, rcov = ~ ar1(Column):ar1(Row),
data = nin89)
sp.asr <- asreml(yield ~ Variety, random = ~ units,
rcov = ~ ar1(Column):ar1(Row), data = nin89)


similar to above model 2, however the correlation is two direction - autoregressive one.

I am not sure how much of these models are possible with open source R packages. Even if solution of any one of these models will be of great help. Even if the bouty of +50 can stimulate to develop such package will be of great help !

See MAYSaseen has provided output from each model and data (as answer) for comparision.

Edits: The following is suggestion I received in mixed model discussion forum: " You might look at the regress and spatialCovariance packages of David Clifford. The former allows fitting of (Gaussian) mixed models where you can specify the structure of the covariance matrix very flexibly (for example, I have used it for pedigree data). The spatialCovariance package uses regress to provide more elaborate models than AR1xAR1, but may be applicable. You may have to correspond with the author about applying it to your exact problem."

## migrated from stackoverflow.comNov 21 '11 at 15:15

This question came from our site for professional and enthusiast programmers.

• I'm pretty sure that models 2-4 are not possible in lme4. Can you (a) tell us why you need to do this in lme4 rather than asreml-R (b) consider posting on r-sig-mixed-models where there is more relevant expertise? – Ben Bolker Nov 1 '11 at 3:16
• basic idea is asreml-R require a license (at least for developed country users), if it is possible in lme4 or other mixed model packages that would be great... – John Nov 1 '11 at 3:52
• I think this is not going to be easy. I think your best bet might be to define a new corStruct in nlme (for anisotropic correlations) ... It would help if you could briefly state (in words or equations) the statistical models corresponding to these ASREML statements, since we are not all familiar with ASREML syntax ... – Ben Bolker Nov 1 '11 at 15:47
• The following is comments in mixed model group: You might look at the regress and spatialCovariance packages of David Clifford. The former allows fitting of (Gaussian) mixed models where you can specify the structure of the covariance matrix very flexibly (for example, I have used it for pedigree data). The spatialCovariance package uses regress to provide more elaborate models than AR1xAR1, but may be applicable. You may have to correspond with the author about applying it to your exact problem. – John Nov 17 '11 at 23:40
• if I get a chance I will try to tackle as much of this as I can, but quite frankly I may not get to it, I've got a lot on my plate. Looking into the packages that David Clifford suggested sounds like a great idea -- maybe you can solve your own problem that way ... I'm pretty sure that model 1 can be done with MCMCglmm, and I'm pretty sure that (other than the spatialCovariance mentioned, which I'm unfamiliar with) the only way to get it done in R is by defining new corStructs -- which is possible, but not trivial. – Ben Bolker Nov 18 '11 at 0:16

You can fit this model with AD Model Builder. AD Model Builder is free software for building general nonlinear models including general nonlinear random effects models. So for example you could fit a negative binomial spatial model where both the mean and over dispersion had an ar(1) x ar(1) structure. I built the code for this example and fit it to the data. If anyone is interested it is probably better to discuss this on the list at http://admb-project.org

Note: There is an R version of ADMB, but the features available in the R package are a subset of the standalone ADMB software.

For this example it is easier to create an ASCII file with the data, read it into the ADMB program, run the program, and then read the parameter estimates etc back into R for whatever you want to do.

You should understand that ADMB is not a collection of packages, but rather a language for writing nonlinear parameter estimation software. As I said before it is better to discuss this on the ADMB list where everyone knows about the software. After it is done and you understand the model you can post the results here. However here is a link to the ML and REML codes I put together for the wheat data.

• Is there R interphase to connect with AD Model Builder? – John Nov 27 '11 at 9:26

Model 0

ASReml-R

rcb0.asr <- asreml(yield~Variety, random=~Rep, data=nin89, na.method.X="include")
summary(rcb0.asr)
$call asreml(fixed = yield ~ Variety, random = ~Rep, data = nin89, na.method.X = "include")$loglik
[1] -454.4691

$nedf [1] 168$sigma
[1] 7.041475

$varcomp gamma component std.error z.ratio constraint Rep!Rep.var 0.1993231 9.882911 8.792829 1.123974 Positive R!variance 1.0000000 49.582368 5.458839 9.082951 Positive attr(,"class") [1] "summary.asreml" summary(rcb0.asr)$varcomp
gamma component std.error  z.ratio constraint
Rep!Rep.var 0.1993231  9.882911  8.792829 1.123974   Positive
R!variance  1.0000000 49.582368  5.458839 9.082951   Positive

> anova(rcb0.asr)
Wald tests for fixed effects

Response: yield

Df Sum of Sq Wald statistic Pr(Chisq)
(Intercept)    1   12001.6        242.054    <2e-16 ***
Variety       55    2387.5         48.152    0.7317
residual (MS)         49.6
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> coef(rcb0.asr)$fixed effect Variety_ARAPAHOE 0.0000 Variety_BRULE -3.3625 Variety_BUCKSKIN -3.8750 Variety_CENTURA -7.7875 Variety_CENTURK78 0.8625 Variety_CHEYENNE -1.3750 Variety_CODY -8.2250 Variety_COLT -2.4375 Variety_GAGE -4.9250 Variety_HOMESTEAD -1.8000 Variety_KS831374 -5.3125 Variety_LANCER -0.8750 Variety_LANCOTA -2.8875 Variety_NE83404 -2.0500 Variety_NE83406 -5.1625 Variety_NE83407 -6.7500 Variety_NE83432 -9.7125 Variety_NE83498 0.6875 Variety_NE83T12 -7.8750 Variety_NE84557 -8.9125 Variety_NE85556 -3.0500 Variety_NE85623 -7.7125 Variety_NE86482 -5.1500 Variety_NE86501 1.5000 Variety_NE86503 3.2125 Variety_NE86507 -5.6500 Variety_NE86509 -2.5875 Variety_NE86527 -7.4250 Variety_NE86582 -4.9000 Variety_NE86606 0.3250 Variety_NE86607 -0.1125 Variety_NE86T666 -7.9000 Variety_NE87403 -4.3125 Variety_NE87408 -3.1375 Variety_NE87409 -8.0625 Variety_NE87446 -1.7625 Variety_NE87451 -4.8250 Variety_NE87457 -5.5250 Variety_NE87463 -3.5250 Variety_NE87499 -9.0250 Variety_NE87512 -6.1875 Variety_NE87513 -2.6250 Variety_NE87522 -4.4375 Variety_NE87612 -7.6375 Variety_NE87613 -0.0375 Variety_NE87615 -3.7500 Variety_NE87619 1.8250 Variety_NE87627 -6.2125 Variety_NORKAN -5.0250 Variety_REDLAND 1.0625 Variety_ROUGHRIDER -8.2500 Variety_SCOUT66 -1.9125 Variety_SIOUXLAND 0.6750 Variety_TAM107 -1.0375 Variety_TAM200 -8.2000 Variety_VONA -5.8375 (Intercept) 29.4375 > coef(rcb0.asr)$random
effect
Rep_1  1.8795997
Rep_2  2.8432659
Rep_3 -0.8712739
Rep_4 -3.8515918


lme4

> rcb0.lmer <- lmer(yield~Variety+(1|Rep), data=nin89)
> print(rcb0.lmer, corr=FALSE)
Linear mixed model fit by REML
Formula: yield ~ Variety + (1 | Rep)
Data: nin89
AIC  BIC logLik deviance REMLdev
1334 1532 -608.9     1456    1218
Random effects:
Groups   Name        Variance Std.Dev.
Rep      (Intercept)  9.8829  3.1437
Residual             49.5824  7.0415
Number of obs: 224, groups: Rep, 4

Fixed effects:
Estimate Std. Error t value
(Intercept)        29.4375     3.8556   7.635
VarietyBRULE       -3.3625     4.9791  -0.675
VarietyBUCKSKIN    -3.8750     4.9791  -0.778
VarietyCENTURA     -7.7875     4.9791  -1.564
VarietyCENTURK78    0.8625     4.9791   0.173
VarietyCHEYENNE    -1.3750     4.9791  -0.276
VarietyCODY        -8.2250     4.9791  -1.652
VarietyCOLT        -2.4375     4.9791  -0.490
VarietyGAGE        -4.9250     4.9791  -0.989
VarietyKS831374    -5.3125     4.9791  -1.067
VarietyLANCER      -0.8750     4.9791  -0.176
VarietyLANCOTA     -2.8875     4.9791  -0.580
VarietyNE83404     -2.0500     4.9791  -0.412
VarietyNE83406     -5.1625     4.9791  -1.037
VarietyNE83407     -6.7500     4.9791  -1.356
VarietyNE83432     -9.7125     4.9791  -1.951
VarietyNE83498      0.6875     4.9791   0.138
VarietyNE83T12     -7.8750     4.9791  -1.582
VarietyNE84557     -8.9125     4.9791  -1.790
VarietyNE85556     -3.0500     4.9791  -0.613
VarietyNE85623     -7.7125     4.9791  -1.549
VarietyNE86482     -5.1500     4.9791  -1.034
VarietyNE86501      1.5000     4.9791   0.301
VarietyNE86503      3.2125     4.9791   0.645
VarietyNE86507     -5.6500     4.9791  -1.135
VarietyNE86509     -2.5875     4.9791  -0.520
VarietyNE86527     -7.4250     4.9791  -1.491
VarietyNE86582     -4.9000     4.9791  -0.984
VarietyNE86606      0.3250     4.9791   0.065
VarietyNE86607     -0.1125     4.9791  -0.023
VarietyNE86T666    -7.9000     4.9791  -1.587
VarietyNE87403     -4.3125     4.9791  -0.866
VarietyNE87408     -3.1375     4.9791  -0.630
VarietyNE87409     -8.0625     4.9791  -1.619
VarietyNE87446     -1.7625     4.9791  -0.354
VarietyNE87451     -4.8250     4.9791  -0.969
VarietyNE87457     -5.5250     4.9791  -1.110
VarietyNE87463     -3.5250     4.9791  -0.708
VarietyNE87499     -9.0250     4.9791  -1.813
VarietyNE87512     -6.1875     4.9791  -1.243
VarietyNE87513     -2.6250     4.9791  -0.527
VarietyNE87522     -4.4375     4.9791  -0.891
VarietyNE87612     -7.6375     4.9791  -1.534
VarietyNE87613     -0.0375     4.9791  -0.008
VarietyNE87615     -3.7500     4.9791  -0.753
VarietyNE87619      1.8250     4.9791   0.367
VarietyNE87627     -6.2125     4.9791  -1.248
VarietyNORKAN      -5.0250     4.9791  -1.009
VarietyREDLAND      1.0625     4.9791   0.213
VarietyROUGHRIDER  -8.2500     4.9791  -1.657
VarietySCOUT66     -1.9125     4.9791  -0.384
VarietySIOUXLAND    0.6750     4.9791   0.136
VarietyTAM107      -1.0375     4.9791  -0.208
VarietyTAM200      -8.2000     4.9791  -1.647
VarietyVONA        -5.8375     4.9791  -1.172
> anova(rcb0.lmer)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
Variety 55 2387.5  43.409  0.8755
> fixef(rcb0.lmer)
(Intercept)      VarietyBRULE   VarietyBUCKSKIN    VarietyCENTURA
29.4375           -3.3625           -3.8750           -7.7875
VarietyCENTURK78   VarietyCHEYENNE       VarietyCODY       VarietyCOLT
0.8625           -1.3750           -8.2250           -2.4375
-4.9250           -1.8000           -5.3125           -0.8750
VarietyLANCOTA    VarietyNE83404    VarietyNE83406    VarietyNE83407
-2.8875           -2.0500           -5.1625           -6.7500
VarietyNE83432    VarietyNE83498    VarietyNE83T12    VarietyNE84557
-9.7125            0.6875           -7.8750           -8.9125
VarietyNE85556    VarietyNE85623    VarietyNE86482    VarietyNE86501
-3.0500           -7.7125           -5.1500            1.5000
VarietyNE86503    VarietyNE86507    VarietyNE86509    VarietyNE86527
3.2125           -5.6500           -2.5875           -7.4250
VarietyNE86582    VarietyNE86606    VarietyNE86607   VarietyNE86T666
-4.9000            0.3250           -0.1125           -7.9000
VarietyNE87403    VarietyNE87408    VarietyNE87409    VarietyNE87446
-4.3125           -3.1375           -8.0625           -1.7625
VarietyNE87451    VarietyNE87457    VarietyNE87463    VarietyNE87499
-4.8250           -5.5250           -3.5250           -9.0250
VarietyNE87512    VarietyNE87513    VarietyNE87522    VarietyNE87612
-6.1875           -2.6250           -4.4375           -7.6375
VarietyNE87613    VarietyNE87615    VarietyNE87619    VarietyNE87627
-0.0375           -3.7500            1.8250           -6.2125
VarietyNORKAN    VarietyREDLAND VarietyROUGHRIDER    VarietySCOUT66
-5.0250            1.0625           -8.2500           -1.9125
VarietySIOUXLAND     VarietyTAM107     VarietyTAM200       VarietyVONA
0.6750           -1.0375           -8.2000           -5.8375
> ranef(rcb0.lmer)
$Rep (Intercept) 1 1.8798700 2 2.8436747 3 -0.8713991 4 -3.8521455  nlme > rcb0.lme <- lme(yield~Variety, random=~1|Rep, data=na.omit(nin89)) > print(rcb0.lme, corr=FALSE) Linear mixed-effects model fit by REML Data: na.omit(nin89) Log-restricted-likelihood: -608.8508 Fixed: yield ~ Variety (Intercept) VarietyBRULE VarietyBUCKSKIN VarietyCENTURA 29.4375 -3.3625 -3.8750 -7.7875 VarietyCENTURK78 VarietyCHEYENNE VarietyCODY VarietyCOLT 0.8625 -1.3750 -8.2250 -2.4375 VarietyGAGE VarietyHOMESTEAD VarietyKS831374 VarietyLANCER -4.9250 -1.8000 -5.3125 -0.8750 VarietyLANCOTA VarietyNE83404 VarietyNE83406 VarietyNE83407 -2.8875 -2.0500 -5.1625 -6.7500 VarietyNE83432 VarietyNE83498 VarietyNE83T12 VarietyNE84557 -9.7125 0.6875 -7.8750 -8.9125 VarietyNE85556 VarietyNE85623 VarietyNE86482 VarietyNE86501 -3.0500 -7.7125 -5.1500 1.5000 VarietyNE86503 VarietyNE86507 VarietyNE86509 VarietyNE86527 3.2125 -5.6500 -2.5875 -7.4250 VarietyNE86582 VarietyNE86606 VarietyNE86607 VarietyNE86T666 -4.9000 0.3250 -0.1125 -7.9000 VarietyNE87403 VarietyNE87408 VarietyNE87409 VarietyNE87446 -4.3125 -3.1375 -8.0625 -1.7625 VarietyNE87451 VarietyNE87457 VarietyNE87463 VarietyNE87499 -4.8250 -5.5250 -3.5250 -9.0250 VarietyNE87512 VarietyNE87513 VarietyNE87522 VarietyNE87612 -6.1875 -2.6250 -4.4375 -7.6375 VarietyNE87613 VarietyNE87615 VarietyNE87619 VarietyNE87627 -0.0375 -3.7500 1.8250 -6.2125 VarietyNORKAN VarietyREDLAND VarietyROUGHRIDER VarietySCOUT66 -5.0250 1.0625 -8.2500 -1.9125 VarietySIOUXLAND VarietyTAM107 VarietyTAM200 VarietyVONA 0.6750 -1.0375 -8.2000 -5.8375 Random effects: Formula: ~1 | Rep (Intercept) Residual StdDev: 3.14371 7.041475 Number of Observations: 224 Number of Groups: 4 > anova(rcb0.lme) numDF denDF F-value p-value (Intercept) 1 165 242.05402 <.0001 Variety 55 165 0.87549 0.7119 > fixef(rcb0.lme) (Intercept) VarietyBRULE VarietyBUCKSKIN VarietyCENTURA 29.4375 -3.3625 -3.8750 -7.7875 VarietyCENTURK78 VarietyCHEYENNE VarietyCODY VarietyCOLT 0.8625 -1.3750 -8.2250 -2.4375 VarietyGAGE VarietyHOMESTEAD VarietyKS831374 VarietyLANCER -4.9250 -1.8000 -5.3125 -0.8750 VarietyLANCOTA VarietyNE83404 VarietyNE83406 VarietyNE83407 -2.8875 -2.0500 -5.1625 -6.7500 VarietyNE83432 VarietyNE83498 VarietyNE83T12 VarietyNE84557 -9.7125 0.6875 -7.8750 -8.9125 VarietyNE85556 VarietyNE85623 VarietyNE86482 VarietyNE86501 -3.0500 -7.7125 -5.1500 1.5000 VarietyNE86503 VarietyNE86507 VarietyNE86509 VarietyNE86527 3.2125 -5.6500 -2.5875 -7.4250 VarietyNE86582 VarietyNE86606 VarietyNE86607 VarietyNE86T666 -4.9000 0.3250 -0.1125 -7.9000 VarietyNE87403 VarietyNE87408 VarietyNE87409 VarietyNE87446 -4.3125 -3.1375 -8.0625 -1.7625 VarietyNE87451 VarietyNE87457 VarietyNE87463 VarietyNE87499 -4.8250 -5.5250 -3.5250 -9.0250 VarietyNE87512 VarietyNE87513 VarietyNE87522 VarietyNE87612 -6.1875 -2.6250 -4.4375 -7.6375 VarietyNE87613 VarietyNE87615 VarietyNE87619 VarietyNE87627 -0.0375 -3.7500 1.8250 -6.2125 VarietyNORKAN VarietyREDLAND VarietyROUGHRIDER VarietySCOUT66 -5.0250 1.0625 -8.2500 -1.9125 VarietySIOUXLAND VarietyTAM107 VarietyTAM200 VarietyVONA 0.6750 -1.0375 -8.2000 -5.8375 > ranef(rcb0.lme) (Intercept) 1 1.8795997 2 2.8432659 3 -0.8712739 4 -3.8515918  Model 1 ASReml-R > rcb.asr <- asreml(yield~Variety, random=~idv(Rep), rcov=~idv(units), data=nin89, na.method.X="include") > summary(rcb.asr)$call
asreml(fixed = yield ~ Variety, random = ~idv(Rep), rcov = ~idv(units),
data = nin89, na.method.X = "include")

$loglik [1] -454.4691$nedf
[1] 168

$sigma [1] 1$varcomp
gamma component std.error  z.ratio constraint
Rep!Rep.var  9.882911  9.882911  8.792823 1.123975   Positive
R!variance   1.000000  1.000000        NA       NA      Fixed
R!units.var 49.582368 49.582368  5.458839 9.082951   Positive

attr(,"class")
[1] "summary.asreml"
> summary(rcb0.asr)$varcomp gamma component std.error z.ratio constraint Rep!Rep.var 0.1993231 9.882911 8.792829 1.123974 Positive R!variance 1.0000000 49.582368 5.458839 9.082951 Positive > anova(rcb.asr) Wald tests for fixed effects Response: yield Terms added sequentially; adjusted for those above Df Sum of Sq Wald statistic Pr(Chisq) (Intercept) 1 242.054 242.054 <2e-16 *** Variety 55 48.152 48.152 0.7317 residual (MS) 1.000 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > coef(rcb.asr)$fixed
effect
Variety_ARAPAHOE    0.0000
Variety_BRULE      -3.3625
Variety_BUCKSKIN   -3.8750
Variety_CENTURA    -7.7875
Variety_CENTURK78   0.8625
Variety_CHEYENNE   -1.3750
Variety_CODY       -8.2250
Variety_COLT       -2.4375
Variety_GAGE       -4.9250
Variety_KS831374   -5.3125
Variety_LANCER     -0.8750
Variety_LANCOTA    -2.8875
Variety_NE83404    -2.0500
Variety_NE83406    -5.1625
Variety_NE83407    -6.7500
Variety_NE83432    -9.7125
Variety_NE83498     0.6875
Variety_NE83T12    -7.8750
Variety_NE84557    -8.9125
Variety_NE85556    -3.0500
Variety_NE85623    -7.7125
Variety_NE86482    -5.1500
Variety_NE86501     1.5000
Variety_NE86503     3.2125
Variety_NE86507    -5.6500
Variety_NE86509    -2.5875
Variety_NE86527    -7.4250
Variety_NE86582    -4.9000
Variety_NE86606     0.3250
Variety_NE86607    -0.1125
Variety_NE86T666   -7.9000
Variety_NE87403    -4.3125
Variety_NE87408    -3.1375
Variety_NE87409    -8.0625
Variety_NE87446    -1.7625
Variety_NE87451    -4.8250
Variety_NE87457    -5.5250
Variety_NE87463    -3.5250
Variety_NE87499    -9.0250
Variety_NE87512    -6.1875
Variety_NE87513    -2.6250
Variety_NE87522    -4.4375
Variety_NE87612    -7.6375
Variety_NE87613    -0.0375
Variety_NE87615    -3.7500
Variety_NE87619     1.8250
Variety_NE87627    -6.2125
Variety_NORKAN     -5.0250
Variety_REDLAND     1.0625
Variety_ROUGHRIDER -8.2500
Variety_SCOUT66    -1.9125
Variety_SIOUXLAND   0.6750
Variety_TAM107     -1.0375
Variety_TAM200     -8.2000
Variety_VONA       -5.8375
(Intercept)        29.4375
> coef(rcb.asr)$random effect Rep_1 1.8795997 Rep_2 2.8432658 Rep_3 -0.8712738 Rep_4 -3.8515916  nlme See the trick > nin89$Int <- 1
> rcb.lme <- lme(yield~Variety, random=list(Int=pdIdent(~Rep-1)), data=na.omit(nin89))
> print(rcb.lme, corr=FALSE)
Linear mixed-effects model fit by REML
Data: na.omit(nin89)
Log-restricted-likelihood: -608.8508
Fixed: yield ~ Variety
(Intercept)      VarietyBRULE   VarietyBUCKSKIN    VarietyCENTURA
29.4375           -3.3625           -3.8750           -7.7875
VarietyCENTURK78   VarietyCHEYENNE       VarietyCODY       VarietyCOLT
0.8625           -1.3750           -8.2250           -2.4375
-4.9250           -1.8000           -5.3125           -0.8750
VarietyLANCOTA    VarietyNE83404    VarietyNE83406    VarietyNE83407
-2.8875           -2.0500           -5.1625           -6.7500
VarietyNE83432    VarietyNE83498    VarietyNE83T12    VarietyNE84557
-9.7125            0.6875           -7.8750           -8.9125
VarietyNE85556    VarietyNE85623    VarietyNE86482    VarietyNE86501
-3.0500           -7.7125           -5.1500            1.5000
VarietyNE86503    VarietyNE86507    VarietyNE86509    VarietyNE86527
3.2125           -5.6500           -2.5875           -7.4250
VarietyNE86582    VarietyNE86606    VarietyNE86607   VarietyNE86T666
-4.9000            0.3250           -0.1125           -7.9000
VarietyNE87403    VarietyNE87408    VarietyNE87409    VarietyNE87446
-4.3125           -3.1375           -8.0625           -1.7625
VarietyNE87451    VarietyNE87457    VarietyNE87463    VarietyNE87499
-4.8250           -5.5250           -3.5250           -9.0250
VarietyNE87512    VarietyNE87513    VarietyNE87522    VarietyNE87612
-6.1875           -2.6250           -4.4375           -7.6375
VarietyNE87613    VarietyNE87615    VarietyNE87619    VarietyNE87627
-0.0375           -3.7500            1.8250           -6.2125
VarietyNORKAN    VarietyREDLAND VarietyROUGHRIDER    VarietySCOUT66
-5.0250            1.0625           -8.2500           -1.9125
VarietySIOUXLAND     VarietyTAM107     VarietyTAM200       VarietyVONA
0.6750           -1.0375           -8.2000           -5.8375

Random effects:
Formula: ~Rep - 1 | Int
Structure: Multiple of an Identity
Rep1    Rep2    Rep3    Rep4 Residual
StdDev: 3.14371 3.14371 3.14371 3.14371 7.041475

Number of Observations: 224
Number of Groups: 1
> anova(rcb.lme)
numDF denDF   F-value p-value
(Intercept)     1   168 242.05402  <.0001
Variety        55   168   0.87549  0.7121
> fixef(rcb.lme)
(Intercept)      VarietyBRULE   VarietyBUCKSKIN    VarietyCENTURA
29.4375           -3.3625           -3.8750           -7.7875
VarietyCENTURK78   VarietyCHEYENNE       VarietyCODY       VarietyCOLT
0.8625           -1.3750           -8.2250           -2.4375
-4.9250           -1.8000           -5.3125           -0.8750
VarietyLANCOTA    VarietyNE83404    VarietyNE83406    VarietyNE83407
-2.8875           -2.0500           -5.1625           -6.7500
VarietyNE83432    VarietyNE83498    VarietyNE83T12    VarietyNE84557
-9.7125            0.6875           -7.8750           -8.9125
VarietyNE85556    VarietyNE85623    VarietyNE86482    VarietyNE86501
-3.0500           -7.7125           -5.1500            1.5000
VarietyNE86503    VarietyNE86507    VarietyNE86509    VarietyNE86527
3.2125           -5.6500           -2.5875           -7.4250
VarietyNE86582    VarietyNE86606    VarietyNE86607   VarietyNE86T666
-4.9000            0.3250           -0.1125           -7.9000
VarietyNE87403    VarietyNE87408    VarietyNE87409    VarietyNE87446
-4.3125           -3.1375           -8.0625           -1.7625
VarietyNE87451    VarietyNE87457    VarietyNE87463    VarietyNE87499
-4.8250           -5.5250           -3.5250           -9.0250
VarietyNE87512    VarietyNE87513    VarietyNE87522    VarietyNE87612
-6.1875           -2.6250           -4.4375           -7.6375
VarietyNE87613    VarietyNE87615    VarietyNE87619    VarietyNE87627
-0.0375           -3.7500            1.8250           -6.2125
VarietyNORKAN    VarietyREDLAND VarietyROUGHRIDER    VarietySCOUT66
-5.0250            1.0625           -8.2500           -1.9125
VarietySIOUXLAND     VarietyTAM107     VarietyTAM200       VarietyVONA
0.6750           -1.0375           -8.2000           -5.8375
> ranef(rcb.lme)
Rep1     Rep2       Rep3      Rep4
1 1.8796 2.843266 -0.8712739 -3.851592


Model 2

ASReml-R

sp1.asr <- asreml(yield~Variety, rcov=~Column:ar1(Row), data=nin89, na.method.X="include")

> summary(sp1.asr)
$call asreml(fixed = yield ~ Variety, rcov = ~Column:ar1(Row), data = nin89, na.method.X = "include")$loglik
[1] -408.1412

$nedf [1] 168$sigma
[1] 7.975127

$varcomp gamma component std.error z.ratio constraint R!variance 1.0000000 63.6026561 11.3182328 5.619486 Positive R!Row.cor 0.7795799 0.7795799 0.0406026 19.200245 Unconstrained attr(,"class") [1] "summary.asreml" > summary(sp1.asr)$varcomp
gamma  component  std.error   z.ratio    constraint
R!variance 1.0000000 63.6026561 11.3182328  5.619486      Positive
R!Row.cor  0.7795799  0.7795799  0.0406026 19.200245 Unconstrained
> anova(sp1.asr)
Wald tests for fixed effects

Response: yield

Df Sum of Sq Wald statistic Pr(Chisq)
(Intercept)    1   24604.3         386.84 < 2.2e-16 ***
Variety       55    7974.4         125.38 2.048e-07 ***
residual (MS)         63.6
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> coef(sp1.asr)$fixed effect Variety_ARAPAHOE 0.0000000 Variety_BRULE -2.4048816 Variety_BUCKSKIN 7.8064972 Variety_CENTURA -1.6997427 Variety_CENTURK78 -1.3829446 Variety_CHEYENNE -1.1113084 Variety_CODY -6.7461911 Variety_COLT -1.7963394 Variety_GAGE -3.4539524 Variety_HOMESTEAD -5.5877510 Variety_KS831374 -0.8589476 Variety_LANCER -2.8418476 Variety_LANCOTA -5.9394801 Variety_NE83404 -3.4112613 Variety_NE83406 -1.9057358 Variety_NE83407 -3.2563922 Variety_NE83432 -5.4594311 Variety_NE83498 0.6446010 Variety_NE83T12 -4.0071361 Variety_NE84557 -4.2005181 Variety_NE85556 1.4836395 Variety_NE85623 -2.7617129 Variety_NE86482 -1.4309381 Variety_NE86501 -2.2287462 Variety_NE86503 -0.4557866 Variety_NE86507 -0.6983418 Variety_NE86509 -3.9215624 Variety_NE86527 0.5294386 Variety_NE86582 -5.4653632 Variety_NE86606 -0.7291575 Variety_NE86607 -0.1265536 Variety_NE86T666 -12.1437291 Variety_NE87403 -7.4623631 Variety_NE87408 -3.3586380 Variety_NE87409 -1.0360336 Variety_NE87446 -4.9030958 Variety_NE87451 -3.2836149 Variety_NE87457 -3.5244583 Variety_NE87463 -3.8427658 Variety_NE87499 -4.6298393 Variety_NE87512 -5.3760809 Variety_NE87513 -5.5656241 Variety_NE87522 -7.6500899 Variety_NE87612 -2.7225851 Variety_NE87613 -0.8793319 Variety_NE87615 -4.0089291 Variety_NE87619 0.7975626 Variety_NE87627 -10.1315147 Variety_NORKAN -7.1804945 Variety_REDLAND 0.6753066 Variety_ROUGHRIDER -0.9637487 Variety_SCOUT66 0.7088916 Variety_SIOUXLAND -1.1998807 Variety_TAM107 -3.7160351 Variety_TAM200 -9.0340942 Variety_VONA -2.7970689 (Intercept) 28.3487457  nlme Working on, yet not figured out. Could be something like this. Still could not figure out how to do rcov=~Column:ar1(Row) with nlme nin89$Int <- 1
sp1.lme <- lme(yield~Variety, random=~1|Int, data=na.omit(nin89))


Model 3

ASReml-R

sp2.asr <- asreml(yield~Variety, rcov=~ar1(Column):ar1(Row), data=nin89, na.method.X="include")

> summary(sp2.asr)
$call asreml(fixed = yield ~ Variety, rcov = ~ar1(Column):ar1(Row), data = nin89, na.method.X = "include")$loglik
[1] -399.3238

$nedf [1] 168$sigma
[1] 6.978728

$varcomp gamma component std.error z.ratio constraint R!variance 1.0000000 48.7026395 7.15527571 6.806536 Positive R!Column.cor 0.4375045 0.4375045 0.08060227 5.427943 Unconstrained R!Row.cor 0.6554798 0.6554798 0.05637709 11.626704 Unconstrained attr(,"class") [1] "summary.asreml" > summary(sp2.asr)$varcomp
gamma  component  std.error   z.ratio    constraint
R!variance   1.0000000 48.7026395 7.15527571  6.806536      Positive
R!Column.cor 0.4375045  0.4375045 0.08060227  5.427943 Unconstrained
R!Row.cor    0.6554798  0.6554798 0.05637709 11.626704 Unconstrained
> anova(sp2.asr)
Wald tests for fixed effects

Response: yield

> coef(sp2.asr)$fixed effect Variety_ARAPAHOE 0.00000000 Variety_BRULE 0.03029321 Variety_BUCKSKIN 8.89207227 Variety_CENTURA -0.68979639 Variety_CENTURK78 0.16461970 Variety_CHEYENNE 0.50267820 Variety_CODY -3.26960093 Variety_COLT -0.51826695 Variety_GAGE -0.95824999 Variety_HOMESTEAD -4.57873078 Variety_KS831374 0.27843476 Variety_LANCER -2.95379384 Variety_LANCOTA -4.67006598 Variety_NE83404 -1.32290865 Variety_NE83406 -1.66351994 Variety_NE83407 -2.64471830 Variety_NE83432 -4.42828427 Variety_NE83498 1.80418738 Variety_NE83T12 -2.11789109 Variety_NE84557 -2.34685080 Variety_NE85556 2.78001120 Variety_NE85623 -1.42164134 Variety_NE86482 -1.63334029 Variety_NE86501 -2.94339063 Variety_NE86503 -0.95747374 Variety_NE86507 0.46223383 Variety_NE86509 -3.27166458 Variety_NE86527 1.86588098 Variety_NE86582 -3.87940069 Variety_NE86606 0.22753741 Variety_NE86607 0.60702026 Variety_NE86T666 -10.27005825 Variety_NE87403 -7.43945904 Variety_NE87408 -3.10433009 Variety_NE87409 1.29746980 Variety_NE87446 -4.15943316 Variety_NE87451 -1.85324718 Variety_NE87457 -2.31156727 Variety_NE87463 -4.47086114 Variety_NE87499 -1.85909637 Variety_NE87512 -4.06473578 Variety_NE87513 -3.99604937 Variety_NE87522 -5.52109215 Variety_NE87612 -1.95543098 Variety_NE87613 -0.83160454 Variety_NE87615 -1.92104271 Variety_NE87619 2.98322047 Variety_NE87627 -7.33205188 Variety_NORKAN -5.78418023 Variety_REDLAND 1.75249392 Variety_ROUGHRIDER -0.97736288 Variety_SCOUT66 2.13126094 Variety_SIOUXLAND -2.54195346 Variety_TAM107 -1.59083563 Variety_TAM200 -6.54229161 Variety_VONA -1.52728371 (Intercept) 27.04285175  nlme Working on, yet not figured out. Could be something like this. Still could not figure out how to do rcov=~ar1(Column):ar1(Row) with nlme nin89$Int <- 1

I could not figured out how to fit Model 2 and 3 with nlme. I'm working on it and will update the answer when get it done. But I've included the output from ASReml-R for Model 2 and 3 for comparison purposes. Kevin has good experience of analyzing such models and Ben Bolker has wonderful authority on Mixed Models. I hope they can help us on Models 2 and 3.