# Mixed/hierarchical/multilevel model for n-of-1/single-subject design: incorporating baseline data?

So I often do little self-experiments where I blind & randomize things; these can be formulated as your normal t-tests, but sometimes the measured metrics have extensive baselines which seem like they could be used for more accurate answers. A bunch of reading upon n-of-1 and single-subject designs suggested that people have been moving to mixed/hierarchical/multilevel models for analyzing such setups (eg. Nelson 2012 "Hierarchical linear modeling versus visual analysis of single subject design data" or "The Aggregation of Single-Case Results using Hierarchical Linear Models").

As I understand it, the idea is to split the subject's data into experiment vs baseline, and treat those as the groups. I'm trying to understand how sensible this is with a recent experiment, so hopefully someone can point out if I go wrong in using lmer here.

We start with a regular linear model which examines purely the experimental data (the numeric Response vs the binary Intervention variables) and ignores the extensive baseline phase before, during, and after the experiment:

R> experiment <- read.csv("http://dl.dropboxusercontent.com/u/85192141/data.csv")
R> summary(lm(Response ~ Intervention, data=experiment))

...
Residuals:
Min      1Q  Median      3Q     Max
-1.0156 -0.8889 -0.0156  0.1111  1.1111

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    3.0156     0.0889    33.9   <2e-16
Intervention  -0.1267     0.1262    -1.0     0.32

Residual standard error: 0.711 on 125 degrees of freedom
(145 observations deleted due to missingness)
Multiple R-squared:  0.008, Adjusted R-squared:  6.73e-05
F-statistic: 1.01 on 1 and 125 DF,  p-value: 0.317

R> confint(lm(Response ~ Intervention, data=experiment))
2.5 % 97.5 %
(Intercept)   2.8397  3.192
Intervention -0.3765  0.123


The estimated coefficient is not statistically-significant: -0.38-0.12. But it's definitely slanted towards being negative. So this is the 'conservative' case, where we ignore the baseline entirely. What's the optimistic case? Well, it seems to me that the optimistic case is when we take the entire baseline and assume it is exactly the same as the 'off'/0 intervention in the experiment, in which case we get a narrower CI (because our estimate of the intercept has halved its standard error):

R> experiment$Intervention[is.na(experiment$Intervention)] <- 0
R> summary(lm(Response ~ Intervention, data=experiment))

...
Residuals:
Min      1Q  Median      3Q     Max
-1.9924 -0.8889  0.0076  1.0076  1.1111

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    2.9924     0.0375   79.88   <2e-16
Intervention  -0.1036     0.1012   -1.02     0.31

Residual standard error: 0.746 on 458 degrees of freedom
Multiple R-squared:  0.00228,   Adjusted R-squared:  0.000101
F-statistic: 1.05 on 1 and 458 DF,  p-value: 0.307

R> confint(lm(Response ~ Intervention, data=experiment))
2.5 %  97.5 %
(Intercept)   2.9188 3.06607
Intervention -0.3025 0.09538


It's narrowed to -0.30-0.10; still not statistically-significant, but closer.

It seems to me that a hierarchical model ought to produce a CI intermediate between the pessimistic and optimistic cases: it loses some power because it's estimating how different the two phases are before it does any combining.

Here is my multilevel model, split between baseline and experimental phases:

library(lme4)
experiment$Phase <- ifelse(is.na(experiment$Intervention), TRUE, FALSE)
model <- lmer(Response ~ Intervention + (1|Phase), data=experiment); summary(model)

...
AIC BIC logLik deviance REMLdev
286 297   -139      273     278
Random effects:
Groups   Name        Variance Std.Dev.
Phase    (Intercept) 0.0106   0.103
Residual             0.5057   0.711
Number of obs: 127, groups: Phase, 1

Fixed effects:
Estimate Std. Error t value
(Intercept)     3.016      0.136    22.2
Intervention   -0.127      0.126    -1.0

Correlation of Fixed Effects:
(Intr)
Interventin -0.461

m <- mcmcsamp((lmer(Response ~ Intervention + (1|Phase), data=experiment)), n = 100000)
HPDinterval(m, prob=0.95)$fixef lower upper (Intercept) -45.3107 56.6558 Intervention -0.3742 0.1191  The estimated CI comes out exactly in the middle, as expected: 1. pessimistic -0.38 0.12 2. hierarchical -0.37 0.11 3. optimistic -0.30 0.10 So, my basic question is: is this a sane approach to take? It's spitting out answers that seem intuitively correct, but that might just be a coincidence. Incidentally, one might be worried about time trends. The randomization/blocking would fix that in the experimental period but not the baseline. Fortunately, that doesn't seem to be an issue: experiment$Time <- 1:nrow(experiment)
summary(lmer(Response ~ Intervention + Time + (1|Phase), data=experiment))

...
AIC BIC logLik deviance REMLdev
298 312   -144      272     288
Random effects:
Groups   Name        Variance Std.Dev.
Phase    (Intercept) 0.0106   0.103
Residual             0.5055   0.711
Number of obs: 127, groups: Phase, 1

Fixed effects:
Estimate Std. Error t value
(Intercept)   3.42517    0.42325    8.09
Intervention -0.12398    0.12621   -0.98
Time         -0.00132    0.00129   -1.02

Correlation of Fixed Effects:
(Intr) Intrvn
Interventin -0.128
Time        -0.947 -0.021

• A couple comments. (1) The results of your current model are not sensible. Notice in the model output that all of the data are contained in a single "group", namely the Phase==F group. What is happening is that all the rows with missing data on the Intervention variable are still being dropped. (2) It is not clear that this is a multilevel problem in the first place. Although you have two "groups" of observations (baseline and non-baseline), considering groups to be a random factor is statistically and conceptually dubious, to say the least. – Jake Westfall May 28 '13 at 5:59