[R-sig-ME] Question of ineraction in LMM model

[Paul Johnson]

On Mon, Jun 25, 2012 at 10:30 PM, chenlei <chenlei at ibcas.ac.cn> wrote:
> Dear lists,
>
> Greetings. I read message for a question in R list on May 29 from Dr. ir. Thierry Onkelinx. It is very helpful for me. It is as following:
>
> -------------------------
>
> Dear Luke,
>
> A workaround is to first create a new variable for the interaction.
>
> TimeTertile <- interaction(Time, Tertile, drop = TRUE)
>
> And then use that variable in the model
>
> m1 <- glmer(y ~ TimeTertile + (1|Individual) + (1|Zone) + Max + Min, family=binomial)
>
> Please not that I use glmer() instead of lmer() which give some more clear code. At the moment, lmer() will call glmer() for you when family is not Gaussian.
>
> Then you can run multicomp
>
> Timetest <- glht(m1, linfct = mcp(TimeTertile = "Tukey"))
>
>
>
> -------------------------
>
> But I encounter problem when I use your method to calculate my data. In my data, I tried to analyzed the effect of different successsional age of forest, environmental variable (represented by elevation) on species richness of forest. We sampled 12 forest plots (variable plot), and each successtional age has 3 forest plots (totally 4 successional age). We tried to analyze the effect of successional age and environment on species richness within each plot. We are interested in comparing species richness of forests at different age. So in my case, one of interactive variables is factor variables (successional age) and another is numeric variable (elevation). Following above-described method, I calculated as following:
>
> elev.age=interaction(elevation,age,drop=T);
>
>          model1=lmer(richness~elev.age+(1|plot));
>
> But the error message is:
>
>            Error in rep.int(c(1, numeric(n)), n - 1L) :
>
>            negative length vectors are not allowed.

When you see errors, we need information about possible causes.

Step one. Please report output from

> sessionInfo()

Step two. Put your data on a server somewhere so we can see it and run
your example code with it.  if the data is huge, make a small subset
we can test to reproduce your trouble.

Now, about your lmer usage:

Better to do it with a data frame, ie, put y, x1, x2 in dat, and run

m1 <- lmer(y ~ x1 + x2 + (1|x2), data = dat)

Generally, I get better results using a data frame.   Depending on the
routine to go find data in your workspace is risky, not just with
lmer, but with any regression function. Using a data frame gives you
security that it is finding what you want.

Although it may not cause an error now in R, it used to be that a
variable named "plot" could be confused with the commonly used
function "plot". When something goes wrong, I always worry about
things like that.  For the same reason, I never name variables "lm",
"rep", "seq", "plot" "glm" and so forth.  Those things used to cause R
to crash, and somewhere deep in my heart I'm still suspicious when
people do that.


>
> I also tried with
>
>         model2=lmer(richness~elev*age+(1|plot);
>
>  comp1<-glht(model2,linfct=(mcp(stage="Tukey",interaction_average=T,covariate_average=T)));
>
> But it gives error message:
>
>        Warning message:
>
> mean(<data.frame>) is deprecated.
>
>        Use colMeans() or sapply(*, mean) instead.
>
That means the author of glht needs to update code to not run
mean(data) any more. You are not causing that, unless your version of
the package is outdated.

You probably should show us the output of the summary of your lmer
model. Whether or not the hypothesis test is reasonable may depend on
the format of the output.  If you are confident the lmer output is
really what you want, then you need to send the lmer output to the
glht author and ask the usage question.  That's a big IF, though.

pj



> I checked it many times, but cannot find the right way. Could anyone please one give me any suggestions or hint? Your help will be highly appreciated!
>
> Best wishes,
>
> Xiangcheng Mi
>
>
>
> ---------------
>
> Xiangcheng Mi,
>
> Institute of Botany, Chinese Academy of Sciences,
>
> Beijing, 100093
>
> China
>
>
>
>
>
>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu