Daji
Reputation: 3
Predicting with lmer: errors and missing p-values
In an experiment, female fish were exposed to two levels of photoperiod (Ambient & Compressed), two levels of temperature (4 & 7). They were in four tanks (two tanks for each photoperiod, one tank for each temperature within photoperiod). There were nine samplings denoted by time_date in the data. Among other responses is "k". My interest is on the effects of photoperiod, temperature and time_date on "k". Challenges faced: Unbalanced design (one photoperiod or temperature level not sampled during a sampling), pseudo-replication (each tank is a treatment (temperature masked within photoperiod)). with some reading, I came across the mixed models. I have tried with lmer (more importantly: I am not sure if am right) and fell into warnings and outputs with no p-values. I appreciate your help. Thank you in advance.
Here is the sample data
fem.fish <- structure(list(time_date = structure(c(8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L,
14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L), .Label = c("30-Jan-18",
"11-Apr-18", "13-Jun-18", "07-Aug-18", "19-Sep-18", "30-Oct-18",
"28-Nov-18", "03-Jan-19", "17-Jan-19", "31-Jan-19", "14-Feb-19",
"28-Feb-19", "14-Mar-19", "27-Mar-19", "10-Apr-19", "24-Apr-19"
), class = "factor"), photo = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Ambient",
"Compress"), class = "factor"), temp = structure(c(2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("4",
"7"), class = "factor"), tank = structure(c(2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("T1",
"T2", "T3", "T4"), class = "factor"), k = c(5.041791145, 5.408503999,
5.535282299, 5.346402317, 5.376649977, 5.072021484, 6.097412109,
4.390658006, 5.13676712, 4.472827193, 5.381892125, 4.882544582,
4.655393586, 5.435528121, 4.985185185, 4.548431822, 5.041791145,
5.408503999, 5.535282299, 5.346402317, 5.376649977, 5.072021484,
6.097412109, 4.390658006, 5.13676712, 4.472827193, 5.381892125,
4.882544582, 4.655393586, 5.435528121, 4.985185185, 4.548431822,
5.517125816, 4.772205603, 5.928149807, 4.152323266, 4.666037968,
4.638984928, 4.044444444, 4.720296599, 5.315500686, 4.967790359,
3.520804755, 4.722326417, 5.051895044, 4.807450844, 5.096461818,
5.28703008, 5.653368614, 6.357164944, 3.979492188, 3.928861374,
5.632685221, 5.264668498, 5.281464786, 5.387205387, 4.332381668,
5.250388878, 4.580237638, 4.650926114, 5.65951009, 4.401587625,
5.194587481, 4.184813255, 4.44738449, 5.829977261, 4.331985587,
4.827988338, 4.022222222, 3.672891297, 5.148148148, 4.068381688,
5.71922963, 4.566763848, 5.330442907, 2.422536369, 5.346580575,
4.971865289, 5.018922289, 5.513702624, 4.432146456, 5.692296224,
4.738120151, 4.896057489, 5.50365439, 5.249023438, 5.737818961,
4.260276996, 5.242507722, 4.580758017, 5.021888504, 5.013662642,
4.308286338, 5.50840192, 4.732342764, 4.672289386, 5.715557782,
3.827088497, 4.632069971, 4.935541824, 4.008746356, 4.963859809,
4.836806618, 4.46244856, 4.839677641, 4.498269896, 4.88357943,
4.984069185, 4.596844478, 5.196200195, 5.165529005, 14.74622771,
5.397084548, 7.983198678, 5.691090246, 5.707491082, 5.187172012,
6.297376093, 4.647178889, 4.282407407, 4.333496094, 4.773656052,
4.770999725, 4.092207407, 3.917638484, 5.193905817, 3.704833984,
5.571239611, 4.226680384, 3.65230095, 4.78515625, 5.603027344,
4.159218067, 4.719370009, 4.437016946, 4.407713499, 4.284050303,
4.676783265, 4.311689337, 4.540625, 4.864470022, 4.668176455,
5.221193416, 4.997084123, 4.112752873, 5.587217586, 6.045051626,
4.605417744, 4.35030714, 5.185252617, 4.752696927, 4.446670562,
4.268256569, 4.30372087, 4.025205761, 5.696474074, 4.068342788,
3.5212701, 4.544646911, 5.212620027, 5.31978738, 4.879910442,
4.606482493, 4.33502906, 5.294067215, 5.770262391, 4.264308136,
4.501028807, 2.944958848, 4.180638577, 4.120435057, 3.833076111,
4.496793003, 4.232167131, 3.783896334, 5.070553936, 4.825776352,
4.643534043, 6.318587106, 5.66205358, 5.194631597, 4.72557037,
4.195096521, 4.956238551, 3.503093444, 5.24857851, 4.792524005,
4.44229595, 5.285131195, 4.335878892, 4.170953361, 4.045779268
)), row.names = c(NA, -192L), class = "data.frame")
What I tried and the first warning
fit1 <- lmer(k ~ 0 + photo*temp*time_date + (1|tank), data = fem.fish, REML = FALSE)
fixed-effect model matrix is rank deficient so dropping 12 columns / coefficients
boundary (singular) fit: see ?isSingular
My summary and another warning on correlation matrix
summary(fit1)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: k ~ 0 + photo * temp * time_date + (1 | tank)
Data: fem.fish
AIC BIC logLik deviance df.resid
551.2 635.9 -249.6 499.2 166
Scaled residuals:
Min 1Q Median 3Q Max
-2.7467 -0.4380 -0.0447 0.3663 9.7226
Random effects:
Groups Name Variance Std.Dev.
tank (Intercept) 0.0000 0.0000
Residual 0.7883 0.8879
Number of obs: 192, groups: tank, 4
Fixed effects:
Estimate Std. Error t value
photoAmbient 5.284e+00 3.139e-01 16.832
photoCompress 4.937e+00 3.139e-01 15.728
temp7 -1.218e-14 4.439e-01 0.000
time_date17-Jan-19 -9.116e-02 4.439e-01 -0.205
time_date31-Jan-19 -9.798e-02 4.439e-01 -0.221
time_date14-Feb-19 1.264e-01 4.439e-01 0.285
time_date28-Feb-19 -3.986e-01 4.439e-01 -0.898
time_date14-Mar-19 3.655e-01 4.439e-01 0.823
time_date27-Mar-19 -3.979e-01 4.439e-01 -0.896
time_date10-Apr-19 -4.122e-01 4.439e-01 -0.929
time_date24-Apr-19 -2.184e-01 4.439e-01 -0.492
photoCompress:temp7 8.874e-15 6.278e-01 0.000
photoCompress:time_date31-Jan-19 -2.957e-01 6.278e-01 -0.471
photoCompress:time_date28-Feb-19 1.575e+00 6.278e-01 2.509
photoCompress:time_date14-Mar-19 -6.073e-01 6.278e-01 -0.967
temp7:time_date17-Jan-19 -4.121e-02 6.278e-01 -0.066
temp7:time_date31-Jan-19 2.382e-01 6.278e-01 0.379
temp7:time_date14-Feb-19 -2.024e-01 6.278e-01 -0.322
temp7:time_date28-Feb-19 -1.441e+00 6.278e-01 -2.295
temp7:time_date14-Mar-19 -1.104e+00 6.278e-01 -1.759
temp7:time_date27-Mar-19 -4.306e-01 6.278e-01 -0.686
temp7:time_date10-Apr-19 -7.885e-01 6.278e-01 -1.256
temp7:time_date24-Apr-19 -5.872e-01 6.278e-01 -0.935
photoCompress:temp7:time_date14-Mar-19 9.077e-01 8.879e-01 1.022
Correlation matrix not shown by default, as p = 24 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it
fit warnings:
fixed-effect model matrix is rank deficient so dropping 12 columns / coefficients
convergence code: 0
boundary (singular) fit: see ?isSingular
My understanding on t-values is not good at all, so I cannot establish whether there are significant effects or even whether the interactions are significant or not.
I will appreciate your suggestions on the modelling (Fitting the right model?) and more of what you find useful
Thank you so much all.
Upvotes: 0
Views: 1129
Answers (2)
Daniel
Reputation: 7832
I have used your data for the following example. I think due to the fact that all your terms are categorical, you get the rank deficient model. I'd suggest you use time as continuous predictor, thereby you get rid of the rank-deficiency warning.
library(lme4)
library(parameters)
library(performance)
levels(fem.fish$time_date) <- 1:nlevels(fem.fish$time_date)
fem.fish$time_date <- as.numeric(fem.fish$time_date)
fit1 <- lmer(
k ~ 1 + photo * temp * time_date + (1 | tank),
data = fem.fish,
REML = FALSE
)
#> boundary (singular) fit: see ?isSingular
The second warning about the singular fit (now the first, and only warning) is because you literally have no variability in your outcome across the different groups (indicated by tank). This means that the random effects model here gives you not much more benefit than a simple linear model.
ranef(fit1)
#> $tank
#> (Intercept)
#> T1 0
#> T2 0
#> T3 0
#> T4 0
#>
#> with conditional variances for "tank"
Finally, you could use the packages parameters and performance to get comprehensive model summaries (including different p-value approximations like Satterthwaite or Kenward-Roger, standardized parameters or (cluster) robust standard errors) or model fit indices (like r2).
parameters::model_parameters(fit1)
#> Parameter | Coefficient | SE | 95% CI | t | df | p
#> -----------------------------------------------------------------------------------------------------
#> (Intercept) | 5.56 | 0.61 | [ 4.36, 6.76] | 9.06 | 182 | < .001
#> photo [Compress] | -1.46 | 1.04 | [-3.50, 0.58] | -1.41 | 182 | 0.160
#> temp [7] | 0.69 | 0.94 | [-1.16, 2.53] | 0.73 | 182 | 0.467
#> time_date | -0.04 | 0.05 | [-0.13, 0.06] | -0.74 | 182 | 0.461
#> photo [Compress] * temp [7] | 0.73 | 1.52 | [-2.24, 3.70] | 0.48 | 182 | 0.631
#> photo [Compress] * time_date | 0.12 | 0.09 | [-0.06, 0.31] | 1.35 | 182 | 0.178
#> temp [7] * time_date | -0.09 | 0.07 | [-0.23, 0.05] | -1.29 | 182 | 0.198
#> (photo [Compress] * temp [7]) * time_date | -0.07 | 0.13 | [-0.33, 0.19] | -0.52 | 182 | 0.604
performance::r2(fit1)
#> Warning: Can't compute random effect variances. Some variance components equal zero.
#> Solution: Respecify random structure!
#> Random effect variances not available. Returned R2 does not account for random effects.
#> # R2 for Mixed Models
#>
#> Conditional R2: NA
#> Marginal R2: 0.088
Now that your time variable is continuous, you may think about a non-linear relationship of the time trend. You could use the spline package to model this, and ggeffects to get effects plots. This works, of course, for the model with linear time trend as well as other curvilinear time trends.
library(ggeffects)
pr <- ggpredict(fit1, c("time_date", "photo", "temp"))
plot(pr)
library(splines)
fit2 <- lmer(
k ~ 1 + photo * temp * bs(time_date) + (1 | tank),
data = fem.fish,
REML = FALSE
)
#> boundary (singular) fit: see ?isSingular
pr <- ggpredict(fit2, c("time_date [all]", "photo", "temp"))
plot(pr)
Hope that helps!
Upvotes: 0
Pablo Vilas
Reputation: 586
Try to import the "lmerTtest" package.
Before fit your model import this package, in this way you will see the p-value and the "*" of significance:
library("lme4")
library("lmerTest")
Upvotes: 2
Related Questions
- R: Predicting with lmer, y ~ . formula error
- Interpretation of an lmer output
- Prediction in a linear mixed model in R
- Issue with using lmer() for a repeated measures ANOVA
- Error in lmer model and lsmeans output
- Computing p-values for a null random effect model in lm4/lmerTest
- lmer output very low P values
- lmer: predictions on population level trigger an error
- Having issues using the lme4 predict function on my mixed models
- Trying to predict with lmer models