### Load lme4 package
require("lme4");
### Load data
dat.long <- read.table("http://sciencerep.org/files/7/the%20cannabis%20show%20-%20data%20in%20long%20format.tsv", header=TRUE, sep = "\t");
### Set 'participant' as factor
dat.long$participant <- factor(dat.long$id);
head(dat.long);
###### Trying to analyse this long dataset properly, specifying that:
### 1) every datapoint is a measurement of a dichotomous
### measure of behavior (usedCannabis_bi);
### 2) there are two measurements per participant (before and after the
#### intervention), which are nested within participants (participant);
### 3) participants are nested within schools (school)
### 4) schools were either control of intervention schools (cannabisShow)
###
### The research question is whether the intervention causes a decrease
### in cannabis use over time (i.e. whether at the second moment, cannabis
### use is lower in the individuals who had the intervention).
###
### Note that in the lmer syntax, the model consists of terms separated by,
### plusses (+) where a term is either 1 (the intercept; implicitly included,
### so needs to be set at 0 to force it out of the model), a variable (or
### interaction between variables), or a random effect term. See the
### explanation at http://glmm.wikidot.com/faq or the book of Douglas Bates
### at http://lme4.r-forge.r-project.org/book/ for more details.
###
### In model 1.null, cannabis use is a function of the intercept (cannabis use
### at baseline) and moment (cannabis use 1 week after the intervention). Both
### the intercept and the effect of moment vary per participant, where
### participants are nested within schools.
###
### In model 1.model, cannabis use is a function of the intercept, moment,
### whether participants received the intervention, and the interaction of
### those last two.
rep_measures.1.null <- lmer(formula = usedCannabis_bi ~
1 + moment + (1 + moment|school/participant),
family=binomial(link = "logit"), data=dat.long);
rep_measures.1.model <- update(rep_measures.1.null, .~. + moment*cannabisShow);
rep_measures.1.null;
rep_measures.1.model;
anova(rep_measures.1.null, rep_measures.1.model);
### This is odd. The full model fits a lot better than the null model.
### However, none of the added terms significantly contributes to the
### prediction of cannabis use . . . ?
### When we force fixed intercepts and slopes for every participant (so that
### intercepts and slopes can only vary per school), suddenly the effect of
### our intervention becomes significant, but the difference in fit between
### the null model and the full model becomes much smaller.
rep_measures.2.null <- lmer(formula = usedCannabis_bi ~
1 + moment + (1 + moment|school),
family=binomial(link = "logit"), data=dat.long);
rep_measures.2.model <- update(rep_measures.2.null, .~. + moment*cannabisShow);
rep_measures.2.null;
rep_measures.2.model;
anova(rep_measures.2.null, rep_measures.2.model);