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 fixed effects, X is an n×p design matrix of full column rank which associates observations with the appropriate combination of fixed effects, u is the q × 1 vector of random effects, Z is the n × q design matrix which associates observations with the appropriate combination of random effects, and e is the n × 1 vector of residual errors.The model (1) is called a linear mixed model or linear mixed effects 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 effects 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 fixed to one. In mixed effects 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 effects: 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."
lme4. Can you (a) tell us why you need to do this inlme4rather thanasreml-R(b) consider posting onr-sig-mixed-modelswhere there is more relevant expertise? $\endgroup$corStructinnlme(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 ... $\endgroup$MCMCglmm, and I'm pretty sure that (other than thespatialCovariancementioned, which I'm unfamiliar with) the only way to get it done in R is by defining newcorStructs -- which is possible, but not trivial. $\endgroup$