1
$\begingroup$

I am trying to see the difference between two variables in mixed models.

Background: The models are for an experiment with a randomized complete block design, but there is a bit of missing data. Yield is being estimated by the fixed effect of 31 different treatments with 1 to 3 replicates in 4 blocks. The additional fixed effects that I am trying to sort through measure several parameters of each treatment (height, diameter, width, height_2--a different method of measuring height--, etc). My basic model starts with the following formula yield ~ Treatment + (Treatment:Block) using the lme4 function lmer in R. I removed the block effect and am only including the interaction because when the models got more complicated I was warned of singular fits, and one of the suggestions I read about on how to deal with that was to simplify the model. I should probably have deleted the interaction and kept only the block effect, but this is beside the point for my question here.

Problem

My problem comes from when I am comparing two models like the following,

  1. a <- lmer(yield ~ Treatment + (1|Treatment:Block) + height, data=df)
  2. b <- lmer(yield ~ Treatment + (1|Treatment:Block) + height_2, data=df)

I want to know if one effect (height or height_2) explain yield better when controlling for treatment and block effects. When I run anova(a,b) the models are refit to ML instead of REML (which is typical), and the table shows Chi Df of 0 and a Pr(>Chisq) of 1. I don't think I am doing this appropriately.
enter image description here

I read somewhere that to see if fixed effects are different make a model including both and another similar model but include an interaction of the two effects like below. From the anova(c,d) I should know whether the variables (height and height_2) have different effects on yield...(I think)

  • . c <- lmer(yield ~ Treatment + (1|Treatment:Block) + height + height_2, data=df)
  • . d <- lmer(yield ~ Treatment + (1|Treatment:Block) + height + height_2 + (height:height_2), data=df)

Is the purpose of the previous comparison correct as I wrote it?

I still need to know: How to test whether one variable is a better predictor compared to the other if my first method is incorrect?

Bonus question: Are my variables (all or some) crossed or nested? Am I coding them correctly in lmer()... I read the answer to this question and I am still a bit fuzzy on how it applies to my situation.

| cite | improve this question | |
$\endgroup$
  • $\begingroup$ Why do you have treatment as a random effect? I would think you should only have Block as fixed effect. $\endgroup$ – svendvn Apr 17 '19 at 21:46
  • $\begingroup$ I have not intended to have treatment as a random effect. Does it become random as a result of the interaction term? The situation is that the treatments acted differently in different Blocks and I need to control for that treatment effect before I can look at the other parameters. When more than 2 of the fixed parameters were added I ran into the singular fit issues and block by itself didn't improve the model, but the interaction did. Why do you say to have only block as a fixed effect? $\endgroup$ – andemexoax Apr 17 '19 at 22:16
  • $\begingroup$ Oh, if the treatment effects differ between the blocks, that is not so good. In that case I don't know how to analyze the data best. $\endgroup$ – svendvn Apr 18 '19 at 0:06
  • $\begingroup$ Oops, I meant that I thought that the treatment should only be fixed! Apoligies $\endgroup$ – svendvn Apr 18 '19 at 0:08
0
$\begingroup$

The first method of comparing the two heights is wrong because the two models are not nested as anova requires. Instead you can just compare the likelihood values of the two models and choose the height variable that produces the highest likelihood (which is equivalent to AIC here because the two models have the same number of free parameters). In your dataset, model A has a higher likelihood (-727.1 compared to -752.7), meaning that height predicts the yield better than heigjt_2.

The second method you suggest compares the height variables by defining a model with both of them and testing if either height variables can be removed from the model (hope I understood that correctly). That is okay, but it can be inconclusive if both of them can be removed from model.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Are you suggesting I do ab_compare<-anova(a,b) then only look at ab_compare$logLik for comparing? If not then how do I get the likelihood value? If yes then how come the likelihood of ab_compare is correct without having the two models nested correctly? $\endgroup$ – andemexoax Apr 17 '19 at 22:31
  • $\begingroup$ Yes that you can do. Anova tries to calculate the chisquare likelihood ratio tests between the two models, and that it can't do because the models are not nested. Anova also reports the likelihood values and yes, those you can use to compare the models even though the models are not nested. (I wrote that the anova was wrong because the chisquare test was impossible) $\endgroup$ – svendvn Apr 18 '19 at 0:02
  • $\begingroup$ How could I compare whether the differences in log likelihood are "significant" through some sort of p<0.05 approach? $\endgroup$ – andemexoax Apr 18 '19 at 1:51
  • $\begingroup$ It is not really possible to get a P value for the comparison here $\endgroup$ – svendvn Apr 18 '19 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.