I have a dataset with growth rate as a response variable (
resp in the example) and temperature, food availability, and salinity as predictor variables (
pred3 in the example). The predictor variables are "continuous" with weekly intervals and have different units. Measurements span weekly (with missing values for some samples) throughout a year (
week in the example; defined from the beginning of the experiment). I have several samples and I want to quantify (over all samples):
- How much each predictor variable explains the variation in growth rate
- The relative effect of each predictor variable on growth rate
I understand that linear mixed models could be a solution for this problem as I have several samples and dependent measurements over time. My question is: What would be the optimal model formulations using
lme4 package for R?
Example data is available here. And here is an overview of it:
library(ggplot2) tmp <- melt(X, id = c("Sample", "weeks")) ggplot(tmp, aes(x = weeks, y = value)) + geom_line() + facet_wrap(Sample ~ variable, scales = "free_y")
I have tried following:
As a solution for point 1:
library("lme4") library("MuMIn") p1 <- lmer(resp ~ pred1 + (1|Sample) + (1|weeks), data = X) p2 <- lmer(resp ~ pred2 + (1|Sample) + (1|weeks), data = X) p3 <- lmer(resp ~ pred3 + (1|Sample) + (1|weeks), data = X) margr2 <- data.frame(Pred = c("pred1", "pred2", "pred3"), marginal.R2 = c(r.squaredGLMM(p1)[], r.squaredGLMM(p2)[], r.squaredGLMM(p3)[])) ggplot(margr2, aes(x = Pred, y = marginal.R2)) + geom_bar(stat = "identity")
Marginal $R^2$ calculated by the method published here should indicate the overall variance explained by each predictor variable as far as I have understood and assuming that my model formulations are correct.
For the relative effect (point 2), I think that I first have to have the predictor variables on a same scale. Only then can I compare them by having them all in the model and removing the intercepts:
Xs <- X Xs[4:6] <- scale(Xs[4:6]) mod <- lmer(resp ~ pred1 + pred2 + pred3 - 1 + (1|weeks) + (1|Sample), data = Xs) cis <- confint(mod)[4:6,] releff <- data.frame(par = rownames(cis), lower = cis[,1], est = fixef(mod), upper = cis[,2])
In order to make the interpretation more intuitive, I scale the effects to maximum absolute value of across confidence intervals (I am only interested in relative effect):
tmp <- c(releff$lower,releff$upper) add <- 100*releff[c("lower", "est", "upper")]/max(abs(tmp)) colnames(add) <- paste0("rel.", colnames(add)) releff <- cbind(releff, add) ggplot(releff, aes(x = par, y = rel.est, ymin = rel.lower, ymax = rel.upper)) + geom_pointrange() + geom_hline(yintercept = 0)
Predictor variables are "significant", where the CIs do not cross the horizontal line (to my understanding). I am not sure whether these approaches make much sense and that's why I am asking for help.