[R-sig-ME] Incorrect output when applying lmer within a loop.

[Ben Bolker]

Anton Wasson <anton.wasson at ...> writes:

>  Hello All, I am a novice user of R and LME4 - I hope I am sending
> this query to an appropriate forum. I am generating erroneous data,
> but I don't know if the reason is misuse of LME4 or the loop
> function.

As it turns out, it's more of an R issue than an lme4 issue.

[snip]
 
> I wrote a loop to run a model (VAR ~ (1|GEN) + (1|OPR) + (1|ROW) +
> (1|RAN)) on each response variable, to extract the variance for each
> factor (and the residual) and to build a dataframe of the variances
> across all the response variables. The names of the variables were
> in the vector "nam".
 
> > nam [1] "X20" "X30" "X40" "X50" "X60" "X70" "X80" "X90" "X100"
>  "X110" "X120" "X130" "X140" "X150" "X160" "X170" 
 "X180" "X190"
>  "X200" "md" "tc" "md90"

Or 

nam <- c(paste0(X,(2:20)*10),c("md","tc","md90"))

[snip]

> x<-1
> for(i in nam) {
>   d<-r[,c("REP","ROW","RAN","GEN","OPR",i)]
>   colnames(d)<-c("REP","ROW","RAN","GEN","OPR","VAR")
>   r.mer<-lmer(VAR ~ (1|GEN) + (1|OPR) + (1|ROW) + (1|RAN), data=d)
>   r.var <- data.frame(summary(r.mer)@REmat)
>   V1[x]<-as.numeric(as.character(r.var[r.var$Groups=="GEN",'Variance']))
>   V2[x]<-as.numeric(as.character(r.var[r.var$Groups=="OPR",'Variance']))
>   V3[x]<-as.numeric(as.character(r.var[r.var$Groups=="RAN",'Variance']))
>   V4[x]<-as.numeric(as.character(r.var[r.var$Groups=="ROW",'Variance']))
>   V5[x]<-as.numeric(as.character(r.var[r.var$Groups=="Residual",'Variance']))
>   x<-x+1
>   }

 It would be cleaner and safer to use unlist(VarCorr(r.mer)) to access the RE
variances,
and sigma(fm1)^2 to get the residual variance.

  You can use 

    reformulate(paste0("(1|",c("GEN","OPR","ROW","RAN"),")"),response="V1")

## V1 ~ (1 | GEN) + (1 | OPR) + (1 | ROW) + (1 | RAN)

to generate formulae with different reponse variables, and use the same
input data set every time.

> df<-data.frame(nam,V1,V2,V3,V4,V5)
> colnames(df)<-c("DEP", "GEN", "OPR", "RAN", "ROW","Residual")
> 
> That script gave the following output for the variances for the factors
> 

[snip]

Your problem is that R stops running the code when you try to run the
example with the error, leaving the last four rows equal to zero as
they were initialized.  It is generally safer to fill in NA in the
data frame that will hold the results, so that you can see what hasn't
been computed.

d <- as.data.frame(matrix(NA,nrow=length(nam),ncol=6),
                   dimnames=list(nam,c("REP","ROW","RAN","GEN","OPR","VAR")))

You could use try() to intercept errors, but you would still have to
check the results (if inherits(model_results,"try-error")) to skip
the row appropriately