Nootropics survey: multilevel analysis

# suggested by http://slatestarcodex.com/2014/02/16/nootropics-survey-results-and-analysis/#comment-40637
library(reshape)
library(lme4)
nootropics <- read.csv("http://dl.dropboxusercontent.com/u/182368464/2014-ssc-nootropicssurvey.csv")

# not going to mess with doses or schedules: just want the ratings
nootropicsResponses <- with(nootropics, data.frame (AdrafinilExperience=AdrafinilExperience, AniracetamExperience=AniracetamExperience, ArmodafinilExperience=ArmodafinilExperience, AshwagandhaExperience=AshwagandhaExperience, BacopaExperience=BacopaExperience, CILTEPExperience=CILTEPExperience, CaffeineExperience=CaffeineExperience, CentrophenoxineExperience=CentrophenoxineExperience, CholineExperience=CholineExperience, ColuracetamExperience=ColuracetamExperience, CreatineExperience=CreatineExperience, DMAEExperience=DMAEExperience, GingkoExperience=GingkoExperience, GinsengExperience=GinsengExperience, HordenineExperience=HordenineExperience, HuperzineExperience=HuperzineExperience, LionapossManeExperience=LionapossManeExperience, MCTOilExperience=MCTOilExperience, ModafinilExperience=ModafinilExperience, NoopeptExperience=NoopeptExperience, OxiracetamExperience=OxiracetamExperience, PhenylpiracetamExperience=PhenylpiracetamExperience, PiracetamExperience=PiracetamExperience, PramiracetamExperience=PramiracetamExperience, RhodiolaExperience=RhodiolaExperience, SulbutiamineExperience=SulbutiamineExperience, TULIPExperience=TULIPExperience, TheanineExperience=TheanineExperience, VinpocetineExperience=VinpocetineExperience, VitaminDExperience=VitaminDExperience ))
# convert each of the '*Experience' columns into a single level of a generic 'Nootropic' factor
long <- reshape(nootropicsResponses, idvar = "Subject", v.names="Response", timevar="Nootropic", times=names(nootropicsResponses), varying=list(names(nootropicsResponses)), direction = "long")
summary(long)
  Nootropic            Response       Subject
 Length:4620        Min.   : 0     Min.   :  1.0
 Class :character   1st Qu.: 3     1st Qu.: 39.0
 Mode  :character   Median : 6     Median : 77.5
                    Mean   : 5     Mean   : 77.5
                    3rd Qu.: 8     3rd Qu.:116.0
                    Max.   :10     Max.   :154.0
                    NA's   :3561

# non-nested design: we're examining a crossover design of Subject x Nootropics, not a nested design like Subject(Nootropics).
# Nootropics is a random effect: each level/nootropics gets its own estimate; and likewise, each Subject gets their own estimate
# see http://www.jstatsoft.org/v20/i02/paper "Estimating the Multilevel Rasch Model: With the lme4 Package" for another example
lmr1 <- lmer(Response ~ (1|Subject) + (1|Nootropic), data=long); summary(lmr1)

...REML criterion at convergence: 4875

Random effects:
 Groups    Name        Variance Std.Dev.
 Subject   (Intercept) 2.98     1.73
 Nootropic (Intercept) 1.24     1.11
 Residual              4.42     2.10
Number of obs: 1059, groups: Subject, 150; Nootropic, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept)    4.995      0.271    18.4


# do Subjects really vary enough to justify 154 random effects rather than 1 fixed effects?
lmr2 <- lmer(Response ~ Subject + (1|Nootropic), data=long)
# yes:
anova(lmr1, lmr2)

...lmr1: Response ~ (1 | Subject) + (1 | Nootropic)
lmr2: Response ~ Subject + (1 | Nootropic)
     Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
lmr1  4 4882 4902  -2437     4874
lmr2  4 5162 5182  -2577     5154     0      0          1

# extract the random effects for both Subjects & Nootropics
rr <- ranef(lmr1, postVar=TRUE)

# visualize the two sets of random-effects
qq <- qqmath(rr)
qq$Nootropic
# http://i.imgur.com/PmrHwqA.png
qq$Subject
# http://i.imgur.com/a0nJiOk.png

