R help - Dec 2006 - How to write a two-way interaction as a random effect in a lmer model?

Dear All,

 

I am working with linear mixed-effects models using the lme4 package in
R. I created a model with the lmer function including some main effects,
a two-way interaction and a random effect. Now I am searching how I
could incorporate an interaction between the random effect and one of
the fixed effects.

 

I tried to express the interaction in:

recap_random3<-lmer(breath~handling+stress+stress:handling+(1|rnr)+(0|rn
r:stress))

however R gives me the following error message

:

Error in eval(expr, envir, enclos) : fl[[2]] must be a factor of length
1744

In addition: Warning messages:

1: numerical expression has 1744 elements: only the first used in:
rnr:stress 

2: numerical expression has 1744 elements: only the first used in:
rnr:stress

 

If I fit this model in SPSS, this gives me output, but I don't know
whether I can trust that.

 

 

Therefore I went to look for an alternative:

After looking at the help function in R for lmer I deduced these models.


Examples from R

(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),
sleepstudy))
anova(fm1, fm2)

 

 

My models would then look like this:

 

recap_random0<-lmer(breath~handling+stress (1|rnr))

recap_random1<-lmer(breath~handling+stress (1|rnr)+(0+stress|rnr))

 

However, I do not know how to interpret the results. Does the model see
stress|rnr as an interaction between stress and rnr, or did it take
stress nested in rnr.

 
> summary(recap_random1)
Linear mixed-effects model fit by REML 

Formula: breath ~ handling + stress + stress:handling + (1 | rnr) + (0 +
stress | rnr) 

  AIC  BIC logLik MLdeviance REMLdeviance

 9748 9781  -4868       9719         9736

Random effects:

 Groups   Name        Variance Std.Dev.

 rnr      (Intercept) 55.6711  7.4613                       # BUT

 rnr      stress      15.2805  3.9090                       # Does this
output line express the interaction??? 

 Residual              7.4567  2.7307                    # If not, how
should the interaction output look like?

number of obs: 1744, groups: rnr, 217; rnr, 217

 

Fixed effects:

                 Estimate Std. Error t value

(Intercept)     28.164471   0.716359   39.32

handling         0.012477   0.012330    1.01

stress           0.636979   0.416011    1.53

 

Correlation of Fixed Effects:

            (Intr) hndlng stress

handling    -0.645              

stress      -0.517  0.667       

 

 

Could anybody please give me an advice how to solve this problem? Which
way is correct to express interaction in the random factors? How should
the output look like?

 

Thank you in advance,

Eva Fucikova

 

********************************************

Msc. Eva Fucikova 
Netherlands Institute of Ecology (NIOO-KNAW) 
PO Box 40 
6666 ZG Heteren 
The Netherlands 
tel.: +31 (0)26 4791248 
Email:E.Fucikova@nioo.knaw.nl 

*******************************************

 


	[[alternative HTML version deleted]]

Hi Eva,
A couple questions:
Are repeated measurements taken on rnr? Is rnr "subject"?
Is stress a continuous variable?

