Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Converting mixed model formula from SAS to R

Tags:

r

sas

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
like image 716
jon Avatar asked Oct 10 '22 13:10

jon


1 Answers

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

like image 86
Aaron left Stack Overflow Avatar answered Oct 13 '22 09:10

Aaron left Stack Overflow