Linguist275_Final(Due3.23.15)
BPeloquin
March 16, 2015
I. Language Detection
library(stringr)
#########################################################
#initial data ------------------------------------------>
str.data = read.delim("~/Desktop/Winter_2014/Linguist275/HW/Final/final.txt",
header = F, sep = "\n", quote = "\"",
dec = ".", fill = TRUE, comment.char = "")
#initial string
str="#aababababaabaababaababaababa#"
#initial data ------------------------------------------>
#########################################################
#########################################################
#helper functions -------------------------------------->
#raw character counts for unconditional probabilites
raw.count = function(str) {
str=unlist(sapply(str, FUN=function(x){strsplit(x, split="")}))
a.count = 0
b.count = 0
hash.count = 0
#skipping the initial #???
for (n in 1:length(str)) {
if(str[n] == 'a') { a.count = a.count + 1 }
else if (str[n] == 'b') { b.count = b.count + 1 }
else { hash.count = hash.count + 1 }
}
return(c(a=a.count, b=b.count, hash=hash.count))
}
raw.count(str)## a b hash
## 17 11 2mless.mle = function(str) {
counts = raw.count(str)
return(counts / sum(counts))
}
#populate a 1st order markov chain transition matrix
#with 3 states
populate.transition = function(str) {
str=unlist(sapply(str, FUN=function(x){strsplit(x, split="")}))
#matrix to hold conditional counts
transition = matrix(
c(0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow=3, ncol=3)
colnames(transition) = c("a", "b", "#")
rownames(transition) = c("a", "b", "#")
#populate matrix with counts
n = 1
next.state = current.state = ""
while (next.state != "#") {
current.state = str[n]
next.state = str[n + 1]
transition[current.state, next.state] =
transition[current.state, next.state] + 1
n = n + 1
}
return(transition)
}
#Calculates relative frequencies across rows
matrix.probs = function(m) {
#Some strings may contain 0 of "a" or "b"
if(sum(m[1,]) != 0) { m[1,] = m[1,] / sum(m[1,]) }
if(sum(m[2,]) != 0) { m[2,] = m[2,] / sum(m[2,]) }
if(sum(m[3,]) != 0) { m[3,] = m[3,] / sum(m[3,]) }
return(m)
}
#Max likelihood estimate markov model
markov.mle = function(str) {
#raw counts
counts = raw.count(str)
#relative probabilites
rel.probs = counts / sum(counts)
#relative probabilities (unconditional) to vars
p.a = rel.probs[1]
p.b = rel.probs[2]
p.hash = rel.probs[3]
m = matrix.probs(populate.transition(str))
#conditional probabilities
p.a.m =
p.a * m["a","a"] + p.b * m["b","a"] + p.hash * m["#","a"]
p.b.m =
p.a * m["a","b"] + p.b * m["b","b"] + p.hash * m["#","b"]
p.hash.m =
p.a * m["a","#"] + p.b * m["b","#"] + p.hash * m["#", "#"]
return(c(a = p.a.m, b = p.b.m, hash = p.hash.m))
}
#Global for later use
mkov.transition = matrix.probs(populate.transition(str))
#helper functions -------------------------------------->
#########################################################a)
MLE estimates:
#########################################################
#model predictions-------------------------------------->
#MLE memoryless
mless.mle = mless.mle(str)
mless.mle## a b hash
## 0.56666667 0.36666667 0.06666667#sanity check
sum(mless.mle)## [1] 1#1st order Markov chain
mkov.mle = markov.mle(str)
mkov.mle## a.a b.a hash.a
## 0.60000000 0.36666667 0.03333333#sanity check
sum(mkov.mle)## [1] 1#model predictions-------------------------------------->
#########################################################MLE estimates using rejection sampling:
#Calculating param estimates using sampling techniques for
#m1 - memoryless and m2 - markov
#m1 (memoryless params estimation using sampling method)
rj.estimates = function(str, n.sims=1000) {
#actual frequencies observed
counts = raw.count(str)
size = nchar(str)
a.cnt = counts[1]
b.cnt = counts[2]
hash.cnt = counts[3]
samples = c()
while (length(samples) < (n.sims * 3)) {
sim.params = runif(3, 0, 1)
sim.params = sim.params / sum(sim.params)
a.samp.cnt = rbinom(n=1, size=size, sim.params[1])
b.samp.cnt = rbinom(n=1, size=size, sim.params[2])
hash.samp.cnt = size - (a.samp.cnt + b.samp.cnt)
if (a.samp.cnt == a.cnt && b.samp.cnt == b.cnt &&
hash.samp.cnt == hash.cnt) {
samples = rbind(samples, sim.params)
}
}
return(c(p.a = mean(samples[,1]),
p.b = mean(samples[,2]), p.hash = mean(samples[,3])))
}
#m2 (markov params estimation using sampling method)
#I'm essentially generating lots of possible strings given the
#parameter estimates from the original string and then looking at the MLE
#for the params
mkov.sampling2 = function(str, i=1000) {
mkov.transition = matrix.probs(populate.transition(str))
results = c()
p.hash = 1
for(n in seq(from=1, to=i)) {
#if current state is p.hash
if(p.hash) {
flip=runif(1,0,1)
if (flip <= mkov.transition[3, 1]) {
p.a = 1
p.b = 0
p.hash = 0
} else if (flip <= 1 - mkov.transition[3,3]) {
p.a = 0
p.b = 1
p.hash = 0
} else {
p.a = 0
p.b = 0
p.hash = 1
}
}
#if current state is p.a
else if (p.a) {
flip = runif(1, 0, 1)
if(flip <= mkov.transition[1,1]) {
p.a = 1
p.b = 0
p.hash = 0
}
else if (flip <= 1 - mkov.transition[1,3]) {
p.a = 0
p.b = 1
p.hash = 0
}
else {
p.a = 0
p.b = 0
p.hash = 1
}
}
#else current state is p.b
else {
flip = runif(1,0,1)
if (flip <= mkov.transition[2,1]) {
p.a = 1
p.b = 0
p.hash = 0
} else if (flip <= 1 - mkov.transition[2,3]) {
p.a = 0
p.b = 1
p.hash = 0
} else {
p.a = 0
p.b = 0
p.hash = 1
}
}
results = cbind(results, c(p.a, p.b, p.hash))
}
return(c(p.a=mean(results[1,]),
p.b = mean(results[2,]), p.hash = mean(results[3,])))
}
#Param estimates using sampling:
mmless.rject=rj.estimates(str)
mmless.rject## p.a p.b p.hash
## 0.5051154 0.3476126 0.1472719mkov.samp.est = mkov.sampling2(str)
mkov.samp.est## p.a p.b p.hash
## 0.580 0.382 0.038b) Entropy Rates
#Entropy for indiiv
entropy.rate = function(params) {
(params[1] * log2(1 / params[1]) +
params[2] * log2(1 / params[2]) +
params[3] * log2(1 / params[3])) / 3
}
#for some reason need to place here again
mless.mle = function(str) {
counts = raw.count(str)
return(counts / sum(counts))
}
#Entropy rates for our four param estimation methods
mmless.mle = mless.mle(str)
mkov.mle = markov.mle(str)
entropy.rate(mmless.mle)## a
## 0.4185122entropy.rate(mkov.mle)## a.a
## 0.3788258entropy.rate(mmless.rject)## p.a
## 0.4781977entropy.rate(mkov.samp.est)## p.a
## 0.3884793c) Classification of strings
#Two methods for estimating posteriors
#1 using dbinom
param.likelihoods1 = function(str, params) {
size = nchar(str)
counts = raw.count(str)
return(dbinom(counts[1], size, params[1]) *
dbinom(counts[2], size, params[2]) *
dbinom(counts[3], size, params[3]))
}
#2 using sampling
param.likelihoods2 = function(str, params, n.sims=1000) {
size = nchar(str)
counts = raw.count(str)
p.a = params[1]
p.b = params[2]
p.hash = params[3]
turns = 0
success = 0
while(success < n.sims) {
a.samp.cnt = rbinom(n=1, size=size, p.a)
b.samp.cnt = rbinom(n=1, size=size, p.b)
hash.samp.cnt = rbinom(n=1, size=size, p.hash)
if (a.samp.cnt == counts[1] &&
b.samp.cnt == counts[2] &&
hash.samp.cnt == counts[3]) {
success = success + 1
}
turns = turns + 1
}
return (success / turns)
}
bayes.factor = function(str, m1, m2) {
return(param.likelihoods1(str, m1) / param.likelihoods1(str, m2))
}
#mless.mle
#mkov.mle
#rj.estimates
#mkov.sampling2#predictions
###:::----------CAUTION SLOW RUN-TIME-------------:::###
results = c()
for (i in 1:nrow(str.data)) {
str = as.character(str.data[i,1])
length = nchar(str)
mless.est = rj.estimates(str)
markov.est = mkov.sampling2(str)
b.factor = bayes.factor(str, mless.est, markov.est)
m.preference = ifelse(b.factor > 1, "m1", "m2")
results = rbind(results, c(str, length, b.factor, m.preference))
}
colnames(results) = c("char.string", "string.len",
"b.factor", "model.preference")
###:::----------CAUTION SLOW RUN-TIME-------------:::###
results## char.string
## [1,] "#abbababababbabaabaabbaababbbaa#"
## [2,] "#abaababaaaaaabbbaabaaaa#"
## [3,] "#aaaaba#"
## [4,] "#babababab#"
## [5,] "#abaabababababbaba#"
## [6,] "#abaaabaababaaaaaaaaaaababbaaabaabbaabababaabababbbabaabaaababaaaaaaaaab#"
## [7,] "#aaaaba#"
## [8,] "#bababaababaabababaabababbabababbaababababbabaaabaabababababababaababababab#"
## [9,] "#abababbababa#"
## [10,] "#aababaaababa#"
## [11,] "#aaaa#"
## [12,] "#b#"
## [13,] "#bbab#"
## [14,] "#abbababa#"
## [15,] "#aaabaaaabaaba#"
## [16,] "#aabbabaabaabaababbababbababa#"
## [17,] "#aabba#"
## [18,] "#abbabaabababbbbababb#"
## [19,] "#aaaababaabaababaaaaababaaabaababababbabbababab#"
## [20,] "#aaaaaaabbaaba#"
## [21,] "#babaabaaabbbababbabbabababbabaaabbaaabbabaabababbababbba#"
## [22,] "#babababababa#"
## [23,] "#aaaa#"
## [24,] "#abaabaaababbbabaaaababbababbaabbabbabababbabaaaabababbbaa#"
## [25,] "#aababaaabaaaaababaaaaaabbaabaaaabbbababaabbbbbaaaaaabbabaaaabaabbabaaaaaaabbaaaaaaabaabaaaba#"
## [26,] "#abbbbaaaabbbaabaaaaabaabbabaaaabababaabaaaaaaaaabababaaaaaa#"
## [27,] "#ababaaaabbbaba#"
## [28,] "#aaaba#"
## [29,] "#bbbabaaaaaaabaaaaababa#"
## [30,] "#baabbababa#"
## [31,] "#aabbabaabbbaaaabbbbbaaab#"
## [32,] "#abaaaaaaababba#"
## [33,] "#abaababaababababa#"
## [34,] "#baaaaab#"
## [35,] "#ab#"
## [36,] "#baaaababababababaaabababababaababaaaabb#"
## [37,] "#aaaababaabaaaabbabbababaaaababbabbabababababb#"
## [38,] "#b#"
## [39,] "#aa#"
## [40,] "#ababaaba#"
## [41,] "#aaaaa#"
## [42,] "#abababababababbababbaaba#"
## [43,] "#aaababaabbaababababbababbabababababbabaaabababaaabaabaabababba#"
## [44,] "#abababaaabaabaabaaabaaaababababbababaabaabababbaababaabababaaababbaababbaabaa#"
## [45,] "#ababababaabbaabbabaabaabababaabababaababbbababaaaaaaababababbbabababbabab#"
## [46,] "#abababbabaababa#"
## [47,] "#a#"
## [48,] "#bbaaabbaaaabaababaaabaaaabbbbabbbaaabbaabaaabbaaaaaaaaaa#"
## [49,] "#baaaaaaabbbaabaaabaabaaababaabaab#"
## [50,] "#abaababababababababaaaababababababaaaaabaabababaaababaaabaababa#"
## string.len b.factor model.preference
## [1,] "32" "0.640487943009293" "m2"
## [2,] "25" "0.563350207916577" "m2"
## [3,] "8" "1.16389686193702" "m1"
## [4,] "11" "1.0314935830525" "m1"
## [5,] "19" "1.03971540114355" "m1"
## [6,] "73" "0.226854688987635" "m2"
## [7,] "8" "1.0083920526819" "m1"
## [8,] "76" "0.163311055045511" "m2"
## [9,] "14" "0.966647964765882" "m2"
## [10,] "14" "0.921021345520204" "m2"
## [11,] "6" "0.676372422060823" "m2"
## [12,] "3" "0.432281794612057" "m2"
## [13,] "6" "1.32549126663473" "m1"
## [14,] "10" "1.38164606434033" "m1"
## [15,] "15" "0.855773927119177" "m2"
## [16,] "30" "0.53376484433247" "m2"
## [17,] "7" "1.48035557227229" "m1"
## [18,] "22" "0.544700423621767" "m2"
## [19,] "48" "0.321960013512217" "m2"
## [20,] "15" "0.750402439933334" "m2"
## [21,] "58" "0.228680081256228" "m2"
## [22,] "14" "0.799924657942208" "m2"
## [23,] "6" "0.614026457387201" "m2"
## [24,] "59" "0.144864830732371" "m2"
## [25,] "94" "0.0821955685184869" "m2"
## [26,] "61" "0.274647007531366" "m2"
## [27,] "16" "0.872571429599374" "m2"
## [28,] "7" "1.22159632011216" "m1"
## [29,] "24" "0.461233510605669" "m2"
## [30,] "12" "1.24359148862928" "m1"
## [31,] "26" "0.488110158929704" "m2"
## [32,] "16" "0.772168553799172" "m2"
## [33,] "19" "0.814200982751354" "m2"
## [34,] "9" "0.995601203642725" "m2"
## [35,] "4" "1.08291964891957" "m1"
## [36,] "41" "0.5795914464702" "m2"
## [37,] "47" "0.29086365826729" "m2"
## [38,] "3" "0.415291225527665" "m2"
## [39,] "4" "0.612979279699205" "m2"
## [40,] "10" "1.0563843620882" "m1"
## [41,] "7" "0.592757774797543" "m2"
## [42,] "26" "0.560948597825319" "m2"
## [43,] "64" "0.183300070240158" "m2"
## [44,] "79" "0.116270263084536" "m2"
## [45,] "75" "0.210765851130452" "m2"
## [46,] "17" "0.910805285669443" "m2"
## [47,] "3" "0.416319338909822" "m2"
## [48,] "58" "0.415890797468166" "m2"
## [49,] "35" "0.345703707115387" "m2"
## [50,] "65" "0.168602497664647" "m2"##################################################
#plotting categories by length------------------->
plot(results[,"string.len"],
results[,"b.factor"], pch=16, col="blue",
main = "Categorizing strings",
ylab="Posterior odds", xlab="character length")
abline(h=1, col="red", lty=2)
#plotting categories by length------------------->
##################################################
m1.strings = results[which(results[,"model.preference"]=="m1"),"char.string"]
m2.strings = results[which(results[,"model.preference"]=="m2"),"char.string"]
aggregate.strings = function(data) {
string = ""
for(i in 1:length(data)) {
string = paste(string, data[i])
}
return(string)
}
m1.strings = aggregate.strings(m1.strings)
m2.strings = aggregate.strings(m2.strings)
#Edit out spaces
m1.strings=str_replace_all(m1.strings, "[^[:alnum:]#]", "")
m1.strings=str_replace_all(m1.strings, "##", "#")
m2.strings=str_replace_all(m2.strings, "[^[:alnum:]#]", "")
m2.strings=str_replace_all(m2.strings, "##", "#")
#Memoryless estimates
rj.estimates(m1.strings)## p.a p.b p.hash
## 0.4806862 0.3685673 0.1507465#Markov string estimates
mkov.sampling2(m2.strings)## p.a p.b p.hash
## 0.490 0.477 0.033matrix.probs(populate.transition(m2.strings))## a b #
## a 0.2666667 0.6666667 0.06666667
## b 0.6666667 0.3333333 0.00000000
## # 1.0000000 0.0000000 0.00000000d)
I’ve aggregated and combined the strings by the posterior and odds ratio. I’m combining the strings into one string and removing duplicate “##”. If these stochastic processes were actually generating these strings I assume the “#” would represent a logical break between two (so that the end “#” for one string would be the the beginning “#” for the next). With this new single string for each model I simply run my markov parameters estimates “mkov.sampling2()” and for the memoryless process “rj.estimates()”. I’ve displayed the transition matrix from my input. I’m a little confused about your distinction between the bayesian posterior odds and the ratio of likelihoods under the two models. Since the two models are equal a priori (uniform prior odds will cancel out) these two metrics should be identical…
II. Engaging the stats-for-psycholinguistics literature: Jaeger (2008)
a)
The fundamental issue with ANOVA is that it compares means of different experimental conditions and based on the sample variance within and between conditions, gives us information about whether to reject a null hypothesis. This is fine for continuous DV’s, however it is problematic for analysis of mean values for categorical variable in which we are likely looking at proportions of success or failures, etc. This is the crux of what Jaeger recommends against – using ANOVA over proportions. I thought Jaeger characterized one element of this problem nicely – ANOVAs for CDA have interpretability issues because they attribute probability mass to events that can never occur, “thereby underestimating the probability mass over events that can occur” (p.3). Furthermore, ANOVA makes assumptions about homogeneity of variance – a crucial problem when working with binomially distributed outcomes because the variance is only equal when p.1 and p.2 are equidistant from .5. This problem is exacerbated as the probabilities get more extreme (closer to 0 and 1). The arcsin-sqrt transformation has traditionally been considered as a method for dealing with some of these problems, however Jaeger shows that this approximation is also inadequate. Logistic regression, by contrast, is simply a linear model with a transformation (logit link function) of proportions into “space that is not bounded by 0 and 1” (log-odds space) (p.4). Additionally, this transformation accounts for the fact that with binomially distributed data, changes in proportions at the extremes should correspond with larger changes in the predictors than at .5. In this way, ordinary logit models better fit the nature of categorical outcomes. Jaeger continues with a host of other reasons why using ordinary logistic regressions is preferable to ANOVA. Among these are:
->Compatibility with multinomial distributions – in part this helps avoid the problems of data exclusion by being able to model the probability of errors.
->They are also more interpretable by providing information about the size and direction of an effect (coefficients in part)
->They are adaptable to imbalanced data
b)
The primary drawback with ordinary logit models is that they cannot account for random effects, which was likely a viable reason to go with ANOVA (despite its deficiencies) in the past. However, Jaeger introduces mixed logit models as a way around this. Without a logit link function, simply applying ANOVA to proportions is like running a linear probability model, which will have many of the same problems as ANOVA. The arcsine transformation is meant to approximate values as transformed by a logit link function (initially this was desirable because they were less costly, computationally), however this is no longer necessary.
c)
The ANOVA over the proportions of correct answers is inaccurate because there appears to be an effect of NP type on the proportion of correct answers for object RC’s compared to subject. However, this is in part because several of the subjects were at ceiling in their accuracy in Subject RC (Pronoun) conditions. ANOVA, even under the arcsine transformation, is unreliable at extremes around 0 or 1. The effect of NP on object RC accuracy appears larger because accuracy for subject RC was already so high. The ANOVA is unable to account for this as it simply considers the size of the difference. Of course, the difference is diminished when the outcome is naturally bounded (proportions) and subjects perform at ceiling. This is basically Jaeger’s first characterization of why ANOVAs are problematic with categorical data – they can assign probability mass to events which can’t exist.
d)
library(lme4)## Loading required package: Matrix
## Loading required package: Rcppd = data.frame(read.delim("~/Desktop/Winter_2014/Linguist275/HW/Final/inbal.tab",
header = T, sep = "\t", quote = "\"",
dec = ".", fill = TRUE, comment.char = ""))
d = subset(d, modality == 1) # comprehension study only
d = subset(d, Correct != "#NULL!" & !is.na(Extraction) & !is.na(NPType))
#convert to factors and relabel
names(d)## [1] "modality" "presentation" "instance" "sentence"
## [5] "itemby4" "itemby2" "child" "age"
## [9] "Extraction" "NPType" "animacymatch" "objectbias"
## [13] "headanimacy" "Response" "Correct" "CorrectRP"
## [17] "Reversal" "Agent" "OnlyObject" "OnlyObjectRP"
## [21] "Resumptive"d$NPType = factor(d$NPType)
levels(d$NPType) = c("lexical", "pronoun")
d$RCType = factor(d$Extraction)
levels(d$RCType) = c("object", "subject")
#for ordinary logit and mixed-effects models
d$Correct = factor(d$Correct)
levels(d$Correct) = c("incorrect", "correct")
d$Condition = paste(d$NPType, d$RCType)
#Verifying Percentage of correct answers
d$Correct = as.numeric(d$Correct)
d.Subject.Pronoun = d[which(d[,"NPType"] == "pronoun" &
d[,"RCType"] == "subject"),]
d.Object.Pronoun = d[which(d[,"NPType"] == "pronoun" &
d[,"RCType"] == "object"),]
d.Subject.Lexical = d[which(d[,"NPType"] == "lexical" &
d[,"RCType"] == "subject"),]
d.Object.Lexical = d[which(d[,"NPType"] == "lexical" &
d[,"RCType"] == "object"),]
length(which(d.Subject.Pronoun$Correct == 2)) /
length(d.Subject.Pronoun$Correct)## [1] 0.8433735length(which(d.Object.Pronoun$Correct == 2)) /
length(d.Object.Pronoun$Correct)## [1] 0.9565217length(which(d.Subject.Lexical$Correct == 2)) /
length(d.Subject.Lexical$Correct)## [1] 0.6890244length(which(d.Object.Lexical$Correct == 2)) /
length(d.Object.Lexical$Correct)## [1] 0.8956044#non-transformed aov() -------------->
#F1 analysis - Subject analysis
summary(aov(as.numeric(Correct) ~ NPType * RCType +
Error(as.factor(child)), data = d))##
## Error: as.factor(child)
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 0.559 0.5587 2.011 0.1716
## RCType 1 0.010 0.0103 0.037 0.8491
## NPType:RCType 1 1.430 1.4304 5.148 0.0345 *
## Residuals 20 5.558 0.2779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 1.97 1.967 17.755 2.86e-05 ***
## RCType 1 4.42 4.423 39.918 4.83e-10 ***
## NPType:RCType 1 0.39 0.395 3.561 0.0596 .
## Residuals 669 74.12 0.111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #same analysis using lmer
anova(lmer(Correct ~ 1 + RCType * NPType + (1 | child), data=d))## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## RCType 1 4.4190 4.4190 39.8877
## NPType 1 1.9540 1.9540 17.6374
## RCType:NPType 1 0.3889 0.3889 3.5108#F2 analysis - Item analysis
summary(aov(as.numeric(Correct) ~ NPType * RCType +
Error(as.factor(itemby4)), data = d))##
## Error: as.factor(itemby4)
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 0.0160 0.01597 0.053 0.829
## RCType 1 0.0248 0.02478 0.082 0.788
## NPType:RCType 1 0.0469 0.04695 0.156 0.713
## Residuals 4 1.2049 0.30122
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 1.91 1.914 16.340 5.89e-05 ***
## RCType 1 4.60 4.598 39.242 6.61e-10 ***
## NPType:RCType 1 0.40 0.402 3.428 0.0645 .
## Residuals 685 80.25 0.117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lmer(Correct ~ 1 + RCType * NPType + (1 | itemby4), data=d))## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## RCType 1 4.4909 4.4909 38.3404
## NPType 1 1.9253 1.9253 16.4376
## RCType:NPType 1 0.3887 0.3887 3.3182#F2 analysis - Combined analysis
summary(aov(as.numeric(Correct) ~ NPType * RCType +
Error(as.factor(child) / as.factor(itemby4)), data = d))##
## Error: as.factor(child)
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 0.559 0.5587 2.011 0.1716
## RCType 1 0.010 0.0103 0.037 0.8491
## NPType:RCType 1 1.430 1.4304 5.148 0.0345 *
## Residuals 20 5.558 0.2779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: as.factor(child):as.factor(itemby4)
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 0.000 0.00000 0.000 0.997
## RCType 1 0.022 0.02224 0.200 0.655
## NPType:RCType 1 0.276 0.27626 2.484 0.117
## Residuals 165 18.354 0.11124
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 2.03 2.027 18.323 2.23e-05 ***
## RCType 1 4.45 4.452 40.237 5.03e-10 ***
## NPType:RCType 1 0.34 0.345 3.118 0.078 .
## Residuals 501 55.43 0.111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#arc-sin transformed aov() - Subject analysis only (as referenced in Jaeger (2006))
d$a.sin.transform = asin(sqrt(as.numeric(d$Correct) / 100))
summary(aov(a.sin.transform ~ NPType * RCType + Error(child), data = d))##
## Error: child
## Df Sum Sq Mean Sq
## NPType 1 0.000539 0.000539
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## NPType 1 0.00336 0.003364 16.390 5.74e-05 ***
## RCType 1 0.00766 0.007659 37.320 1.67e-09 ***
## NPType:RCType 1 0.00066 0.000662 3.224 0.073 .
## Residuals 691 0.14182 0.000205
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The results of my ANOVA show that the interaction is significant, but not nearly to the degree observed in Jaeger (2006)… I’m finding significance at .1 suggesting a fairly high likelihood of Type 1 error. My F-value for RCType is inflated compared to Jaeger’s which would indicate that this anova is actually doing a better job assessing the size of the effect from this variable, perhaps pulling from the interaction.
#Ordinary logit with interaction
#recode as factor
d$Correct = factor(d$Correct)
levels(d$Correct) = c("incorrect", "correct")
m.ordinary.logit = glm(Correct ~ RCType * NPType,
family=binomial, data=d)
summary(m.ordinary.logit)##
## Call:
## glm(formula = Correct ~ RCType * NPType, family = binomial, data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5042 0.2982 0.4696 0.5837 0.8631
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.14931 0.24242 8.866 < 2e-16 ***
## RCTypesubject -1.35375 0.29534 -4.584 4.57e-06 ***
## NPTypepronoun 0.94173 0.43521 2.164 0.0305 *
## RCTypesubject:NPTypepronoun -0.05375 0.51329 -0.105 0.9166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 587.02 on 695 degrees of freedom
## Residual deviance: 535.04 on 692 degrees of freedom
## AIC: 543.04
##
## Number of Fisher Scoring iterations: 5#mixed logistic regression
library(lme4)
mr.1 = lmer(Correct ~ 1 + RCType * NPType +
(1 + RCType * NPType | child), data=d)
mr.2 = lmer(Correct ~ 1 + RCType * NPType +
(1 + as.numeric(RCType) | child), data=d)
mr.3 = lmer(Correct ~ 1 + RCType * NPType +
(1 + as.numeric(RCType) | child) + (1 | itemby4), data=d)
summary(mr.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Correct ~ 1 + RCType * NPType + (1 + RCType * NPType | child)
## Data: d
##
## REML criterion at convergence: 480.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.93175 0.05601 0.21768 0.46547 1.71506
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## child (Intercept) 0.007437 0.08624
## RCTypesubject 0.007736 0.08796 0.25
## NPTypepronoun 0.002063 0.04542 -1.00 -0.22
## RCTypesubject:NPTypepronoun 0.002151 0.04637 0.72 -0.48 -0.75
## Residual 0.108673 0.32966
## Number of obs: 696, groups: child, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.89525 0.03014 62.88
## RCTypesubject -0.20803 0.03979 -5.23
## NPTypepronoun 0.06114 0.03570 1.71
## RCTypesubject:NPTypepronoun 0.09516 0.05095 1.87
##
## Correlation of Fixed Effects:
## (Intr) RCTyps NPTypp
## RCTypesbjct -0.432
## NPTypepronn -0.707 0.396
## RCTypsb:NPT 0.468 -0.662 -0.690summary(mr.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Correct ~ 1 + RCType * NPType + (1 + as.numeric(RCType) | child)
## Data: d
##
## REML criterion at convergence: 482.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.90881 0.08094 0.25874 0.47210 1.68366
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## child (Intercept) 0.0001652 0.01285
## as.numeric(RCType) 0.0046405 0.06812 -1.00
## Residual 0.1096192 0.33109
## Number of obs: 696, groups: child, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.89535 0.02702 70.15
## RCTypesubject -0.20795 0.03827 -5.43
## NPTypepronoun 0.06098 0.03461 1.76
## RCTypesubject:NPTypepronoun 0.09512 0.05027 1.89
##
## Correlation of Fixed Effects:
## (Intr) RCTyps NPTypp
## RCTypesbjct -0.431
## NPTypepronn -0.644 0.455
## RCTypsb:NPT 0.443 -0.661 -0.689summary(mr.3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Correct ~ 1 + RCType * NPType + (1 + as.numeric(RCType) | child) +
## (1 | itemby4)
## Data: d
##
## REML criterion at convergence: 480.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95055 0.03911 0.25042 0.48734 1.72317
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## child (Intercept) 0.0001845 0.01358
## as.numeric(RCType) 0.0047023 0.06857 -1.00
## itemby4 (Intercept) 0.0011587 0.03404
## Residual 0.1085686 0.32950
## Number of obs: 696, groups: child, 24; itemby4, 8
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.89561 0.02946 64.34
## RCTypesubject -0.21099 0.03832 -5.51
## NPTypepronoun 0.06104 0.03445 1.77
## RCTypesubject:NPTypepronoun 0.09694 0.05033 1.93
##
## Correlation of Fixed Effects:
## (Intr) RCTyps NPTypp
## RCTypesbjct -0.390
## NPTypepronn -0.588 0.453
## RCTypsb:NPT 0.403 -0.661 -0.685Like Jaeger, my ordinary logistic and mixed logistic models are correctly identifying significant contributions from the RCType and NPType predictors, but not an interaction. As with my ANOVAs, however, my numbers are different than Jaeger’s. At this point I’m not sure what I’m doing wrong. I know Jaeger used a different regression library and functions (lrm from Design library) and I wonder if some of these differences come from the implementations of the functions.