I want to fit a mixed model using nlme package in R which is equivalent to following SAS codes:
proc mixed data = one;
class var1 var2 year loc rep;
model yld = var1 * var2;
random loc year(loc) rep*year(loc);
EDITS: Explanation of what is experiment about
the same combination of var1 and var2 were tested in replicates (rep- replicates are numbered 1:3). The replicates (rep) is considered random. This set of experiment is repeated over locations (loc) and years (year). Although replicates are numbered 1:3 within each location and year for covinience because they do not have any name, replication 1 within a location and a year doesnot have correlation replication 1 within other location and other year
I tried the following codes:
require(nlme)
fm1 <- lme(yld ~ var1*var2, data = one, random = loc + year / loc + rep * year / loc)
Is my codes correct?
EDITS: data and model based on suggestions you can download the example data file from the following link: https://sites.google.com/site/johndatastuff/mydata1.csv
data$var1 <- as.factor(data$var1)
data$var2 <- as.factor(data$var2)
data$year <- as.factor(data$year)
data$loc <- as.factor(data$loc)
data$rep <- as.factor(data$rep)
following suggestions from the comments below:
fm1 <- lme(yld ~ var1*var2, data = data, random = ~ loc + year / loc + rep * year / loc)
Error in getGroups.data.frame(dataMix, groups) :
Invalid formula for groups
EXPECTED BASED ON SAS OUTPUT
Type 3 tests of fixed effects
var1*var2 14 238 F value 16.12 Pr >F = < 0.0001
Covariance parameters:
loc = 0, year(loc) = 922161, year*rep(loc) = 2077492, residual = 1109238
I tried the following model, I still getting some errors:
Edits: Just for information I tried the following model
require(lme4)
fm1 <- lmer(yld ~ var1*var2 + (1|loc) + (1|year / loc) + (1|rep : (year / loc)),
data = data)
Error in rep:`:` : NA/NaN argument
In addition: Warning message:
In rep:`:` : numerical expression has 270 elements: only the first used
Thanks for the more detailed information. I stored the data in d
to avoid confusion with the data
function and parameter; the commands works either way but this avoiding data
is generally considered good practice.
Note that the interaction is hard to fit because of the lack of balance between var
and var2
; for reference here's the crosstabs:
> xtabs(~var1 + var2, data=d)
var2
var1 1 2 3 4 5
1 18 18 18 18 18
2 0 18 18 18 18
3 0 0 18 18 18
4 0 0 0 18 18
5 0 0 0 0 18
Normally to just fit the interaction (and no main effects) you'd use :
instead of *
, but here it works best to make a single factor, like this:
d$var12 <- factor(paste(d$var1, d$var2, sep=""))
Then with nlme
, try
fm1 <- lme(yld ~ var12, random = ~ 1 | loc/year/rep, data = d)
anova(fm1)
and with lme4
, try
fm1 <- lmer(yld ~ var12 + (1 | loc/year/rep), data=d)
anova(fm1)
Also note that because nlme
and lme4
have overlap in their function names you need to only load one at time into your R session; to switch you need to close R and restart. (Other ways exist but that's the simplest to explain.)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With