1 library(languageR)
   2 
   3 ## ---------------------------------------------------
   4 # Lexical decision time data
   5 #
   6 # Q: To what extent are lexical decisions influenced by
   7 #    (a) properties of the participants, (b) words (items) 
   8 #    used in the experiment, (c) design properties of the 
   9 #    experiment, or (d) more general factors?
  10 ## ---------------------------------------------------
  11 data(lexdec)
  12 library(lme4, keep.source = FALSE)
  13 library(arm)
  14 
  15 lm.simp <- lm(RT ~ 1, data=lexdec)
  16 display(lm.simp)
  17 
  18 lm.1 <- lm(RT ~ meanWeight + Frequency, data=lexdec)
  19 summary(lm.1)
  20 vif(lm.1)
  21 
  22 lm.2 <- lm(RT ~ 1 + Length + meanWeight + Frequency + FamilySize + SynsetCount
  23 		   , data=lexdec)
  24 summary(lm.2)
  25 
  26 # comparison of models
  27 anova(lm.2, lm.1, lm.simp)
  28 
  29 lexdec$SubjFreqRatio <- lexdec$SubjFreq - lexdec$Frequency
  30 lm.full <- lm(RT ~ 1 + meanWeight + Frequency + Length + DerivEntropy + SubjFreqRatio
  31 			, data= lexdec)
  32 summary(lm.full)
  33 vif(lm.full)
  34 
  35 cor(lexdec[,c('meanWeight','Frequency','Length','DerivEntropy','SubjFreqRatio')])
  36 
  37 lm.full <- lm(RT ~ 1 + Frequency + SubjFreqRatio
  38 			, data= lexdec)
  39 summary(lm.full)
  40 
  41 ## ---------------------------------------------------
  42 # investigate the factors in the data set and build a good model 
  43 # ONLY USING LINGUISTIC FACTORS
  44 ## ---------------------------------------------------
  45 
  46 ## ---------------------------------------------------
  47 # Do subjects differ?
  48 ## ---------------------------------------------------
  49 
  50 ll.1 <- lmList(RT ~ 1 + Frequency + SubjFreqRatio | Subject, data=lexdec)
  51 ll.1
  52 str(ll.1)
  53 
  54 # a little function to calculate the standard error of a 
  55 # sample (in a vector)
  56 se <- function(x) {
  57 	sqrt(var(x)) / sqrt(length(x))
  58 }
  59 
  60 # coef and standard error of coef from each model
  61 paste("[", sapply(ll.1, coef), ", ", sapply(ll.1, se.coef), "]", sep="")
  62 
  63 # coef and standard error across model
  64 paste("[", mean(sapply(ll.1, coef)[2,]), ", ", se(sapply(ll.1, se.coef)[,2]), "]", sep="")
  65 
  66 
  67 # Trellis plot of individual participant coefficients, xyplot
  68 trellis.device(color=F)
  69 xyplot(RT ~ Frequency | Subject,
  70 	data=lexdec, 
  71 	main="By-subject LMs",
  72 	ylab="log reaction times",
  73 	xlab="",
  74 #	ylim= c(-200,200),
  75 	panel=function(x, y){
  76 		panel.xyplot(x, y, col=1)
  77 		panel.lmline(x, y, lty=4, col="blue", lwd=3)
  78 	}	
  79 )
  80 
  81 old.prompt = grid::grid.prompt(T)
  82 qqmath( ~ RT | Subject, data= lexdec, layout=c(4,4),
  83        prepanel = prepanel.qqmathline,
  84        panel = function(x, ...) {
  85           panel.qqmathline(x, col= 2, lty=1, ...)
  86           panel.qqmath(x, ...)
  87        }
  88 )
  89 grid::grid.prompt(old.prompt)
  90 
  91 xylowess.fnc(RT ~ Frequency | Subject, data=lexdec)
  92 
  93 
  94 ## ---------------------------------------------------
  95 # multilevel (aka mixed) model
  96 ## ---------------------------------------------------
  97 
  98 lmer.1 <- lmer(RT ~ 1 + Frequency + SubjFreqRatio +
  99 		    (1 | Subject)
 100 		   , data=lexdec)
 101 summary(lmer.1)
 102 
 103 # the random effects
 104 ranef(lmer.1)
 105 var(ranef(lmer.1)[[1]])
 106 
 107 
 108 ## ---------------------------------------------------
 109 # R2
 110 ## ---------------------------------------------------
 111 summary(lm.full)
 112 cor(fitted(lm.full), lexdec$RT)^ 2
 113 
 114 cor(fitted(lmer.1), lexdec$RT) ^ 2
 115 
 116 
 117 ## ---------------------------------------------------
 118 # use R2 comparisons to investigate
 119 # item effects
 120 # design effects
 121 # effects of subject differences
 122 ## ---------------------------------------------------
 123 
 124 lmer.item <- lmer(RT ~ 1 + Frequency + SubjFreqRatio +
 125 		    (1 | Word) + (1 | Subject)
 126 		   , data=lexdec)
 127 cor(fitted(lmer.item), lexdec$RT) ^ 2
 128 
 129 
 130 anova(lmer.1, lmer.item)
 131 
 132 
 133 
 134 
 135 lmer.2 <- lmer(RT ~ 1 + Frequency + SubjFreqRatio +
 136 		    (1 | Word) + (1 + Frequency | Subject)
 137 		   , data=lexdec)
 138 lmer.2
 139 
 140 
 141 
 142 
 143 
 144 
 145 
 146 
 147 ## ---------------------------------------------------
 148 # full models
 149 ## ---------------------------------------------------
 150 
 151 lmer.full = lmer(RT ~ 1 + Correct + Trial + PrevType * meanWeight + 
 152    Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), 
 153    data = lexdec)
 154 anova(lmer.full, lmer.1)
 155 
 156 p <- pvals.fnc(lmer.full, withMCMC=T)
 157 p$summary
 158 
 159 
 160 
 161 
 162 lmer.4 <- lmer(RT ~ 1 + Correct + PrevCorrect + Trial + PrevType + meanWeight + 
 163    Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), 
 164    data = lexdec)
 165 lmer.4
 166 anova(lmer.full, lmer.1)