attachment:LSA13-Lecture1-GLM.R of LSA2013Regression

   1 ## @knitr Options, include=F, cache=F, message=F
   2 opts_chunk$set(size='footnotesize', eval=T, prompt=F, message=F, tidy=F, fig.width=4, fig.height=4)
   3 
   4 
   5 ## @knitr stargazerexamples, echo=F, results='asis', cache=T
   6 library(languageR)
   7 library(stargazer)
   8 data(lexdec)
   9 lexdec$CorrectResponse = with(lexdec, Correct == "correct")
  10 g1a = lm(RT ~ Frequency * NativeLanguage + Trial, data=lexdec)
  11 g1b = lm(RT ~ Frequency * NativeLanguage, data=lexdec)
  12 g2 = glm(CorrectResponse ~ Frequency * NativeLanguage, data=lexdec, family=binomial)
  13 stargazer(g1a,g1b,g2,
  14   align=TRUE, 
  15   dep.var.labels=c("(logged) RT\n","Correct response?"),
  16   covariate.labels=c("Intercept", "Word frequency (logged)","Native language","Trial position","Word frequency (logged):Native language"),
  17   omit.stat=c("LL","ser","f"),
  18   intercept.top=T, 
  19   title="Example Stargazer table generated from R"
  20 )
  21 
  22 
  23 ## @knitr Preliminaries, echo=T, results='markup', cache=T, message=F
  24 version
  25 ls()
  26 
  27 
  28 ## @knitr linearcombination, echo=F, results='markup', cache=T
  29 d = read.csv("scripts/discourse-info.tab", sep="\t", header=T)
  30 d = d[sample(1:nrow(d), 8),]
  31 print(d)
  32 
  33 
  34 ## @knitr linearcombination2, echo=F, results='markup', cache=T
  35 print(d)
  36 
  37 
  38 ## @knitr LoadData1, echo=F, results='markup', cache=T, message=F
  39 library(languageR)
  40 data(lexdec)
  41 # randomly sampling rows from lexdec
  42 head(lexdec[sample(x=1:nrow(lexdec), size=20, replace=F),c(1:5, 9:10)])
  43 
  44 
  45 ## @knitr lexdec, echo=T, eval=F, cache=T, message=F
  46 library(languageR)
  47 data(lexdec)
  48 
  49 # Try some of the following:
  50 help(lexdec)      # learn about this R object
  51 ?lexdec           # the same
  52 class(lexdec)     # what type of R object is lexdec?
  53 nrow(lexdec)      # number of rows (cases)
  54 str(lexdec)       # works on almost all R objects
  55 summary(lexdec)   # summary of each variable in the data.frame
  56 summary(lexdec$RT)
  57 mean(lexdec$RT)
  58 var(lexdec$RT)
  59 with(lexdec, mean(RT))
  60 
  61 # create a new data.frame that is a subset of lexdec
  62 new = subset(lexdec, RT > 7)   
  63 nrow(new)
  64 summary(new$RT)
  65 
  66 
  67 ## @knitr IllustrateData1, echo=T, results='markup', cache=T, message=F
  68 summary(lexdec[,c('RT','Frequency')])
  69 
  70 
  71 ## @knitr LinearModelExample1, echo=T, results='markup', cache=T, message=F
  72 # glm() is R's call to a GLM
  73 # 1, in an R-formula, is a specific symbol for the intercept
  74 # Frequency and RT are variables in the data.frame lexdec
  75 # data tells glm which data.frame to use
  76 # family tells glm which distributions the outcome is assumed to have
  77 m = glm(RT ~ 1 + Frequency, data = lexdec, family = gaussian)
  78 
  79 
  80 ## @knitr LinearModelExample1b, echo=T, results='hide', cache=T, message=F
  81 library(languageR)
  82 data(lexdec)
  83 glm(RT ~ 1 + Frequency, data = lexdec, family = gaussian)
  84 glm(RT ~ Frequency, data = lexdec, family = gaussian)
  85 glm(RT ~ Frequency, data = lexdec)
  86 
  87 
  88 ## @knitr LinearModelExample1c, echo=T, eval=F, cache=T, message=F
  89 glm(RT ~ - 1 + Frequency, data = lexdec)
  90 
  91 
  92 ## @knitr LinearModelExample2, echo=T, results='markup', cache=T, message=F, size='scriptsize'
  93 summary(m)$coefficients[,1:2]
  94 sqrt(summary(m)$dispersion)
  95 
  96 
  97 ## @knitr LinearModelInterceptExample1, echo=T, results='markup', cache=T, message=F
  98 m0 = glm(RT ~ 1, data = lexdec, family = gaussian)
  99 
 100 
 101 ## @knitr LinearModelInterceptExample2, echo=T, results='markup', cache=T, message=F
 102 summary(m0)$coefficients[,1:2]
 103 
 104 
 105 ## @knitr mean, echo=T, results='markup', cache=T, message=F
 106 options(digits=6)
 107 mean(lexdec$RT)
 108 
 109 
 110 ## @knitr LinearModelInterceptVisualization1, echo=T, results='markup', cache=T, fig.height=2, fig.width=3
 111 par(cex=.7, mar=c(4,4,0.1,0.1))
 112 plot(x = lexdec$Frequency, 
 113      y = lexdec$RT, 
 114      ylab = "(log-transformed) RT", 
 115      xlab = "(log-transformed) word frequency", 
 116      type = "n"
 117 )
 118 points(x = lexdec$Frequency, y = lexdec$RT, pch=1, cex=.1, col = "grey")
 119 abline(m0, col = "magenta")
 120 
 121 
 122 ## @knitr LinearModelInterceptVisualization2, echo=F, results='markup', cache=T, fig.height=1.5, fig.width=3
 123 par(cex=.65, mar=c(4,4,0.5,0.1))
 124 plot(x = exp(lexdec$Frequency), 
 125      y = lexdec$RT, 
 126      ylab = "(log-transformed) RT", 
 127      xlab = "word frequency", 
 128      type = "n"
 129 )
 130 points(x = exp(lexdec$Frequency), y = lexdec$RT, pch=1, cex=.1, col = "grey")
 131 o = order(exp(lexdec$Frequency))
 132 lines(x = exp(lexdec$Frequency[o]), y = predict(m0)[o], col = "magenta")
 133 
 134 plot(x = exp(lexdec$Frequency), 
 135      y = exp(lexdec$RT), 
 136      ylab = "RT", 
 137      xlab = "word frequency", 
 138      type = "n"
 139 )
 140 points(x = exp(lexdec$Frequency), y = exp(lexdec$RT), pch=1, cex=.1, col = "grey")
 141 lines(x = exp(lexdec$Frequency[o]), y = exp(predict(m0)[o]), col = "magenta")
 142 
 143 
 144 ## @knitr LinearModelVisualization1, echo=F, results='markup', cache=T, fig.height=2, fig.width=3
 145 par(cex=.7, mar=c(4,4,1,0.1))
 146 plot(x = lexdec$Frequency, 
 147      y = lexdec$RT, 
 148      ylab = "(log-transformed) RT", 
 149      xlab = "(log-transformed) word frequency", 
 150      type = "n"
 151 )
 152 points(x = lexdec$Frequency, y = lexdec$RT, pch=1, cex=.1, col = "grey")
 153 abline(m, col = "green")
 154 
 155 
 156 ## @knitr LinearModelVisualization2, echo=F, results='markup', cache=T, fig.height=1.5, fig.width=3
 157 par(cex=.65, mar=c(4,4,0.5,0.1))
 158 plot(x = exp(lexdec$Frequency), 
 159      y = lexdec$RT, 
 160      ylab = "(log-transformed) RT", 
 161      xlab = "word frequency", 
 162      type = "n"
 163 )
 164 points(x = exp(lexdec$Frequency), y = lexdec$RT, pch=1, cex=.1, col = "grey")
 165 lines(x = exp(lexdec$Frequency[o]), y = predict(m)[o], col = "green")
 166 
 167 plot(x = exp(lexdec$Frequency), 
 168      y = exp(lexdec$RT), 
 169      ylab = "RT", 
 170      xlab = "word frequency", 
 171      type = "n"
 172 )
 173 points(x = exp(lexdec$Frequency), y = exp(lexdec$RT), pch=1, cex=.1, col = "grey")
 174 lines(x = exp(lexdec$Frequency[o]), y = exp(predict(m)[o]), col = "green")
 175 
 176 
 177 ## @knitr LinearModelExample3, echo=F, results='markup', cache=T, message=F
 178 summary(m)$coefficients
 179 
 180 
 181 ## @knitr LinearModelExample4, echo=T, results='markup', cache=T, message=F
 182 summary(m)
 183 
 184 
 185 ## @knitr LinearModelTwoPredictorsExample1, echo=T, results='markup', cache=T, message=F
 186 summary(glm(RT ~ Frequency + Trial, data = lexdec))$coefficients
 187 
 188 
 189 ## @knitr LinearModelTwoPredictorsVisualziation1, echo=F, results='markup', cache=T, message=F, fig.width=3.5, fig.height=3.5
 190 library(scatterplot3d)
 191 m2 = lm(RT ~ Frequency + Trial, data = lexdec)
 192 par(cex = .8, cex.lab = .6)
 193 s = scatterplot3d(x=lexdec$Frequency, y=lexdec$Trial, 
 194   z=lexdec$RT,
 195 	xlab="Word Frequency\n(log-transformed)", 
 196 	ylab="Trial number",
 197 	zlab="Response latency\n(in log-transformed msecs)", 
 198 	main="Predicting Lexical Decision RTs"
 199 )
 200 s$plane3d(m2, col= "purple", lwd=2, lty.box = "solid")
 201 
 202 
 203 ## @knitr LinearModelThreePredictorsExample1, echo=T, results='markup', cache=T, message=F
 204 s = summary(glm(RT ~ Frequency + Trial + NativeLanguage, data = lexdec))
 205 s$coefficients
 206 
 207 
 208 ## @knitr LinearModelThreePredictorsVisualziation1, echo=F, results='markup', cache=T, message=F, fig.width=3.5, fig.height=3
 209 m3 = lm(RT ~ Frequency + Trial + NativeLanguage, data = lexdec)
 210 
 211 par(cex = .8, cex.lab = .6)
 212 s = scatterplot3d(x=lexdec$Frequency, y=lexdec$Trial, 
 213   z=lexdec$RT,
 214   xlab="Word Frequency\n(log-transformed)", 
 215 	ylab="Trial number",
 216 	zlab="Response latency\n(in log-transformed msecs)", 
 217 	main="Predicting Lexical Decision RTs",
 218   sub="Native Speakers (red) and\nNon-Native Speakers (purple)"
 219 )
 220 s$plane3d(coef(m3)[1], coef(m3)[2], coef(m3)[3], 
 221 	col= "red", lwd=2, 
 222 	lty.box = "solid")
 223 s$plane3d(coef(m3)[1] + coef(m3)[4], coef(m3)[2], coef(m3)[3], 
 224 	col= "purple", lwd=2, 
 225 	lty.box = "solid")
 226 
 227 
 228 ## @knitr LinearModelInteractionExample1, echo=T, results='markup', cache=T, message=F
 229 m4 = glm(RT ~ 
 230            Frequency + 
 231            Trial + 
 232            NativeLanguage + 
 233            Frequency:NativeLanguage, 
 234          data = lexdec
 235 )
 236 summary(m4)$coefficients
 237 
 238 
 239 ## @knitr LinearModelInteractionVisualziation1, echo=F, results='markup', cache=T, message=F, fig.width=3.5, fig.height=3
 240 m4 = lm(RT ~ Frequency + Trial + NativeLanguage + Frequency:NativeLanguage, data = lexdec)
 241 
 242 par(cex = .8, cex.lab = .6)
 243 s = scatterplot3d(x=lexdec$Frequency, y=lexdec$Trial, 
 244   z=lexdec$RT,
 245   xlab="Word Frequency\n(log-transformed)", 
 246   ylab="Trial number",
 247 	zlab="Response latency\n(in log-transformed msecs)", 
 248 	main="Predicting Lexical Decision RTs",
 249   sub="Interaction with Native Speakers (red) and\nNon-Native Speakers (purple)"
 250 )
 251 s$plane3d(coef(m4)[1], coef(m4)[2], coef(m4)[3], 
 252 	col= "red", lwd=2, 
 253 	lty.box = "solid")
 254 s$plane3d(coef(m4)[1] + coef(m4)[4], coef(m4)[2] + coef(m4)[5], coef(m4)[3], 
 255 	col= "purple", lwd=2, 
 256 	lty.box = "solid")
 257 
 258 
 259 ## @knitr LinearModelInteractionExample2, echo=F, results='markup', cache=T, message=F
 260 summary(m4)$coefficients
Attachment 'LSA13-Lecture1-GLM.R'

Attached Files
