Want more features on Pastebin? Sign Up, it's FREE!
Public Pastes
Guest

Nootropics survey: multilevel analysis

By: a guest on Feb 19th, 2014  |  syntax: R  |  size: 7.96 KB  |  views: 99  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print  |  QR code  |  clone
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. # suggested by http://slatestarcodex.com/2014/02/16/nootropics-survey-results-and-analysis/#comment-40637
  2. library(reshape)
  3. library(lme4)
  4. nootropics <- read.csv("http://dl.dropboxusercontent.com/u/182368464/2014-ssc-nootropicssurvey.csv")
  5.  
  6. # not going to mess with doses or schedules: just want the ratings
  7. 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 ))
  8. # convert each of the '*Experience' columns into a single level of a generic 'Nootropic' factor
  9. long <- reshape(nootropicsResponses, idvar = "Subject", v.names="Response", timevar="Nootropic", times=names(nootropicsResponses), varying=list(names(nootropicsResponses)), direction = "long")
  10. summary(long)
  11.   Nootropic            Response       Subject
  12.  Length:4620        Min.   : 0     Min.   :  1.0
  13.  Class :character   1st Qu.: 3     1st Qu.: 39.0
  14.  Mode  :character   Median : 6     Median : 77.5
  15.                     Mean   : 5     Mean   : 77.5
  16.                     3rd Qu.: 8     3rd Qu.:116.0
  17.                     Max.   :10     Max.   :154.0
  18.                     NA's   :3561
  19.  
  20. # non-nested design: we're examining a crossover design of Subject x Nootropics, not a nested design like Subject(Nootropics).
  21. # Nootropics is a random effect: each level/nootropics gets its own estimate; and likewise, each Subject gets their own estimate
  22. # see http://www.jstatsoft.org/v20/i02/paper "Estimating the Multilevel Rasch Model: With the lme4 Package" for another example
  23. lmr1 <- lmer(Response ~ (1|Subject) + (1|Nootropic), data=long); summary(lmr1)
  24.  
  25. ...REML criterion at convergence: 4875
  26.  
  27. Random effects:
  28.  Groups    Name        Variance Std.Dev.
  29.  Subject   (Intercept) 2.98     1.73
  30.  Nootropic (Intercept) 1.24     1.11
  31.  Residual              4.42     2.10
  32. Number of obs: 1059, groups: Subject, 150; Nootropic, 30
  33.  
  34. Fixed effects:
  35.             Estimate Std. Error t value
  36. (Intercept)    4.995      0.271    18.4
  37.  
  38.  
  39. # do Subjects really vary enough to justify 154 random effects rather than 1 fixed effects?
  40. lmr2 <- lmer(Response ~ Subject + (1|Nootropic), data=long)
  41. # yes:
  42. anova(lmr1, lmr2)
  43.  
  44. ...lmr1: Response ~ (1 | Subject) + (1 | Nootropic)
  45. lmr2: Response ~ Subject + (1 | Nootropic)
  46.      Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
  47. lmr1  4 4882 4902  -2437     4874
  48. lmr2  4 5162 5182  -2577     5154     0      0          1
  49.  
  50. # extract the random effects for both Subjects & Nootropics
  51. rr <- ranef(lmr1, postVar=TRUE)
  52.  
  53. # visualize the two sets of random-effects
  54. qq <- qqmath(rr)
  55. qq$Nootropic
  56. # http://i.imgur.com/PmrHwqA.png
  57. qq$Subject
  58. # http://i.imgur.com/a0nJiOk.png
  59.  
  60. # view nootropics sorted by coefficient size
  61. noots <- rr$Nootropic
  62. noots[order(noots$'(Intercept)'), , drop = FALSE]
  63.                           (Intercept)
  64. GingkoExperience             -1.85454
  65. DMAEExperience               -1.84701
  66. GinsengExperience            -1.13313
  67. CholineExperience            -1.08120
  68. HordenineExperience          -0.90128
  69. VinpocetineExperience        -0.87606
  70. HuperzineExperience          -0.72931
  71. LionapossManeExperience      -0.72782
  72. VitaminDExperience           -0.51186
  73. MCTOilExperience             -0.46230
  74. BacopaExperience             -0.41610
  75. CentrophenoxineExperience    -0.33918
  76. CILTEPExperience             -0.26260
  77. PiracetamExperience          -0.10523
  78. SulbutiamineExperience       -0.02941
  79. OxiracetamExperience          0.08505
  80. AshwagandhaExperience         0.09323
  81. AniracetamExperience          0.10354
  82. NoopeptExperience             0.16937
  83. CreatineExperience            0.29738
  84. PramiracetamExperience        0.36177
  85. RhodiolaExperience            0.44556
  86. PhenylpiracetamExperience     0.50122
  87. AdrafinilExperience           0.59758
  88. TheanineExperience            0.63526
  89. ColuracetamExperience         0.69391
  90. TULIPExperience               0.80629
  91. CaffeineExperience            1.58431
  92. ModafinilExperience           2.39321
  93. ArmodafinilExperience         2.50935
  94.  
  95. # what does the alternate version look like?
  96. rr2 <- ranef(lmr2, postVar=TRUE)
  97. noots2 <- rr2$Nootropic
  98. noots2[order(noots2$'(Intercept)'), , drop = FALSE]
  99.                           (Intercept)
  100. GingkoExperience            -1.399810
  101. GinsengExperience           -1.188626
  102. CholineExperience           -0.978048
  103. DMAEExperience              -0.954214
  104. VitaminDExperience          -0.901817
  105. LionapossManeExperience     -0.662886
  106. VinpocetineExperience       -0.611086
  107. MCTOilExperience            -0.457433
  108. HuperzineExperience         -0.437104
  109. BacopaExperience            -0.320166
  110. HordenineExperience         -0.241617
  111. CentrophenoxineExperience   -0.233822
  112. PiracetamExperience         -0.148546
  113. AshwagandhaExperience       -0.135738
  114. SulbutiamineExperience      -0.015350
  115. CILTEPExperience             0.003327
  116. CreatineExperience           0.012964
  117. TULIPExperience              0.144545
  118. RhodiolaExperience           0.207083
  119. AniracetamExperience         0.326791
  120. PramiracetamExperience       0.354544
  121. OxiracetamExperience         0.374706
  122. NoopeptExperience            0.407759
  123. TheanineExperience           0.515546
  124. PhenylpiracetamExperience    0.568206
  125. AdrafinilExperience          0.594847
  126. ColuracetamExperience        0.778954
  127. CaffeineExperience           1.234636
  128. ArmodafinilExperience        1.456277
  129. ModafinilExperience          1.706078
  130.  
  131. # What did correcting for subject-specific ratings buy us?
  132. # How much did nootropic rankings or estimates change from
  133. # the alternative?
  134. subjectDelta <- rr$Nootropic - rr2$Nootropic
  135. sum(subjectDelta$'(Intercept)')
  136. [1] 1.767e-12
  137. subjectDelta[order(subjectDelta$'(Intercept)'), , drop = FALSE]
  138.                           (Intercept)
  139. DMAEExperience              -0.892796
  140. HordenineExperience         -0.659664
  141. GingkoExperience            -0.454732
  142. HuperzineExperience         -0.292204
  143. OxiracetamExperience        -0.289652
  144. CILTEPExperience            -0.265925
  145. VinpocetineExperience       -0.264976
  146. NoopeptExperience           -0.238387
  147. AniracetamExperience        -0.223253
  148. CentrophenoxineExperience   -0.105358
  149. CholineExperience           -0.103155
  150. BacopaExperience            -0.095930
  151. ColuracetamExperience       -0.085046
  152. PhenylpiracetamExperience   -0.066984
  153. LionapossManeExperience     -0.064930
  154. SulbutiamineExperience      -0.014057
  155. MCTOilExperience            -0.004862
  156. AdrafinilExperience          0.002729
  157. PramiracetamExperience       0.007223
  158. PiracetamExperience          0.043313
  159. GinsengExperience            0.055499
  160. TheanineExperience           0.119715
  161. AshwagandhaExperience        0.228967
  162. RhodiolaExperience           0.238472
  163. CreatineExperience           0.284411
  164. CaffeineExperience           0.349678
  165. VitaminDExperience           0.389960
  166. TULIPExperience              0.661741
  167. ModafinilExperience          0.687129
  168. ArmodafinilExperience        1.053072
  169. # so the biggest difference was DMAE falls by +1 rating
  170. # and Armodafinil looks better by +1 (and Modafinil improves too)
clone this paste RAW Paste Data