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

[Nicolay Cunha]

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)

#-----------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?

#############################################################################################
# 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.
#############################################################################################


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]]