[R-sig-ME] Doubt about model structure and extraction of slopes and intercepts

[Phillip Alday]

On 12/3/20 8:48 pm, Nicolay Cunha wrote:
> Hi everyone,
>
> I am trying to fit a mixed model with a data set that was extracted from a
> grid. We grid a certain area and counted the frequency of occurrence of three
> levels of a factor variable in each grid (fac = fac1, fac2 and fac3), I am
> basically interested in testing the relationship of the frequency of
> occurrence of each level of “fac” with external variables (X1 and X2, for
> instance). The issue is, I am not interested in knowing if the frequency of
> occurrence of each level of “fac” differs between each other (fac1 ¹ fac2 ¹
> fac3) because it may be influenced by sampling bias, but I am interested if
> X1 and X2 explain the frequency of occurrence of fac1, fac2 and fac3. In
> this case, my sample unit is each level of fac per cell of the grid, and
> grid is my random factor (in my original data, there is one additional
> complication. In my case, “fac”, is a trait and each sample is a species
> individual occurring in the grid, so I will also have to correct for this,
> but I think that it can be put aside at the moment).
>
>
> Below is a reproducible example:
>
> ############################################################
>
> set.seed(457)
> # generating the data
> dat <- data.frame(X1 = runif(150,-2,2), X2 = runif(150,-2,2), fac = gl(n =
> 3,k = 50))
> modmat <- model.matrix(~X1*X2*fac,dat)
> betas <- runif(12,-2,2)
> dat$y <- rnorm(150, modmat%*%betas, 1 ) # frequency of occurrence
> dat$rand <- rep(c(1:5), each = 30) # random effect
>
> library(lme4)
> m1 <- lmer(y ~ fac/X1 + fac/X2 + (1|rand), data = dat, REML = FALSE)
> library(car)
> Anova(m1)
> summary(m1)



Thanks for having a MWE!

> #-----------First question--------------#
> # Does it make sense to write the model the way it is – by nesting fac with
> the predictor variables, to assess the relationship I am interested in?
> # Is there a better(or correct) way for doing this?

If I undertstand correctly, you're actually dealing with counts, right?
Wouldn't a Poisson model make more sense?

The nesting effectively works out to a particular contrast coding.
Contrast coding in mixed models works the same way as in classical
linear regression, so any resources you can find for that will also help
you here.

The problem I see here is that you need to think about whether you care
about absolute or relative differences. Right now, your interactions are
for absolute differences, e.g. you're correcting each fac for X1, but if
the baseline counts for each fac are drastically different, this may not
be the comparison you want.

> #############################################################################################
> # plotting the results
>
> with(dat, plot(X1, predict(m1), pch=21, col = "black", bg = c("green",
> "purple", "salmon")[dat$fac]))
> X1_seq = seq(min(dat$X1), max(dat$X1), by=0.001)
> ypred_X1_fac1 = (fixef(m1)[1])  + (fixef(phylo_lmm_LF4)[4])*(X1_seq)
> lines(X1_seq, ypred_X1_fac1, col = "green", lty = "solid", lwd = 1.5)
> ypred_X1_fac2 = (fixef(m1)[2])  + (fixef(phylo_lmm_LF4)[5])*(X1_seq)
> lines(X1_seq, ypred_X1_fac2, col = "purple", lty = "solid", lwd = 1.5)
> ypred_X1_fac3 = (fixef(m1)[3])  + (fixef(phylo_lmm_LF4)[6])*(X1_seq)
> lines(X1_seq, ypred_X1_fac3, col = "salmon", lty = "solid", lwd = 1.5)
>
> #-----------Second question--------------#
> # How to extract the slopes and intercepts for all combinations between X1,
> X2 and fac in model this way?
> # I tried plotting this but failed miserably and I am sure that I am
> committing some terrible mistake.
> #############################################################################################


Check out the effects and emmeans packages.

Phillip


>
> I do thank any help, advise and orientation with the model building and
> plotting the predictions.
>
> Cheers,
> Nicolay
>
>
> ------------------------------------------------------------------------
> Grupo de Ecología de la Polinización
> (https://sites.google.com/view/ecopol/home)
> INIBIOMA, CONICET-Universidad Nacional del Comahue
> Quintral 1250
> 8400 San Carlos de Bariloche
> Rio Negro, Argentina
> ------------------------------------------------------------------------
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models