Converting mixed model formula from SAS to R

Question

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

What happened when you tried it? Can you make a small reproducible example? — Aaron left Stack Overflow, Nov 07 '11 at 02:55
Did you try adding a ~ in the formula? fm1 <- lme(yld ~ var1*var2, data = one, random = ~ loc + year / loc + rep * year / loc) — Max Gordon, Nov 07 '11 at 09:03
Hmm... I haven't used lme/lmer but anything more advanced than lmer(yld ~ var1*var2 +(1|year), data=data) fails. — Max Gordon, Nov 07 '11 at 14:15
The bottom failure is because of "rep : ( year" the colon is placed up against a set of brackets. This is not what R expects. I don't know about SAS but it seems you might want something like: lmer(yld ~ var1*var2 + (1|loc) + (1|loc/year) + (1|loc/year:rep)), data = data)?? — , Nov 07 '11 at 15:23
I couldn't get the above to work when I tried it. You should probably label your replicates explicitly, perhaps by "data$rly <- with(data, interaction(rep, loc, year))" and then try "lmer1<-lmer(yld ~ var1*var2 + (1|loc) + (1|year)+(1|rlc), data = data)" — , Nov 07 '11 at 15:33
@John Don't cross-post over SE sites, please. As a rule of thumb, pure R questions should go on SO, stats questions belong here; otherwise, select one site and flag your question for mod attention. I'm moving this on SO, it will be merged with your other question to keep the history of comments. (Thanks Aaron for noticing that!) — chl, Nov 07 '11 at 16:24
(For future reference, @chl's "here" refers to stats.stackexchange.com, where this question originated, not stackoverflow, where it ended up.) — Aaron left Stack Overflow, Nov 07 '11 at 16:45
my apologies, I tried to move, the question as i was not able to get solutions....but did not know...then just leveled as moved ...thank you for handling the situation — jon, Nov 07 '11 at 17:12

Aaron left Stack Overflow · Accepted Answer · 2011-11-07T16:43:26.967

4

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.)

edited Nov 07 '11 at 16:43

answered Nov 07 '11 at 15:25

Aaron left Stack Overflow

36,704
7
77
142

thanks; if you look at the SAS model, my interest is mainly on the interaction term not on main effect, so even I do not worry fitting main effect model for var1 and var2 – jon Nov 07 '11 at 15:35
unbalance is purposeful, the expected product combination between 1-2 is same as 2-1, similarly 1-3 is same as 3-1. That is why you find 0 in lower triangle of the you Xtab. – jon Nov 07 '11 at 16:15
1

I believe that `detach(package:lme4)` would be enough (if, e.g. `lme4` has been loaded after `nlme` and we now want to use `nlme::lmList`), but you already know this, I guess. – chl Nov 07 '11 at 17:46
@chl: Actually, I'd forgotten how to do it another way (though knew it was possible), so thanks! – Aaron left Stack Overflow Nov 09 '11 at 18:18

Converting mixed model formula from SAS to R

1 Answers1