See below.
On Dec 11, 2006, at 11:04 AM, Fucikova, Eva wrote:
> Dear All,
>
>
>
> I am working with linear mixed-effects models using the lme4  
> package in
> R. I created a model with the lmer function including some main  
> effects,
> a two-way interaction and a random effect. Now I am searching how I
> could incorporate an interaction between the random effect and one of
> the fixed effects.
>
>
>
> I tried to express the interaction in:
>
> recap_random3<-lmer(breath~handling+stress+stress:handling+(1|rnr) 
> +(0|rn
> r:stress))
>
lmer will want rnr:stress to be a factor, but (0|factor) doesn't make  
sense. To test for random variation in slopes, you want (stress |  
rnr), (or (1|rnr) + (0+stress | rnr) for uncorrelated slopes).
> however R gives me the following error message
>
> :
>
> Error in eval(expr, envir, enclos) : fl[[2]] must be a factor of  
> length
> 1744
>
> In addition: Warning messages:
>
> 1: numerical expression has 1744 elements: only the first used in:
> rnr:stress
>
> 2: numerical expression has 1744 elements: only the first used in:
> rnr:stress
>
>
>
> If I fit this model in SPSS, this gives me output, but I don't know
> whether I can trust that.
>
>
>
>
>
> Therefore I went to look for an alternative:
>
> After looking at the help function in R for lmer I deduced these  
> models.
>
>
> Examples from R
>
> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
> (fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject),
> sleepstudy))
> anova(fm1, fm2)
>
>
>
>
>
> My models would then look like this:
>
>
>
> recap_random0<-lmer(breath~handling+stress (1|rnr))
>
> recap_random1<-lmer(breath~handling+stress (1|rnr)+(0+stress|rnr))
>
Or rather with a plus sign
>
> recap_random0<-lmer(breath~handling+stress + (1|rnr))
>
> recap_random1<-lmer(breath~handling+stress + (1|rnr)+(0+stress|rnr))
or
>
> recap_random0<-lmer(breath~handling+stress (stress | rnr))
>
>
> However, I do not know how to interpret the results. Does the model  
> see
> stress|rnr as an interaction between stress and rnr, or did it take
> stress nested in rnr.
>
>
>
>> summary(recap_random1)
>
> Linear mixed-effects model fit by REML
>
> Formula: breath ~ handling + stress + stress:handling + (1 | rnr) +  
> (0 +
> stress | rnr)
>
>   AIC  BIC logLik MLdeviance REMLdeviance
>
>  9748 9781  -4868       9719         9736
>
> Random effects:
>
>  Groups   Name        Variance Std.Dev.
>
>  rnr      (Intercept) 55.6711  7.4613                       # BUT
>
>  rnr      stress      15.2805  3.9090                       # Does  
> this
> output line express the interaction???
>
If the model statement was correct, these would be variances  
associated with the intercepts and slopes of different levels of rnr.  
Thus, I would call that rnr by stress interaction.

Cheers,
Hank
>  Residual              7.4567  2.7307                    # If not, how
> should the interaction output look like?
>
> number of obs: 1744, groups: rnr, 217; rnr, 217
>
>
>
> Fixed effects:
>
>                  Estimate Std. Error t value
>
> (Intercept)     28.164471   0.716359   39.32
>
> handling         0.012477   0.012330    1.01
>
> stress           0.636979   0.416011    1.53
>
>
>
> Correlation of Fixed Effects:
>
>             (Intr) hndlng stress
>
> handling    -0.645
>
> stress      -0.517  0.667
>
>
>
>
>
> Could anybody please give me an advice how to solve this problem?  
> Which
> way is correct to express interaction in the random factors? How  
> should
> the output look like?
>
>
>
> Thank you in advance,
>
> Eva Fucikova
>
>
>
> ********************************************
>
> Msc. Eva Fucikova
> Netherlands Institute of Ecology (NIOO-KNAW)
> PO Box 40
> 6666 ZG Heteren
> The Netherlands
> tel.: +31 (0)26 4791248
> Email:E.Fucikova at nioo.knaw.nl
>
> *******************************************
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting- 
> guide.html
> and provide commented, minimal, self-contained, reproducible code.


Hank Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/
"E Pluribus Unum"

If you send an attachment, please try to send it in a format anyone  
can read, such as PDF, text, Open Document Format, HTML, or RTF.  
Please try not to send me MS Word or PowerPoint attachments-
Why? See:  http://www.gnu.org/philosophy/no-word-attachments.html

Hi folks,
I am not providing a small replicable example, because I assume the  
problem is related to my quirky data (~650 obs.).
I am using the latest lme4, matrix and coda, and R 2.4.0.

I frequently get the following error message for this particular lmer  
model. I do not get this message for the vast majority of my uses of  
mcmcsamp.

I was wondering if it might be characteristic of particular kinds of  
problems.

 > modb <- lmer(log(basal+1) ~ nutrient*amd +  (1|rack) +  (1| 
gen:amd), data=dat.b2)

 > mod.mc <- mcmcsamp(modb, n=10^4)

Error: Matrix is not pd after safe_pd_matrix!
Error in t(.Call(mer_MCMCsamp, object, saveb, n, trans, verbose)) :
	error in evaluating the argument 'x' in selecting a method for  
function 't'

Cheers,
Hank

Dr. Hank Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/~stevenmh/
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/

"E Pluribus Unum"