# view nootropics sorted by coefficient size
noots <- rr$Nootropic
noots[order(noots$'(Intercept)'), , drop = FALSE]
                          (Intercept)
GingkoExperience             -1.85454
DMAEExperience               -1.84701
GinsengExperience            -1.13313
CholineExperience            -1.08120
HordenineExperience          -0.90128
VinpocetineExperience        -0.87606
HuperzineExperience          -0.72931
LionapossManeExperience      -0.72782
VitaminDExperience           -0.51186
MCTOilExperience             -0.46230
BacopaExperience             -0.41610
CentrophenoxineExperience    -0.33918
CILTEPExperience             -0.26260
PiracetamExperience          -0.10523
SulbutiamineExperience       -0.02941
OxiracetamExperience          0.08505
AshwagandhaExperience         0.09323
AniracetamExperience          0.10354
NoopeptExperience             0.16937
CreatineExperience            0.29738
PramiracetamExperience        0.36177
RhodiolaExperience            0.44556
PhenylpiracetamExperience     0.50122
AdrafinilExperience           0.59758
TheanineExperience            0.63526
ColuracetamExperience         0.69391
TULIPExperience               0.80629
CaffeineExperience            1.58431
ModafinilExperience           2.39321
ArmodafinilExperience         2.50935

# what does the alternate version look like?
rr2 <- ranef(lmr2, postVar=TRUE)
noots2 <- rr2$Nootropic
noots2[order(noots2$'(Intercept)'), , drop = FALSE]
                          (Intercept)
GingkoExperience            -1.399810
GinsengExperience           -1.188626
CholineExperience           -0.978048
DMAEExperience              -0.954214
VitaminDExperience          -0.901817
LionapossManeExperience     -0.662886
VinpocetineExperience       -0.611086
MCTOilExperience            -0.457433
HuperzineExperience         -0.437104
BacopaExperience            -0.320166
HordenineExperience         -0.241617
CentrophenoxineExperience   -0.233822
PiracetamExperience         -0.148546
AshwagandhaExperience       -0.135738
SulbutiamineExperience      -0.015350
CILTEPExperience             0.003327
CreatineExperience           0.012964
TULIPExperience              0.144545
RhodiolaExperience           0.207083
AniracetamExperience         0.326791
PramiracetamExperience       0.354544
OxiracetamExperience         0.374706
NoopeptExperience            0.407759
TheanineExperience           0.515546
PhenylpiracetamExperience    0.568206
AdrafinilExperience          0.594847
ColuracetamExperience        0.778954
CaffeineExperience           1.234636
ArmodafinilExperience        1.456277
ModafinilExperience          1.706078

# What did correcting for subject-specific ratings buy us?
# How much did nootropic rankings or estimates change from
# the alternative?
subjectDelta <- rr$Nootropic - rr2$Nootropic
sum(subjectDelta$'(Intercept)')
[1] 1.767e-12
subjectDelta[order(subjectDelta$'(Intercept)'), , drop = FALSE]
                          (Intercept)
DMAEExperience              -0.892796
HordenineExperience         -0.659664
GingkoExperience            -0.454732
HuperzineExperience         -0.292204
OxiracetamExperience        -0.289652
CILTEPExperience            -0.265925
VinpocetineExperience       -0.264976
NoopeptExperience           -0.238387
AniracetamExperience        -0.223253
CentrophenoxineExperience   -0.105358
CholineExperience           -0.103155
BacopaExperience            -0.095930
ColuracetamExperience       -0.085046
PhenylpiracetamExperience   -0.066984
LionapossManeExperience     -0.064930
SulbutiamineExperience      -0.014057
MCTOilExperience            -0.004862
AdrafinilExperience          0.002729
PramiracetamExperience       0.007223
PiracetamExperience          0.043313
GinsengExperience            0.055499
TheanineExperience           0.119715
AshwagandhaExperience        0.228967
RhodiolaExperience           0.238472
CreatineExperience           0.284411
CaffeineExperience           0.349678
VitaminDExperience           0.389960
TULIPExperience              0.661741
ModafinilExperience          0.687129
ArmodafinilExperience        1.053072
# so the biggest difference was DMAE falls by +1 rating
# and Armodafinil looks better by +1 (and Modafinil improves too)