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
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 <- read.csv("http://dl.dropboxusercontent.com/u/85192141/data.csv") 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:
- pessimistic -0.38 0.12
- hierarchical -0.37 0.11
- 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