5
$\begingroup$

I wish to analyse a simple lab experiment. I have 8 fish. Four are fed on diet A, and four on diet B. I measure their Nitrogen (N) over 5 time periods (so 5 repeated measures per fish). I wish to know 3 pieces of information:

1) On diet A, does N change with time (i.e. a true slope different from zero)?

2) On diet B, does N change with time (i.e. a true slope different from zero)?

3) Do slopes of A and B differ from one another?

Thus I have run a simple mixed model in R, of the form:

Nitrogen ~ Time * Diet + (1|replicate)

Replicate is used as a random effect to control for repeated mesaures per individual fish.

I have compared this model to a simpler model without an interaction term, using analysis of deviance:

Nitrogen ~ Time + Diet + (1|replicate)

The model with the interaction explained signifcantly more variance and thus the better model.

Based on this interaction model, I have the following output table with the parameter estimates, and I wondered if I can use this table to answer my three questions? (I have calculated confidence intervals for these parameter estimates to determine if effects are real).

Fixed effects:                              
             Estimate  Std. Error   t value   +95% CI     -95% CI
(Intercept)  15.8624     0.2332    68.0300    16.3194     15.4053
time          0.0069     0.0009     7.7600     0.0086      0.0051
dietB        -0.1948     0.3298    -0.5900     0.4515     -0.8411
time:dietB   -0.0066     0.0013    -5.2800    -0.0042     -0.0091

I understand the "Intercept" is value of N when time = 0 for diet A (R software works alphabetically so dietA before dietB), while "dietB" + Intercept is intercept for dietB when time = 0.

However, to answer my 3 questions I need to better understand how to interpret the slope information and how to report it. My current understanding is:

"time" gives me a real slope value (0.0069) to report for diet A, and having calculated confidence intervals, I can see they do not cross zero in this instance meaning diet A is a true slope (i.e. different from a slope of zero). So this is Question 1 answered, I hope?

"time:dietB" tells me that slope of dietB is 0.0069("time") + -0.0066, which gives me a slope value to report for diet B (0.0069 + - 0.0066 = 0.0003). As I currently understand it, the Std. Error for "time:dietB" is associated with this comparative slope value (-0.0066), and not on slope B's real value (0.0003), and consequently, confidence intervals also refer to this comparative value (-0.0066), and thus as confidence intervals do not pass through zero, this means that slope B is different from slope A. So this means Question 3 is answered, I hope?

If this is the case, then how do I answer my Question 2 - that slope B is different from a slope of zero? Can I get confidence intervals to answer this question, and/or do I need them?! Or is this information already contained within this summary table? Or instead, do I need to conduct a subsequent / entirely different analysis?

$\endgroup$
1
$\begingroup$

I think it's always difficult to interpret the size of interaction effects, but it can be done in a different way to make it easier.

If you create two new variables where the first is 0 for diet B and equal to time for diet A, and the other variable is the oppoosite, you assess the effect of the slope for each diet more directly:

TimeA <- ifelse(Diet=="A", Time, 0)
TimeB <- ifelse(Diet=="B", Time, 0)

Nitrogen ~ Diet + TimeA + TimeB + (1|replicate)

You will now get separate estimates for the slope for time A and time B, with standard errors so that you can calculate confidence intervals to report. You will see that the estimates for Diet are the same, and the estimate for TimeA is the same as the estimate for Time in your model. The difference is that in this model, you will get one estimate for the slope of time for Diet B, and it will probably be non-significant in your model.

Note that this assumes that the baseline for Time is 0. If not, the results will be more difficult to interpret.

$\endgroup$
  • $\begingroup$ Thanks. This method usefully provides a standard error for the second slope estimate, that can be convereted into a confidence interval to assess difference in slope from zero. The std error estimate using this method is the same as re-leveling the factor and re-running the model. (All other parameter estimates and their std. error also remian unchanged). $\endgroup$ – user2890989 Oct 11 '15 at 11:28
3
$\begingroup$

Whenever I have interactions I try to use the contrast to investigate them closer. Looking at your data I would say that:

1) On diet A, does N change with time (i.e. a true slope different from zero)?

Yes, lower CI > 0

2) On diet B, does N change with time (i.e. a true slope different from zero)?

No, the slope will be close to 0 with the time interaction.

3) Do slopes of A and B differ from one another?

I would recommend testing the two models against each other. I usually do an ANOVA test, e.g. model1 = time + diet and model2 = time * diet with anova(model1, model2).

Without having your data I've experimented a little with the builtin Orthodont dataset:

Basic models

# The contrast package doesn't seem to work with lme4
# so for the contrast I use nlme
library(nlme)
library(contrast)

