[R-sig-ME] Hierarchical modelling using lmer v lme and including weights

[ONKELINX, Thierry]

Dear Tom,

> library(psych);  library(lme4);  library(nlme)

Loading lme4 and nlme together is not a good idea. It will lead to conflicts due to function with the same name.

> #Q1 - CAN I GET P VALUES? (Treat0.1 - looks 'significant' (P<0.05) but only just

FAQ 7.35: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-are-p_002dvalues-not-displayed-when-using-lmer_0028_0029_003f

> #model_1 shows difference variances for different levels of Treat, try varIdent
> vf1=varIdent(form=~1|Treat)
> model_2=lmer(LogG~Treat+(1|Tank)+(1|Pot), weights=vf1,
> data=data4);summary(model_2)#OR model_2=lmer(LogG~Treat+(1|Tank/Pot),
> weights=vf1, data=data4);summary(model_2) # Q2 - model_2 fails, why the
> error message - '...variable lengths differ'?

varIdent() is from the nlme package. It does not work with lme4

> #Different approach
> #USE lme
> model_3=lme(LogG~Treat,random=~1|Tank+1|Pot,weights=vf1,data=data4);
> summary(model_3)
> #Q3 - Model_3 'runs' but get error message '....+ not meaningful for factors'.
> Why?

Crossed random effects are rather hard to code in lme. You'll need one of the pdClasses. pdBlocked() I think. I haven't done that myself in lme() so I can't help you with that.
Futhermore unless you moved the pots between tanks (thus pot A from tank 1 is the same as pot A from tank 2) you should use nested random effects instead of crossed random effects.

> #respecify random component - looks more encouraging:
> model_3=lme(LogG~Treat,random=~1|Tank/Pot,weights=vf1,data=data4);
> summary(model_3)
> #Q3.1 - Should I assume from these results that the variances components
> indicate that most of the variance is occurring at the residual level # i.e. the level
> of individual measurements within each pot and that, therefore, I could usefully
> additionally replicate these measurements?

IMHO that is correct

> 
> 
> #Q3.2 can I get appropriate degrees of freedom with the parameter estimates
> (for fixed effects)?  This would reassure that correct analysis is being conducted?

summary(model_3)$tTable

No. It up to you to decide if ou have an acceptable model or not.

Best regards,

Thierry