1
$\begingroup$

I have a question regarding repeated measures and GLMs:

Suppose i have counted the abundance of some species in lakes at different time points - every lake received a different treatment at time = 0.

My data looks like this:

df <- structure(list(y = c(1, 3, 5, 1, 4, 1, 4, 1, 1, 0, 5, 2, 3, 2, 
                           3, 2, 4, 4, 3, 2, 1, 3, 8, 1, 5, 4, 6, 3, 5, 0, 1, 2, 0, 2, 6, 
                           1, 7, 3, 3, 2, 11, 0, 0, 1, 0, 1, 3, 0, 10, 6, 6, 2, 9, 0, 0, 
                           2, 0, 0, 3, 1, 10, 7, 4, 3, 12, 0, 0, 1, 0, 2, 4, 0, 8, 5, 3, 
                           4, 8, 1, 3, 5, 0, 5, 4, 2, 3, 4, 4, 2, 7, 1, 8, 4, 3, 7, 5, 7, 
                           4, 7, 3, 4, 7, 2, 7, 5, 3, 3, 6, 12, 7, 7, 1, 5, 20, 4, 10, 4, 
                           3, 4, 14, 15, 4, 7, 3, 2, 14, 1, 8, 8, 1, 3, 9, 15), 
                     time = structure(c(1L, 
                             1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
                             2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
                             3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
                             5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 
                             6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
                             7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 
                             9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 
                             10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 
                             11L, 11L, 11L, 11L, 11L, 11L, 11L), 
                                      .Label = c("-4", "-1", "0.1", "1", "2", "4", "8", "12", 
                                                 "15", "19", "24"), class = "factor"), 
                     treatment = structure(c(2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 
                                             3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 
                                             2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 
                                             3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 
                                             4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 
                                             3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 
                                             2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 
                                             3L, 1L, 5L, 4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 
                                             4L, 2L, 5L, 3L, 1L, 4L, 2L, 1L, 1L, 3L, 1L, 5L, 4L, 2L, 5L, 
                                             3L, 1L, 4L), 
                                           .Label = c("0", "0.1", "0.9", "6", "44"), 
                                           class = "factor"), 
                     plots = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                                         11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
                                         1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 
                                         3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 
                                         5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 
                                         7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
                                         9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                                         11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
                                         1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 
                                         3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L), 
                                       .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12"), 
                                       class = "factor")), 
                .Names = c("y","time", "treatment", "plots"), 
                row.names = c(NA, -132L), class = "data.frame")
str(df)
# y : Counts
# time : sampling time
# treatment: treatment applied 
# plots: every plot/lake forms a series, treatment has been applied to different plots

And here is the time-course of the counts plotted - the lines are just a smoother... enter image description here There seems to be some interaction between treatment and time (counts drop downafter t = 0, but then recover).

I am mainly interested in the treatment and treatment:time interaction. The time effect is of minor interest - since it is known that there is time course over the year...

I could fit a negative-binomial GLM to this data:

require(MASS)
# negbin GLM with interaction
mod_nb <- glm.nb(y ~ time * treatment, data = df)
(mod_nb_aov <- anova(mod_nb))

Q1: Does this take the repeated measure of the same lakes into account? (since i have time a fixed factor in the model)

Since I am not sure about Q1, I thought that I could use restricted permutations to take this into account (permuted lakes, keeping timely neighbored observations together).

This could be done with the permute-package quite easily - something along these lines:

require(permute)
# permute complete time-series, but not within series
control <- how(within = Within(type = 'none'),
               plots = Plots(strata = df$plots, type = 'free'),
               nperm=200)
permutations <- shuffleSet(nrow(df), control = control)

out <- NULL
for(i in 1:nrow(permutations)){
  df_perm <- df[permutations[i, ] , c('time', 'treatment')]
  out[[i]] <-  glm.nb(df$y ~ time * treatment, data = df_perm)
}

Q2-x: But the I wondered if this is not redundant? - Keeping time and restricting permutations. Maybe I should fit a model y ~ time:treatment + treatment?

I know that this could be done using mixed models, however is there also a way via restricted permutations?

I hope I have clearly described the problem... Any thoughts are appreciated.

$\endgroup$
2
$\begingroup$

Since you said in chat that you can't use a GLMM, I have no idea how to deal with the repeated measures. However, I suggest to use something like a cosinor model:

df$time <- as.numeric(as.character(df$time))

library(MASS)
mod_nb <- glm.nb(y ~ time * treatment, data = df)
period=28
mod_nb1 <- glm.nb(y ~ time * cos(time/period*2*pi)*sin(time/period*2*pi)*treatment, data = df)
#Warning message:
#  glm.fit: fitted rates numerically 0 occurred 

AIC(mod_nb, mod_nb1)
#        df      AIC
#mod_nb  11 575.2750
#mod_nb1 41 561.6936
summary(mod_nb1)

set.seed(42)
df_pred <- stack(simulate(mod_nb1, nsim=10000))
df_pred$time <- df$time
df_pred$treatment <- df$treatment

library(plyr)
df_pred <- ddply(df_pred, .(time, treatment), summarize, y=mean(values))


library(ggplot2)
ggplot(df, aes(x=time, y=y, colour=treatment, group=treatment)) + geom_point() +
  geom_line(data=df_pred)

enter image description here

Of course, you'd need to test if you should simplify the model.

$\endgroup$
1
$\begingroup$

In answer to Q1, your approach does not take into account any issue arising from the repeated measures. Not sure about the permutation approach to deal with repeated measures, but you might look at the bootstrap (either parametric or non-parametric) to address any potential issue arising from the repeated measures. The basic approach (using the boot package) is to set up a resampling scheme which keeps the repeated measures structure in your data, fit the model n times to n samples of your data and assess your statistics of interest (coefficients, their s.e., AIC, Deviance etc...). What is your reason for not using glmms?

$\endgroup$
  • $\begingroup$ jupp, the permutation approach is similar to the bootstrap approach - keeps the repeated measures together. $\endgroup$ – EDi Dec 13 '13 at 11:58

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.