# Do a mixed model
fit <- lme(distance ~ age * Sex, data = Orthodont, random = ~ 1 | Subject)
summary(fit)

Gives:

Linear mixed-effects model fit by REML
 Data: Orthodont 
       AIC      BIC    logLik
  445.7572 461.6236 -216.8786

Random effects:
 Formula: ~1 | Subject
        (Intercept) Residual
StdDev:    1.816214 1.386382

Fixed effects: distance ~ age * Sex 
                  Value Std.Error DF   t-value p-value
(Intercept)   16.340625 0.9813122 79 16.651810  0.0000
age            0.784375 0.0775011 79 10.120823  0.0000
SexFemale      1.032102 1.5374208 25  0.671321  0.5082
age:SexFemale -0.304830 0.1214209 79 -2.510520  0.0141
 Correlation: 
              (Intr) age    SexFml
age           -0.869              
SexFemale     -0.638  0.555       
age:SexFemale  0.555 -0.638 -0.869

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.59804400 -0.45461690  0.01578365  0.50244658  3.68620792 

Number of Observations: 108
Number of Groups: 27 

Now look at the contrast - note that the b argument is the reference

contrast(fit, 
         b = list(age = 0, Sex = "Male"),
         a = list(age = 0, Sex = "Female"))

gives:

lme model parameter contrast

 Contrast     S.E.     Lower    Upper    t  df Pr(>|t|)
 1.032102 1.537421 -2.017365 4.081569 0.67 102   0.5035

Now look how they behave through different ages:

contrast(fit, 
         b = list(age = 0:6, Sex = "Male"),
         a = list(age = 0, Sex = "Female"))

gives:

lme model parameter contrast

   Contrast     S.E.     Lower      Upper     t  df Pr(>|t|)
  1.0321023 1.537421 -2.017365  4.0815692  0.67 102   0.5035
  0.2477273 1.495837 -2.719258  3.2147126  0.17 102   0.8688
 -0.5366477 1.457191 -3.426979  2.3536839 -0.37 102   0.7134
 -1.3210227 1.421723 -4.141004  1.4989584 -0.93 102   0.3550
 -2.1053977 1.389676 -4.861814  0.6510188 -1.52 102   0.1329
 -2.8897727 1.361292 -5.589890 -0.1896558 -2.12 102   0.0362
 -3.6741477 1.336805 -6.325693 -1.0226023 -2.75 102   0.0071

Now lets plot this

# One can supply both sexes into the b/a arguments but I find it
# often confusing which is which and I prefer to do it by hand
cntrsts <- 
  lapply(levels(Orthodont$Sex), 
         function(sex) {
           contrast(fit, 
                    b = list(age = seq(0, 10, by=.1), Sex = sex),
                    a = list(age = 0, Sex = "Female"))

         })
names(cntrsts) <- levels(Orthodont$Sex)

# Create the dataset for ggplot
library(magrittr)
cntr_data <- 
  lapply(names(cntrsts),
         function(sex)
         data.frame(
           Age = cntrsts[[sex]]$age,
           Sex = sex,
           Effect = cntrsts[[sex]]$Contrast,
           Lower = cntrsts[[sex]]$Lower,
           Upper = cntrsts[[sex]]$Upper)) %>%
  do.call(rbind, .)

# Do the plot
library(ggplot2)
ggplot(cntr_data, aes(y = Effect, x = Age, 
                      group = Sex, fill = Sex)) +
  facet_grid(~ Sex) +
  geom_ribbon(aes(ymax = Upper, ymin = Lower, color = Sex)) + 
  geom_line() +
  theme_bw()

enter image description here

The ANOVA-test

# I use lme4 for this as the lme gives a warning that
#   "REML comparisons are not meaningful."
library(lme4)
fit_plain <- lmer(distance ~ age + Sex + (1 | Subject), data = Orthodont)
fit_int <- update(fit_plain, .~.+age:Sex)
# Note that this gives a warning
anova(fit_plain,
      fit_int)

gives:

Data: Orthodont
Models:
fit_plain: distance ~ age + Sex + (1 | Subject)
fit_int: distance ~ age + Sex + (1 | Subject) + age:Sex
          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit_plain  5 444.86 458.27 -217.43   434.86                           
fit_int    6 440.64 456.73 -214.32   428.64 6.2174      1    0.01265 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
$\endgroup$
  • $\begingroup$ Thanks for this, it provides some helpful insight into how to think about and approach this problem. I had already tried the anova model comparison, but investigating the contrast in lower and upper parameter estimates is revealing. $\endgroup$ – user2890989 Oct 11 '15 at 11:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.