attachment:LSA13-Lecture2-GLMM.R of LSA2013Regression

   1 ## @knitr Options, include=F, cache=F, message=F
   2 opts_chunk$set(size='footnotesize', eval=T, cache=T, prompt=F, message=F, tidy=F, fig.width=4, fig.height=4)
   3 
   4 
   5 ## @knitr LinearMixedModelSimulation1, echo=T, results='markup', cache=T
   6 ## simulate some data
   7 sigma.b <- 125       # inter-subject variation larger than
   8 sigma.e <- 40        # intra-subject, inter-trial variation
   9 alpha <- 500
  10 beta <- 12
  11 M <- 6                            # number of participants
  12 n <- 50                           # trials per participant
  13 b <- rnorm(M, 0, sigma.b)         # individual differences
  14 nneighbors <- rpois(M*n,3) + 1    # generate num. neighbors
  15 subj <- rep(1:M,n)
  16 RT <- alpha + beta * nneighbors + # simulate RTs!
  17   b[subj] + rnorm(M*n,0,sigma.e)  #
  18 
  19 
  20 ## @knitr LinearMixedModel1, echo=T, results='markup', cache=T
  21 library(languageR)
  22 library(lme4)
  23 data(lexdec)
  24 l1 = lmer(RT ~ Frequency + (1 | Subject), data=lexdec)
  25 
  26 
  27 ## @knitr LinearMixedModel2, echo=F, results='markup', cache=T, dependson='LinearMixedModel1'
  28 l1
  29 
  30 
  31 ## @knitr LinearMixedModelMCMC, echo=T, results='markup', cache=T, dependson='LinearMixedModel1'
  32 library(languageR)
  33 p = pvals.fnc(l1, nsim = 10000, addPlot=F)
  34 print(p)
  35 
  36 
  37 ## @knitr LinearMixedModelR2, echo=T, results='markup', dependson='LinearMixedModel1'
  38 cor(fitted(l1), lexdec$RT)^2
  39 
  40 
  41 ## @knitr LinearMixedModelMCMC2, echo=F, results='hide', cache=T, dependson='LinearMixedModel1', fig.width=3, fig.height=3
  42 library(languageR)
  43 pvals.fnc(l1, nsim = 10000)
  44 
  45 
  46 ## @knitr Ranef, echo=T, results='markup', dependson='LinearMixedModel1'
  47 mean(ranef(l1)$Subject[[1]])
  48 head(ranef(l1)$Subject[1], 5)
  49 
  50 
  51 ## @knitr RanefPredictedVsObserved, echo=F, results='markup', fig.width=3, fig.height=3, dependson='LinearMixedModel1'
  52 # un-log RTs
  53 par(cex=.7)
  54 o = exp(as.vector(tapply(lexdec$RT, lexdec$Subject, mean)))
  55 plot(o, type = "n", xlab = "Subject Index", 
  56   ylab = "Response latency in msecs", 
  57 	main = "Subject as random effect")
  58 # place by-subject mean on plot
  59 text(o, as.character(unique(lexdec$Subject)), col = "blue")
  60 # add overall mean
  61 abline(mean(exp(lexdec$RT)),0, col = "gray", lwd=2)
  62 
  63 # un-log model's predictions
  64 p = as.numeric(unlist(exp(coef(l1)$Subject[1] + mean(lexdec$Frequency)*coef(l1)$Subject[2])))
  65 # place model's predictions on plot
  66 text(p, as.character(unique(lexdec$Subject)), col = "red")
  67 legend(x=2, y=850, 
  68        legend=c("Predicted", "Observed"), 
  69        col=c("red", "blue"), pch=1
  70 )
  71 
  72 
  73 ## @knitr LinearMixedModel3, echo=T, results='markup'
  74 l2 = lmer(RT ~ 1 + (1 | Subject) + (1 | Word), 
  75           data = lexdec)
  76 head(ranef(l2)$Subject, 5)
  77 head(ranef(l2)$Word, 5)
  78 
  79 
  80 ## @knitr RanefPredictedVsObserved2, echo=F, results='markup', fig.width=5, fig.height=3.8
  81 # un-log RTs
  82 par(cex=.7)
  83 s = exp(as.vector(tapply(lexdec$RT, lexdec$Subject, mean)))
  84 o = exp(as.vector(tapply(lexdec$RT, lexdec$Word, mean)))
  85 plot(o, ylim = c(min(s,o),max(s,o)), type = "n", xlab = "Word Index", 
  86   ylab = "Response latency in msecs", 
  87   main = "Word as random effect")
  88 # place by-subject mean on plot
  89 text(o, as.character(unique(lexdec$Word)), col = "blue")
  90 # add overall mean
  91 abline(mean(exp(lexdec$RT)),0, col = "gray", lwd=2)
  92 
  93 # un-log model's predictions
  94 p = as.numeric(unlist(exp(coef(l2)$Word[1])))
  95 # place model's predictions on plot
  96 segments(x0=1:length(o), y0=o, x1=1:length(p), 
  97   y1=p, col="red", lwd=2)
  98 points(p, pch=ifelse(p > o, "^", ifelse(p < o, "v", ".")), col="red", lwd=2)
  99 legend(x=2, y=850, 
 100        legend=c("Predicted", "Observed"), 
 101        col=c("red", "blue"), pch=1
 102 )
 103 
 104 
 105 ## @knitr RandomSlopes1, echo=T, results='markup'
 106 l3 = lmer(RT ~ 1 + Frequency + (1 + Frequency | Subject) + (1 | Word), 
 107             data = lexdec)
 108 
 109 
 110 ## @knitr RandomSlopes1b, echo=T, results='markup'
 111 print(l3, corr=F)
 112 
 113 
 114 ## @knitr RandomSlopes2, echo=T, results='markup'
 115 head(ranef(l3)$Subject, 5)
 116 head(ranef(l3)$Word, 5)
 117 
 118 
 119 ## @knitr RandomSlopes4, echo=T, eval=FALSE
 120 ## dotplot(ranef(l3))
 121 
 122 
 123 ## @knitr RandomSlopes5, echo=T, eval=FALSE
 124 ## qqmath(ranef(l3))
 125 
 126 
 127 ## @knitr RandomCorrelations1, echo=T, results='markup', fig.width=2, fig.height=2
 128 par(cex=.5)
 129 plot(ranef(l3)$Subject, main="Random Effect Correlation")
 130 
 131 
 132 ## @knitr RandomCorrelations1b, echo=T, results='markup', fig.width=2, fig.height=2
 133 lexdec$cFrequency = lexdec$Frequency - mean(lexdec$Frequency)
 134 l3b = lmer(RT ~ 1 + cFrequency + (1 + cFrequency | Subject),
 135           data= lexdec)
 136 print(l3b, corr=F)
 137 
 138 
 139 ## @knitr RandomCorrelations2, echo=T, results='markup', fig.width=2, fig.height=2
 140 lexdec$HighFrequency = 
 141     ifelse(lexdec$Frequency > median(lexdec$Frequency), 1, 0)
 142 l4 = lmer(RT ~ 1 + HighFrequency + (1 + HighFrequency | Subject),
 143           data= lexdec)
 144 print(l4, corr=F)
 145 
 146 
 147 ## @knitr RandomCorrelations2b, echo=F, results='markup', fig.width=2, fig.height=2
 148 par(cex=.5)
 149 plot(ranef(l4)$Subject, main="Random Effect Correlation")
 150 
 151 
 152 ## @knitr RandomCorrelations3, echo=T, results='markup', fig.width=2, fig.height=2
 153 lexdec$HighFrequency = ifelse(lexdec$Frequency > median(lexdec$Frequency), 
 154           .5, -.5)
 155 l4b = lmer(RT ~ 1 + HighFrequency + (1 + HighFrequency | Subject),
 156           data= lexdec)
 157 print(l4b, corr=F)
 158 par(cex=.5)
 159 plot(ranef(l4b)$Subject, main="Random Effect Correlation")
 160 
 161 
 162 ## @knitr RandomCorrelations3b, echo=F, results='markup', fig.width=2, fig.height=2
 163 par(cex=.5)
 164 plot(ranef(l4b)$Subject, main="Random Effect Correlation")
 165 
 166 
 167 ## @knitr RandomEffectExtraction, echo=T, results='markup'
 168 # use vcov(model) for variance-covariance matrix of fixed effects
 169 # these determiner the standard errors of the fixed effects.
 170 vcov(l4b)
 171 
 172 # get standard errors of fixed effects:
 173 sqrt(diag(vcov(l4b)))
 174 
 175 
 176 ## @knitr RandomEffectExtraction2, echo=T, results='markup'
 177 # use VarCorr(model) for variance-covariance matrix of random effects
 178 # unlike the fixed effect variance-covariance matrix, only the *structure*
 179 # of the VarCorr is known, the values are parameters that are estimated
 180 # when we fit the model. 
 181 VarCorr(l4b)
 182 
 183 # for further details, see 
 184 # https://stat.ethz.ch/pipermail/r-help/2006-July/109308.html
 185 
 186 
 187 ## @knitr RandomEffectExtraction3, echo=T, results='hide', dependson='LinearMixedModel3'
 188 s = mcmcsamp(l2, 50000)
 189 
 190 # a package for output analysis of MCMC
 191 library(coda)
 192 HPDinterval(s)
 193 
 194 
 195 ## @knitr L2repeated, echo=F, results='markup'
 196 l2
 197 
 198 
 199 ## @knitr RandomEffectExtraction4, echo=F, results='markup', dependson='RandomEffectExtraction3'
 200 HPDinterval(s)
 201 
 202 # The ST slot contains estimates that determine the variance-covariance
 203 # matrix of the random effects (these are the parameters of the Cholesky 
 204 # decomposition; see next slide)
 205 # (to understand what this all means, see # http://sciencemeanderthal.
 206 # wordpress.com/2012/06/28/cholesky-decomposition-of-variance-covariance-
 207 # matrices-in-the-classic-twin-study/)
 208 
 209 
 210 ## @knitr RandomEffectExtraction5, echo=T, results='markup', dependson='RandomEffectExtraction3'
 211 HPDinterval(VarCorr(s, type = "varcov"))
 212 
 213 # (given in the order they appear in l2: variance for Word intercept, 
 214 # Subject intercept, Residual)
Attachment 'LSA13-Lecture2-GLMM.R'

Attached Files
