# How to find the best fixed effect in a mixed model?

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.

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.

• Why do you have treatment as a random effect? I would think you should only have Block as fixed effect. – svendvn Apr 17 at 21:46
• 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? – andemexoax Apr 17 at 22:16
• 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. – svendvn Apr 18 at 0:06
• Oops, I meant that I thought that the treatment should only be fixed! Apoligies – svendvn Apr 18 at 0:08

• 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? – andemexoax Apr 17 at 22:31