# I could not fit a linear mixed effects model with time as a random effect, but the time effect should influence the observed results

I have data which has an independent variable (temperature), dependent variable (gasflux) and an additional time variable as given below.

y= as.data.frame(structure(
list(time = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2,  3, 4, 5, 6, 7, 8, 9, 10),
dependent = c(10, 5, 2.5, 1.3, 0.6,  0.3, 0.2, 0.1, 0.05, 0.025, 14, 7,
3.5, 1.8, 0.9, 0.4, 0.2, 0.1, 0.05, 0.025, 16, 8, 4, 2, 1,
0.5, 0.3, 0.1, 0.05, 0.025, 20, 10, 5, 2.5, 1.3, 0.6, 0.3,
0.2, 0.1, 0.05),
independent = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4)),
.Names = c("time", "dependent", "independent"),
class = c("tbl_df", "tbl", "data.frame"),
row.names = c(NA, -40L),
spec = structure(list(cols = structure(list(
time = structure(list(), class = c("collector_double", "collector")),
dependent = structure(list(), class = c("collector_double", "collector")),
independent = structure(list(), class = c("collector_double", "collector"))),
.Names = c("time", "dependent", "independent")),
default = structure(list(), class = c("collector_guess", "collector"))),
.Names = c("cols", "default"), class = "col_spec")))


The data look like this:

My objective is to find the effect of the independent variable on the dependent variable, but the effect changes with time. I tried linear regression of independent variable with the dependent variable, it didn't give significant effect.

I wanted to try a linear mixed effect, but I couldn't fit time as a random effect. Any suggestions?

• Why do you want time as a random effect? – conjugateprior Feb 27 '19 at 21:09
• I thought that could help me give significant effect of independent variable in dependent variable. – Kathiravan Meeran Feb 27 '19 at 21:55
• @conjugateprior is there any other statistical method or test which show me the relation between independent and dependent variable relation? – Kathiravan Meeran Feb 27 '19 at 21:58
• Your data are not linear. You need to linearize them with a link function (e.g., log) ...you likely also have auto correlation and may need an ar-1 covariance structure, or use non-linear regression...exponential decay function should fit. – OliverFishCode Feb 28 '19 at 4:57
• In this case it obvious there is an effect, but in general you should not go looking to produce a significant effect. let the data drive the outcome, don't steer results toward what you desire. – OliverFishCode Feb 28 '19 at 5:01

library(glmmTMB)
library(DHARMa)
options(scipen = 999)


Using GLMMTMB (Generalized Linear Mixed Models based on Template model builder) using the Gamma series of continuous distribution and a log link, the independent variable is statistically significant. Your data while continuous were likely not normally distributed, thus the use of the Gamma distribution. I have also evaluated the residuals using the DHARMa package. Ignore the right panel in the residual plot as it is an effect of the GLMMTMB package, but based on the left panel the residuals are good to go.

Code:

y = as.data.frame(structure(list(time = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8,
9, 10), dependent = c(10, 5, 2.5, 1.3, 0.6, 0.3, 0.2, 0.1, 0.05, 0.025, 14, 7,
3.5, 1.8, 0.9, 0.4, 0.2, 0.1, 0.05, 0.025, 16, 8, 4, 2, 1, 0.5, 0.3, 0.1, 0.05,
0.025, 20, 10, 5, 2.5, 1.3, 0.6, 0.3, 0.2, 0.1, 0.05), independent = c(0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4)), .Names = c("time", "dependent",
"independent"), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -40L),
spec = structure(list(cols = structure(list(time = structure(list(), class =
c("collector_double",
"collector")), dependent = structure(list(), class = c("collector_double",
"collector")), independent = structure(list(), class = c("collector_double",
"collector"))), .Names = c("time", "dependent", "independent")), default = structure(list(),
class = c("collector_guess", "collector"))), .Names = c("cols", "default"),
class = "col_spec")))

model = glmmTMB(dependent ~ independent + (1|time),data = y, family = Gamma(link = "log"))
summary(model)

res = residuals(model, type = 'pearson')
resid_check = simulateResiduals(model)# build sim. residuals
plot(resid_check , rank=T)#check residuals
hist(res)


Output:

Family: Gamma  ( log )
Formula:          dependent ~ independent + (1 | time)
Data: y

AIC      BIC   logLik deviance df.resid
1.8      8.5      3.1     -6.2       36

Random effects:

Conditional model:
Groups Name        Variance Std.Dev.
time   (Intercept) 3.865    1.966
Number of obs: 40, groups:  time, 10

Dispersion estimate for Gamma family (sigma^2): 0.0226

Conditional model:
Estimate Std. Error z value            Pr(>|z|)
(Intercept) -0.75449    0.62317  -1.211               0.226
independent  0.16155    0.01649   9.796 <0.0000000000000002


DHARMA residual evaluation plot (Ignore the right panel, a result of using GLMMTMB)

Histogram of Pearson Residuals

Goodluck!

• @KathiravanMeeran I hope this helps – OliverFishCode Feb 28 '19 at 18:04
• Your answer explained everything clearly and it works. Amazing!! – Kathiravan Meeran Mar 1 '19 at 13:40

A simple log transformation seems to describe this data well

m1 <- lm(log(dependent) ~ independent + time, data = y)

summary(m1)\$r.squared # 0.995


Consequently, you might build the log link into your model, as @devon-oliver does in their answer, e.g. like this:

m2 <- glm(dependent ~ independent + time, data = y,

Confidence intervals for predictions are a little trickier here though. The model improves slightly if independent is treated as discrete, but this might be noise, or a sign that the relationship has structure that would be clearer with more variation.