Non-psychometric MLE
Olsson et al.
07/11/2020
Non-psychometric MLE chromatic background.
Load Package for Mixed-Effects Modelling
library('readr')
library('boot')
library('lme4')## Loading required package: Matrixlibrary('MASS')
library('emmeans')Read in data.
col.dta <- read_delim('colour_discriminate_short_format.txt', delim='\t')## Parsed with column specification:
## cols(
## Colour.difference = col_double(),
## pcorr = col_double(),
## corr = col_double(),
## incorr = col_double(),
## ind = col_character(),
## background = col_character(),
## sex = col_character(),
## stimuli = col_character(),
## batch = col_character()
## )data <- col.dta
data$stimuli <- factor(data$stimuli)
data$target <- with(data, ifelse(background == stimuli, 'same', 'diff'))
data$target <- as.factor(data$target)
data$background <- factor(data$background)
data$ind <- factor(data$ind)
data$batch <- as.factor(data$batch)
data$chick <- paste0(data$ind, data$batch)Useful Functions
resplot <- function(mod){
#are residuals normally distributed
hist(residuals(mod), prob = T, xlab = formula(mod), main = paste('Residuals for',(formula(mod)[3])))
lines(density(residuals(mod)), col = 'red')
lines(seq(min(residuals(mod)),max(residuals(mod)), length.out = 10^3), dnorm(seq(min(residuals(mod)),max(residuals(mod)), length.out = 10^3), 0, sd(residuals(mod))), col = 'blue')
legend('topright', legend = c('kernel density', 'fitted normal'), lty = 1, col = c('red', 'blue'))
boxplot(residuals(mod),
add = T, axes = F, horizontal = T, cex = 0.5, outline = T, border = rgb(0,0.1,0,0.7), at = par('yaxp')[2]*0.1,
pars = list(boxwex = par('yaxp')[2]*0.3, staplewex = par('yaxp')[2]*0.5, outwex = par('yaxp')[2]*0.05))#
} Fit Logistic Regression
Logistic regression fits a relationship between a binomial variable (yes|no, correct|incorrect, present|absent, 1|0) and any other predictor variable. Variance in binomial variable is constrained by the fact that there cannot be > 100% 1s or 0s in a dataset, so relationships always follows an “S” curve. Fit a mixed-effects (both controlled [‘fixed’] and unpredictable [‘random’] effects) logistic regression predicting ‘response’ (a correct or incorrect) and taking into account different intercepts (starting biases: “…1|ind”) and slopes (learning rates: “(Colour.difference-1|…”) of different inds, using the Generlized Linear Mixed Effects Regression (glmer) command
N.B. There are not enough observations to estimate random effects of chick and batch across Colour.difference:target:background
As it happens, batches are shared across backgrounds whereas chicks are unique to each condition. It might make more sense to look at random effects of chick on colour difference alone.
nonpsych1 <- glmer(cbind(corr,incorr)~Colour.difference*background*target*sex +
(1+Colour.difference|chick)+
(1+Colour.difference*background|batch),
data = data, family = binomial(link = 'logit'))Needs a little help converging, increase tolerance.
print(.Machine$double.eps * 10^8)## [1] 2.220446e-08nonpsych2 <- glmer(cbind(corr,incorr)~Colour.difference*background*target*sex +
(1+Colour.difference|chick)+
(1+Colour.difference*background|batch),
data = data, family = binomial(link = 'logit'),
control = glmerControl(tol = .Machine$double.eps * 10^8))Fit a logistic regression using no information about experience, but controlling for the same random effects, as a test for an effect of experience on response (learning)
nonpsych0 <- glmer(cbind(corr,incorr)~1 + (1|chick)+ (1|batch),
data = data, family = binomial(link = 'logit'))Compare the variance in the data explained by experience
anova(nonpsych0, nonpsych1, nonpsych2, test = 'Chisq') ## Data: data
## Models:
## nonpsych0: cbind(corr, incorr) ~ 1 + (1 | chick) + (1 | batch)
## nonpsych1: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych1: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych1: background | batch)
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych0 3 1941.5 1950.8 -967.74 1935.47
## nonpsych1 29 1001.9 1091.8 -471.93 943.86 991.61 26 <2e-16 ***
## nonpsych2 29 1001.9 1091.8 -471.93 943.86 0.00 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model taking experience into account deviates less from the data recorded (there is an effect of experience).
anova(nonpsych0, nonpsych2, test = 'Chisq')## Data: data
## Models:
## nonpsych0: cbind(corr, incorr) ~ 1 + (1 | chick) + (1 | batch)
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych0 3 1941.5 1950.8 -967.74 1935.47
## nonpsych2 29 1001.9 1091.8 -471.93 943.86 991.61 26 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For more complex models, Akaike Information Criteria (AIC) may be compared.
AIC(nonpsych0, nonpsych1, nonpsych2) ## df AIC
## nonpsych0 3 1941.473
## nonpsych1 29 1001.862
## nonpsych2 29 1001.863The model with fixed effects has lower AIC, indicating a better fit
resplot(nonpsych2)#looks good
shapiro.test(residuals(nonpsych2)) # even passes a Shapiro-Wilk test##
## Shapiro-Wilk normality test
##
## data: residuals(nonpsych2)
## W = 0.98809, p-value = 0.1793no unfitted variables
sum(is.na(coef(summary(nonpsych2))[,'Estimate']))## [1] 0sum(!is.na(coef(summary(nonpsych2))[,'Estimate']))## [1] 16Get the model summary for nonpsych2.
summary(nonpsych2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(corr, incorr) ~ Colour.difference * background * target *
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch)
## Data: data
## Control: glmerControl(tol = .Machine$double.eps * 10^8)
##
## AIC BIC logLik deviance df.resid
## 1001.9 1091.8 -471.9 943.9 135
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.3666 -0.9807 -0.0180 0.8258 2.7466
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## chick (Intercept) 4.132e-01 0.6428160
## Colour.difference 2.074e-01 0.4553683 -0.94
## batch (Intercept) 3.575e-08 0.0001891
## Colour.difference 4.084e-07 0.0006391 -0.41
## backgroundorange 4.207e-02 0.2051139 0.49 0.59
## Colour.difference:backgroundorange 6.075e-03 0.0779422 0.50 0.58
##
##
##
##
##
##
## 1.00
## Number of obs: 164, groups: chick, 32; batch, 4
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) -0.249389 0.332256
## Colour.difference 0.896342 0.226750
## backgroundorange -0.003663 0.522932
## targetsame 0.226848 0.469968
## sexmale -0.339054 0.542146
## Colour.difference:backgroundorange 0.146706 0.351972
## Colour.difference:targetsame 0.311576 0.332136
## backgroundorange:targetsame -0.762745 0.714792
## Colour.difference:sexmale 0.025165 0.368431
## backgroundorange:sexmale 0.177818 0.771854
## targetsame:sexmale 0.361086 0.760550
## Colour.difference:backgroundorange:targetsame 0.197734 0.507653
## Colour.difference:backgroundorange:sexmale -0.151522 0.534113
## Colour.difference:targetsame:sexmale -0.086605 0.535967
## backgroundorange:targetsame:sexmale 0.295686 1.090496
## Colour.difference:backgroundorange:targetsame:sexmale 0.031113 0.777085
## z value Pr(>|z|)
## (Intercept) -0.751 0.453
## Colour.difference 3.953 7.72e-05 ***
## backgroundorange -0.007 0.994
## targetsame 0.483 0.629
## sexmale -0.625 0.532
## Colour.difference:backgroundorange 0.417 0.677
## Colour.difference:targetsame 0.938 0.348
## backgroundorange:targetsame -1.067 0.286
## Colour.difference:sexmale 0.068 0.946
## backgroundorange:sexmale 0.230 0.818
## targetsame:sexmale 0.475 0.635
## Colour.difference:backgroundorange:targetsame 0.390 0.697
## Colour.difference:backgroundorange:sexmale -0.284 0.777
## Colour.difference:targetsame:sexmale -0.162 0.872
## backgroundorange:targetsame:sexmale 0.271 0.786
## Colour.difference:backgroundorange:targetsame:sexmale 0.040 0.968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## convergence code: 0
## Model failed to converge with max|grad| = 0.0541608 (tol = 0.002, component 1)Look for candidate variables to remove
anova(nonpsych2)[order(anova(nonpsych2)$`F value`),]## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Colour.difference:background:target:sex 1 0.002 0.002 0.0017
## Colour.difference:target:sex 1 0.041 0.041 0.0412
## background 1 0.113 0.113 0.1133
## Colour.difference:background:sex 1 0.169 0.169 0.1691
## Colour.difference:background:target 1 0.280 0.280 0.2804
## Colour.difference:sex 1 0.323 0.323 0.3234
## background:sex 1 0.542 0.542 0.5418
## background:target:sex 1 0.606 0.606 0.6065
## Colour.difference:background 1 0.783 0.783 0.7828
## sex 1 0.913 0.913 0.9128
## background:target 1 2.352 2.352 2.3516
## target:sex 1 3.669 3.669 3.6688
## Colour.difference:target 1 4.021 4.021 4.0213
## target 1 29.238 29.238 29.2379
## Colour.difference 1 129.104 129.104 129.1042Remove top level interaction
nonpsych.a<- update(nonpsych2,.~.-Colour.difference:background:target:sex)#4way
formula(nonpsych.a)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + Colour.difference:sex + background:sex +
## target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## Colour.difference:target:sex + background:target:sexRemove Colour.difference:target:sex
nonpsych.b<- update(nonpsych.a,.~.-Colour.difference:target:sex)#3way
formula(nonpsych.b)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + Colour.difference:sex + background:sex +
## target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## background:target:sexRemove Colour.difference:background:sex
nonpsych.c<- update(nonpsych.b,.~.-Colour.difference:background:sex )#3way
formula(nonpsych.c)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + Colour.difference:sex + background:sex +
## target:sex + Colour.difference:background:target + background:target:sexRemove Colour.difference:background:target
nonpsych.d<- update(nonpsych.c,.~.-Colour.difference:background:target )#3way
formula(nonpsych.d)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + Colour.difference:sex + background:sex +
## target:sex + background:target:sexRemove background:target:sex
nonpsych.e<- update(nonpsych.d,.~.-background:target:sex )#3way
formula(nonpsych.e)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + Colour.difference:sex + background:sex +
## target:sexRemove Colour.difference:sex
nonpsych.f<- update(nonpsych.e,.~.-Colour.difference:sex )#2way
formula(nonpsych.f)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + background:sex + target:sexRemove background:sex
nonpsych.g<- update(nonpsych.f,.~.-background:sex )#2way
formula(nonpsych.g)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## background | batch) + Colour.difference:background + Colour.difference:target +
## background:target + target:sexIf we remove Colour.difference:background, we also have to remove (1+Colour.difference*background|batch) and replace it with + (1+Colour.difference+background|batch)
Otherwise all variance attributed to the interaction can be still be attributed to its random effect. I find that implausible.
nonpsych.h<- update(nonpsych.g,.~.-Colour.difference:background -
(1+Colour.difference*background|batch) +
+ (1+Colour.difference+background|batch) )#2way
formula(nonpsych.h)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## background | batch) + Colour.difference:target + background:target +
## target:sexnonpsych.i<- update(nonpsych.h,.~.-background:target)#2way
formula(nonpsych.i)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## background | batch) + Colour.difference:target + target:sexnonpsych.j<- update(nonpsych.i, ~.-Colour.difference:target )#2way
formula(nonpsych.j)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## background | batch) + target:sexnonpsych.k<- update(nonpsych.j,.~.-target:sex)#2way
formula(nonpsych.k)## cbind(corr, incorr) ~ Colour.difference + background + target +
## sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## background | batch)If we remove background, we also have to remove (1+Colour.difference+background|batch) and replace it with + (1+Colour.difference|batch)
nonpsych.l<- update(nonpsych.k,.~.-background -
(1+Colour.difference+background|batch) +
(1+Colour.difference|batch) )#1way
formula(nonpsych.l)## cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## Colour.difference | chick) + (1 + Colour.difference | batch)nonpsych.m<- update(nonpsych.l,.~.-sex)#1way
formula(nonpsych.m)## cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## chick) + (1 + Colour.difference | batch)nonpsych.n<- update(nonpsych.m,.~.-target)#1way
formula(nonpsych.n)## cbind(corr, incorr) ~ Colour.difference + (1 + Colour.difference |
## chick) + (1 + Colour.difference | batch)If we remove Colour.difference, we also have to remove (1+Colour.difference|batch) and (1+Colour.difference|chick) and replace them with (1|chick) + (1|batch)
nonpsych.o<- update(nonpsych.n,.~.-Colour.difference -
(1+Colour.difference|chick) +
(1+Colour.difference|batch) -
(1|chick) + (1|batch) )formula(nonpsych.o)#same as nonpsych0, no fixed effects left to remove## cbind(corr, incorr) ~ (1 + Colour.difference | batch) + (1 |
## batch)anova(nonpsych2, nonpsych.a, nonpsych.b, nonpsych.c, nonpsych.d, nonpsych.e, nonpsych.f, nonpsych.g, nonpsych.h, nonpsych.i, nonpsych.j, nonpsych.k, nonpsych.l, nonpsych.m, nonpsych.n, nonpsych.o)## Data: data
## Models:
## nonpsych.o: cbind(corr, incorr) ~ (1 + Colour.difference | batch) + (1 |
## nonpsych.o: batch)
## nonpsych.n: cbind(corr, incorr) ~ Colour.difference + (1 + Colour.difference |
## nonpsych.n: chick) + (1 + Colour.difference | batch)
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## nonpsych.l: cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## nonpsych.l: Colour.difference | chick) + (1 + Colour.difference | batch)
## nonpsych.k: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.k: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.k: background | batch)
## nonpsych.j: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.j: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.j: background | batch) + target:sex
## nonpsych.i: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.i: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.i: background | batch) + Colour.difference:target + target:sex
## nonpsych.h: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.h: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.h: background | batch) + Colour.difference:target + background:target +
## nonpsych.h: target:sex
## nonpsych.g: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.g: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.g: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.g: background:target + target:sex
## nonpsych.f: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.f: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.f: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.f: background:target + background:sex + target:sex
## nonpsych.e: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.e: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.e: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.e: background:target + Colour.difference:sex + background:sex +
## nonpsych.e: target:sex
## nonpsych.d: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.d: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.d: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.d: background:target + Colour.difference:sex + background:sex +
## nonpsych.d: target:sex + background:target:sex
## nonpsych.c: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.c: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.c: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.c: background:target + Colour.difference:sex + background:sex +
## nonpsych.c: target:sex + Colour.difference:background:target + background:target:sex
## nonpsych.b: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.b: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.b: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.b: background:target + Colour.difference:sex + background:sex +
## nonpsych.b: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.b: background:target:sex
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.o 5 1127.27 1142.8 -558.63 1117.27
## nonpsych.n 8 989.11 1013.9 -486.56 973.11 144.1549 3 < 2.2e-16 ***
## nonpsych.m 9 976.35 1004.2 -479.17 958.35 14.7681 1 0.0001216 ***
## nonpsych.l 10 977.99 1009.0 -478.99 957.99 0.3568 1 0.5503130
## nonpsych.k 14 983.33 1026.7 -477.66 955.33 2.6619 4 0.6158920
## nonpsych.j 15 982.87 1029.4 -476.44 952.87 2.4527 1 0.1173220
## nonpsych.i 16 981.39 1031.0 -474.70 949.39 3.4825 1 0.0620223 .
## nonpsych.h 17 980.22 1032.9 -473.11 946.22 3.1688 1 0.0750589 .
## nonpsych.g 22 989.55 1057.7 -472.77 945.55 0.6770 5 0.9842072
## nonpsych.f 23 991.07 1062.4 -472.53 945.07 0.4782 1 0.4892392
## nonpsych.e 24 992.92 1067.3 -472.46 944.92 0.1463 1 0.7020980
## nonpsych.d 25 994.31 1071.8 -472.15 944.31 0.6126 1 0.4338281
## nonpsych.c 26 996.01 1076.6 -472.01 944.01 0.2978 1 0.5852428
## nonpsych.b 27 997.89 1081.6 -471.95 943.89 0.1188 1 0.7303385
## nonpsych.a 28 999.86 1086.7 -471.93 943.86 0.0285 1 0.8658964
## nonpsych2 29 1001.86 1091.8 -471.93 943.86 0.0010 1 0.9746008
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(nonpsych2, nonpsych.a, nonpsych.b, nonpsych.c, nonpsych.d, nonpsych.e, nonpsych.f, nonpsych.g, nonpsych.h, nonpsych.i, nonpsych.j, nonpsych.k, nonpsych.l, nonpsych.m, nonpsych.n, nonpsych.o)[order(row.names(anova(nonpsych2, nonpsych.a, nonpsych.b, nonpsych.c, nonpsych.d, nonpsych.e, nonpsych.f, nonpsych.g, nonpsych.h, nonpsych.i, nonpsych.j, nonpsych.k, nonpsych.l, nonpsych.m, nonpsych.n, nonpsych.o))),]## Data: data
## Models:
## nonpsych.o: cbind(corr, incorr) ~ (1 + Colour.difference | batch) + (1 |
## nonpsych.o: batch)
## nonpsych.n: cbind(corr, incorr) ~ Colour.difference + (1 + Colour.difference |
## nonpsych.n: chick) + (1 + Colour.difference | batch)
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## nonpsych.l: cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## nonpsych.l: Colour.difference | chick) + (1 + Colour.difference | batch)
## nonpsych.k: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.k: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.k: background | batch)
## nonpsych.j: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.j: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.j: background | batch) + target:sex
## nonpsych.i: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.i: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.i: background | batch) + Colour.difference:target + target:sex
## nonpsych.h: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.h: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.h: background | batch) + Colour.difference:target + background:target +
## nonpsych.h: target:sex
## nonpsych.g: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.g: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.g: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.g: background:target + target:sex
## nonpsych.f: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.f: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.f: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.f: background:target + background:sex + target:sex
## nonpsych.e: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.e: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.e: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.e: background:target + Colour.difference:sex + background:sex +
## nonpsych.e: target:sex
## nonpsych.d: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.d: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.d: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.d: background:target + Colour.difference:sex + background:sex +
## nonpsych.d: target:sex + background:target:sex
## nonpsych.c: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.c: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.c: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.c: background:target + Colour.difference:sex + background:sex +
## nonpsych.c: target:sex + Colour.difference:background:target + background:target:sex
## nonpsych.b: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.b: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.b: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.b: background:target + Colour.difference:sex + background:sex +
## nonpsych.b: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.b: background:target:sex
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.a 28 999.86 1086.7 -471.93 943.86 0.0285 1 0.8658964
## nonpsych.b 27 997.89 1081.6 -471.95 943.89 0.1188 1 0.7303385
## nonpsych.c 26 996.01 1076.6 -472.01 944.01 0.2978 1 0.5852428
## nonpsych.d 25 994.31 1071.8 -472.15 944.31 0.6126 1 0.4338281
## nonpsych.e 24 992.92 1067.3 -472.46 944.92 0.1463 1 0.7020980
## nonpsych.f 23 991.07 1062.4 -472.53 945.07 0.4782 1 0.4892392
## nonpsych.g 22 989.55 1057.7 -472.77 945.55 0.6770 5 0.9842072
## nonpsych.h 17 980.22 1032.9 -473.11 946.22 3.1688 1 0.0750589 .
## nonpsych.i 16 981.39 1031.0 -474.70 949.39 3.4825 1 0.0620223 .
## nonpsych.j 15 982.87 1029.4 -476.44 952.87 2.4527 1 0.1173220
## nonpsych.k 14 983.33 1026.7 -477.66 955.33 2.6619 4 0.6158920
## nonpsych.l 10 977.99 1009.0 -478.99 957.99 0.3568 1 0.5503130
## nonpsych.m 9 976.35 1004.2 -479.17 958.35 14.7681 1 0.0001216 ***
## nonpsych.n 8 989.11 1013.9 -486.56 973.11 144.1549 3 < 2.2e-16 ***
## nonpsych.o 5 1127.27 1142.8 -558.63 1117.27
## nonpsych2 29 1001.86 1091.8 -471.93 943.86 0.0010 1 0.9746008
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We use the anova function to compare AIC and model deviance.
anova(nonpsych2, nonpsych.a, nonpsych.b, nonpsych.c, nonpsych.d, nonpsych.e, nonpsych.f, nonpsych.g, nonpsych.h, nonpsych.i, nonpsych.j, nonpsych.k, nonpsych.l, nonpsych.m, nonpsych.n, nonpsych.o)[rev(order(AIC(nonpsych2, nonpsych.a, nonpsych.b, nonpsych.c, nonpsych.d, nonpsych.e, nonpsych.f, nonpsych.g, nonpsych.h, nonpsych.i, nonpsych.j, nonpsych.k, nonpsych.l, nonpsych.m, nonpsych.n, nonpsych.o)$AIC)),]## Data: data
## Models:
## nonpsych.o: cbind(corr, incorr) ~ (1 + Colour.difference | batch) + (1 |
## nonpsych.o: batch)
## nonpsych.n: cbind(corr, incorr) ~ Colour.difference + (1 + Colour.difference |
## nonpsych.n: chick) + (1 + Colour.difference | batch)
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## nonpsych.l: cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## nonpsych.l: Colour.difference | chick) + (1 + Colour.difference | batch)
## nonpsych.k: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.k: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.k: background | batch)
## nonpsych.j: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.j: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.j: background | batch) + target:sex
## nonpsych.i: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.i: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.i: background | batch) + Colour.difference:target + target:sex
## nonpsych.h: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.h: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.h: background | batch) + Colour.difference:target + background:target +
## nonpsych.h: target:sex
## nonpsych.g: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.g: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.g: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.g: background:target + target:sex
## nonpsych.f: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.f: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.f: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.f: background:target + background:sex + target:sex
## nonpsych.e: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.e: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.e: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.e: background:target + Colour.difference:sex + background:sex +
## nonpsych.e: target:sex
## nonpsych.d: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.d: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.d: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.d: background:target + Colour.difference:sex + background:sex +
## nonpsych.d: target:sex + background:target:sex
## nonpsych.c: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.c: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.c: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.c: background:target + Colour.difference:sex + background:sex +
## nonpsych.c: target:sex + Colour.difference:background:target + background:target:sex
## nonpsych.b: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.b: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.b: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.b: background:target + Colour.difference:sex + background:sex +
## nonpsych.b: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.b: background:target:sex
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych2 29 1001.86 1091.8 -471.93 943.86 0.0010 1 0.9746008
## nonpsych.o 5 1127.27 1142.8 -558.63 1117.27
## nonpsych.n 8 989.11 1013.9 -486.56 973.11 144.1549 3 < 2.2e-16 ***
## nonpsych.m 9 976.35 1004.2 -479.17 958.35 14.7681 1 0.0001216 ***
## nonpsych.l 10 977.99 1009.0 -478.99 957.99 0.3568 1 0.5503130
## nonpsych.k 14 983.33 1026.7 -477.66 955.33 2.6619 4 0.6158920
## nonpsych.j 15 982.87 1029.4 -476.44 952.87 2.4527 1 0.1173220
## nonpsych.i 16 981.39 1031.0 -474.70 949.39 3.4825 1 0.0620223 .
## nonpsych.h 17 980.22 1032.9 -473.11 946.22 3.1688 1 0.0750589 .
## nonpsych.a 28 999.86 1086.7 -471.93 943.86 0.0285 1 0.8658964
## nonpsych.d 25 994.31 1071.8 -472.15 944.31 0.6126 1 0.4338281
## nonpsych.e 24 992.92 1067.3 -472.46 944.92 0.1463 1 0.7020980
## nonpsych.f 23 991.07 1062.4 -472.53 945.07 0.4782 1 0.4892392
## nonpsych.g 22 989.55 1057.7 -472.77 945.55 0.6770 5 0.9842072
## nonpsych.c 26 996.01 1076.6 -472.01 944.01 0.2978 1 0.5852428
## nonpsych.b 27 997.89 1081.6 -471.95 943.89 0.1188 1 0.7303385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(nonpsych2, nonpsych.a)## Data: data
## Models:
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.a 28 999.86 1086.7 -471.93 943.86
## nonpsych2 29 1001.86 1091.8 -471.93 943.86 0.001 1 0.974620200420 reviewer’s comments ———————————
Previous model selection procedure: 1. If m1 and m2 are not significantly different, keep m1. 2. If model 1 and m2 are significantly different and m1 has lower AIC, keep m1. 3. If m1 and m2 are significantly different and m2 has lower AIC, choose m2. NEW model selection procedure: 1. If m1 and m2 are not significantly different, choose m2. 2. If m1 and m2 are significantly different and m1 has lower AIC, keep m1. 3. If m1 and m2 are significantly different and m2 has lower AIC, choose m2.
anova(nonpsych2, nonpsych.a)## Data: data
## Models:
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## nonpsych2: cbind(corr, incorr) ~ Colour.difference * background * target *
## nonpsych2: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych2: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.a 28 999.86 1086.7 -471.93 943.86
## nonpsych2 29 1001.86 1091.8 -471.93 943.86 0.001 1 0.9746Rule 1, drop Colour.difference:background:target:sex
anova(nonpsych.a, nonpsych.b)## Data: data
## Models:
## nonpsych.b: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.b: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.b: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.b: background:target + Colour.difference:sex + background:sex +
## nonpsych.b: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.b: background:target:sex
## nonpsych.a: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.a: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.a: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.a: background:target + Colour.difference:sex + background:sex +
## nonpsych.a: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.a: Colour.difference:target:sex + background:target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.b 27 997.89 1081.6 -471.95 943.89
## nonpsych.a 28 999.86 1086.7 -471.93 943.86 0.0285 1 0.8659Rule 1, drop Colour.difference:target:sex
anova(nonpsych.b,nonpsych.c)## Data: data
## Models:
## nonpsych.c: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.c: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.c: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.c: background:target + Colour.difference:sex + background:sex +
## nonpsych.c: target:sex + Colour.difference:background:target + background:target:sex
## nonpsych.b: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.b: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.b: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.b: background:target + Colour.difference:sex + background:sex +
## nonpsych.b: target:sex + Colour.difference:background:target + Colour.difference:background:sex +
## nonpsych.b: background:target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.c 26 996.01 1076.6 -472.01 944.01
## nonpsych.b 27 997.89 1081.6 -471.95 943.89 0.1188 1 0.7303Rule 1, drop Colour.difference:background:sex
anova(nonpsych.c,nonpsych.d)## Data: data
## Models:
## nonpsych.d: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.d: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.d: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.d: background:target + Colour.difference:sex + background:sex +
## nonpsych.d: target:sex + background:target:sex
## nonpsych.c: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.c: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.c: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.c: background:target + Colour.difference:sex + background:sex +
## nonpsych.c: target:sex + Colour.difference:background:target + background:target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.d 25 994.31 1071.8 -472.15 944.31
## nonpsych.c 26 996.01 1076.6 -472.01 944.01 0.2978 1 0.5852Rule 1, drop Colour.difference:background:target
anova(nonpsych.d,nonpsych.e)## Data: data
## Models:
## nonpsych.e: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.e: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.e: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.e: background:target + Colour.difference:sex + background:sex +
## nonpsych.e: target:sex
## nonpsych.d: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.d: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.d: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.d: background:target + Colour.difference:sex + background:sex +
## nonpsych.d: target:sex + background:target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.e 24 992.92 1067.3 -472.46 944.92
## nonpsych.d 25 994.31 1071.8 -472.15 944.31 0.6126 1 0.4338Rule 1, drop background:target:sex
anova(nonpsych.e,nonpsych.f)## Data: data
## Models:
## nonpsych.f: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.f: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.f: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.f: background:target + background:sex + target:sex
## nonpsych.e: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.e: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.e: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.e: background:target + Colour.difference:sex + background:sex +
## nonpsych.e: target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.f 23 991.07 1062.4 -472.53 945.07
## nonpsych.e 24 992.92 1067.3 -472.46 944.92 0.1463 1 0.7021Rule 1, drop Colour.difference:sex
anova(nonpsych.f,nonpsych.g)## Data: data
## Models:
## nonpsych.g: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.g: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.g: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.g: background:target + target:sex
## nonpsych.f: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.f: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.f: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.f: background:target + background:sex + target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.g 22 989.55 1057.7 -472.77 945.55
## nonpsych.f 23 991.07 1062.4 -472.53 945.07 0.4782 1 0.4892Rule 1, drop background:sex
anova(nonpsych.g,nonpsych.h)## Data: data
## Models:
## nonpsych.h: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.h: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.h: background | batch) + Colour.difference:target + background:target +
## nonpsych.h: target:sex
## nonpsych.g: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.g: sex + (1 + Colour.difference | chick) + (1 + Colour.difference *
## nonpsych.g: background | batch) + Colour.difference:background + Colour.difference:target +
## nonpsych.g: background:target + target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.h 17 980.22 1032.9 -473.11 946.22
## nonpsych.g 22 989.55 1057.7 -472.77 945.55 0.677 5 0.9842Rule 1, drop Colour.difference:background
anova(nonpsych.h,nonpsych.i)## Data: data
## Models:
## nonpsych.i: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.i: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.i: background | batch) + Colour.difference:target + target:sex
## nonpsych.h: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.h: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.h: background | batch) + Colour.difference:target + background:target +
## nonpsych.h: target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.i 16 981.39 1031.0 -474.70 949.39
## nonpsych.h 17 980.22 1032.9 -473.11 946.22 3.1688 1 0.07506 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rule 1, drop background:target N.B. AIC now increasing!
anova(nonpsych.i,nonpsych.j)## Data: data
## Models:
## nonpsych.j: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.j: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.j: background | batch) + target:sex
## nonpsych.i: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.i: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.i: background | batch) + Colour.difference:target + target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.j 15 982.87 1029.4 -476.44 952.87
## nonpsych.i 16 981.39 1031.0 -474.70 949.39 3.4825 1 0.06202 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rule 1, I suppose, drop Colour.difference:target N.B. AIC now increasing!
anova(nonpsych.j,nonpsych.k)## Data: data
## Models:
## nonpsych.k: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.k: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.k: background | batch)
## nonpsych.j: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.j: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.j: background | batch) + target:sex
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.k 14 983.33 1026.7 -477.66 955.33
## nonpsych.j 15 982.87 1029.4 -476.44 952.87 2.4527 1 0.1173Rule 1, I suppose, drop target:sex N.B. AIC now increasing!
anova(nonpsych.k,nonpsych.l)## Data: data
## Models:
## nonpsych.l: cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## nonpsych.l: Colour.difference | chick) + (1 + Colour.difference | batch)
## nonpsych.k: cbind(corr, incorr) ~ Colour.difference + background + target +
## nonpsych.k: sex + (1 + Colour.difference | chick) + (1 + Colour.difference +
## nonpsych.k: background | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.l 10 977.99 1009.0 -478.99 957.99
## nonpsych.k 14 983.33 1026.7 -477.66 955.33 2.6619 4 0.6159Rule 1, drop background
anova(nonpsych.l,nonpsych.m)## Data: data
## Models:
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## nonpsych.l: cbind(corr, incorr) ~ Colour.difference + target + sex + (1 +
## nonpsych.l: Colour.difference | chick) + (1 + Colour.difference | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.m 9 976.35 1004.2 -479.17 958.35
## nonpsych.l 10 977.99 1009.0 -478.99 957.99 0.3568 1 0.5503Rule 1, drop sex
anova(nonpsych.m,nonpsych.n)## Data: data
## Models:
## nonpsych.n: cbind(corr, incorr) ~ Colour.difference + (1 + Colour.difference |
## nonpsych.n: chick) + (1 + Colour.difference | batch)
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.n 8 989.11 1013.9 -486.56 973.11
## nonpsych.m 9 976.35 1004.2 -479.17 958.35 14.768 1 0.0001216 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rule 2, keep target
nonpsych.m1 <- update(nonpsych.m, .~. - Colour.difference -
(1 + Colour.difference | chick) -
(1 + Colour.difference | batch)
+ (1 | chick) +
(1 | batch)
)
anova(nonpsych.m,nonpsych.m1)## Data: data
## Models:
## nonpsych.m1: cbind(corr, incorr) ~ target + (1 | chick) + (1 | batch)
## nonpsych.m: cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## nonpsych.m: chick) + (1 + Colour.difference | batch)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## nonpsych.m1 4 1927.51 1939.9 -959.76 1919.51
## nonpsych.m 9 976.35 1004.2 -479.17 958.35 961.17 5 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rule 2, keep Colour.difference So final model, according to reviewer 2
formula(nonpsych.m)## cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## chick) + (1 + Colour.difference | batch)hypothesis testing —————————
ls.comps <- lsmeans(nonpsych.m, list(pairwise ~ target))
summary(ls.comps)[2]## $`pairwise differences of target`
## contrast estimate SE df z.ratio p.value
## diff - same -0.579 0.13 Inf -4.454 <.0001
##
## Results are given on the log odds ratio (not the response) scale.20200806 reviewer’s comments ———————————
L662: what models do you compare. You should compare against the final model, i.e. the model that has only significant terms in it. if you use the full model, you underestimate your probabilities as you underestimate your residual variance. Please make clear what your comparison is.
formula(nonpsych.m)#final model, with only "significant terms"## cbind(corr, incorr) ~ Colour.difference + target + (1 + Colour.difference |
## chick) + (1 + Colour.difference | batch)extract only vital info
ifo <- c("deviance","Chisq","Df","Pr(>Chisq)")terms to test # . sex
anova(nonpsych.m, update(nonpsych.m, .~. +sex))[ifo]## deviance Chisq Df Pr(>Chisq)
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + sex) 957.99 0.3568 1 0.5503. background
anova(nonpsych.m, update(nonpsych.m, .~. +background))[ifo]## deviance Chisq Df Pr(>Chisq)
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + background) 958.32 0.0292 1 0.8644. same:sex
anova(nonpsych.m, update(nonpsych.m, .~. +target*sex))[ifo]## deviance Chisq Df Pr(>Chisq)
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + target * sex) 955.99 2.3526 2 0.3084. Colour.difference:target
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*target))[ifo]## deviance Chisq Df
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + Colour.difference * target) 955.23 3.1199 1
## Pr(>Chisq)
## nonpsych.m
## update(nonpsych.m, . ~ . + Colour.difference * target) 0.07734 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1. background:target
anova(nonpsych.m, update(nonpsych.m, .~. +background*target))[ifo]## deviance Chisq Df Pr(>Chisq)
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + background * target) 956.35 1.9927 2 0.3692. Colour.difference:background
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*background))[ifo]## deviance Chisq Df
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + Colour.difference * background) 957.68 0.6674 2
## Pr(>Chisq)
## nonpsych.m
## update(nonpsych.m, . ~ . + Colour.difference * background) 0.7163. background:sex
anova(nonpsych.m, update(nonpsych.m, .~. +background*sex))[ifo]## deviance Chisq Df Pr(>Chisq)
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + background * sex) 957.30 1.0469 3 0.7899. Colour.difference:sex
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*sex))[ifo]## deviance Chisq Df
## nonpsych.m 958.35
## update(nonpsych.m, . ~ . + Colour.difference * sex) 957.86 0.4841 2
## Pr(>Chisq)
## nonpsych.m
## update(nonpsych.m, . ~ . + Colour.difference * sex) 0.785. background:target:sex
anova(nonpsych.m, update(nonpsych.m, .~. +background*target*sex)
)[2,ifo]## deviance Chisq Df
## update(nonpsych.m, . ~ . + background * target * sex) 952.28 6.0654 6
## Pr(>Chisq)
## update(nonpsych.m, . ~ . + background * target * sex) 0.4159. Colour.difference:background:target
unlist(
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*background*target)
)[2,ifo]
)## deviance Chisq Df Pr(>Chisq)
## 952.0507817 6.2954922 5.0000000 0.2785187. Colour.difference:background:sex
unlist(
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*background*sex)
)[2,ifo]
)## deviance Chisq Df Pr(>Chisq)
## 956.3960780 1.9501959 6.0000000 0.9242216. Colour.difference:target:sex
unlist(
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*target*sex))[2,ifo]
)## deviance Chisq Df Pr(>Chisq)
## 952.3884536 5.9578203 5.0000000 0.3103448. Colour.difference:background:same:sex
unlist(
anova(nonpsych.m, update(nonpsych.m, .~. +Colour.difference*background*target*sex))[2,ifo]
)## deviance Chisq Df Pr(>Chisq)
## 947.5137150 10.8325589 13.0000000 0.6248424Setting up the fitted function
xseq <- unique(data$Colour.difference) # imaginary stimulus levels to fit to. Change the high and low values, the two values after seq(, to fit to your datas stimulus levels
newdata <- expand.grid(Colour.difference = xseq,
background = unique(data$background),
target = unique(data$target),
sex = unique(data$sex),
chick = unique(data$chick),
batch = unique(data$batch))
nonpsych1.pred <- predict(nonpsych2,
newdata = data.frame(Colour.difference = data$Colour.difference,
background = data$background,
target = data$target, sex = data$sex),
type="response",
re.form = NA)# ##Setting up the real data from the function
nonpsych2.pred <- predict(nonpsych2, newdata = newdata, type="response")#,se.fit=TRUE) ## Setting up the imaginary data from the function
plot(data$Colour.difference, nonpsych1.pred)
Setting up the plot
bktg <- data.frame(bk = levels(data$background)[c(1,2,2,1)],
tg = rev(levels(newdata$target))[c(1,2,1,2)])
for(i in 1:dim(bktg)[1]){
bk <- bktg$bk[i]
tg <- bktg$tg[i]
plot(NULL,
xlab="Colour difference",
ylab="Proportion correct",
xlim=c(0, 6),
ylim=c(0.4, 1),
main = paste('background =',bk, ', target =',tg)
)
if(bk == 'green' & tg == 'same'){
legend('bottomright',
legend = c('Male', 'Female', paste('Batch',c('A','B','C','D'))),
col = c(2*(2:1), rep(1, 4)),
lty = c(1,1,1:4),
pch = c(24,22,rep(NA,4)),
cex = 0.7)
}
for(ck in levels(newdata$chick)){
if(sum(data$background == bk & data$chick == ck & data$target == tg)){
sx <- unique(subset(data, chick == ck)$sex)
bc <- unique(subset(data, chick == ck)$batch)
xx <- subset(newdata,
background == bk & chick == ck &
target == tg & sex == sx & batch == bc
)
yy <- nonpsych2.pred[
newdata$background == bk & newdata$chick == ck &
newdata$target == tg & newdata$sex == sx & newdata$batch == bc
]
ORDER <- order(xx$Colour.difference)
lines(xx$Colour.difference[ORDER],
yy[ORDER],
col = 2*which(unique(subset(data, chick == ck)$sex) == levels(newdata$sex)),
lty = which(unique(subset(data, chick == ck)$batch) == levels(newdata$batch)),
lwd = 0.5
)
}
}
points(data$Colour.difference[data$background == bk & data$target == tg],
data$pcorr[data$background == bk & data$target == tg],
pch = c(22,24)[as.numeric(data$sex[data$background == bk & data$target == tg])],
bg = 'white',
col = rgb(0,0,0,0.3),
lwd = 2)
}
The following chunk defines the colour series.
cls<-c("purple4","slateblue3","slateblue2","red3","green3",
"slateblue1","pink3","orange3","navajowhite4","gray50",
"gray70","gray30","darkblue","navajowhite2","orange4",
"steelblue","gray10","purple3","magenta4","slateblue4",
"green2","blue2","darkred","darkgreen","orange2",
"seagreen","salmon4","navajowhite1","navajowhite3","yellow3",
"blue3","magenta3")#shorten as.factor for use in level sorting #N.B. unique(x) may have been more efficient than levels(as.factor(x))
AF <- function(x){as.factor(x)}#panels to plot in vertical and horizontal direction
hw <- c(2,2)
par(mfrow = c(hw), mai = .75*c(.8,1,.5,0))
bktg <- data.frame(bk = levels(data$background)[c(1,2,2,1)],
tg = rev(levels(newdata$target))[c(1,2,1,2)]
)
for(i in 1:dim(bktg)[1]){
bk <- bktg$bk[i]
tg <- bktg$tg[i]
plot(NULL,
xlab="Colour difference",
ylab="Proportion correct",
xlim=c(0, 6),
ylim=c(0.4, 1),
main = paste('background =',bk, ', target =',tg)
)
legend(2,0.7,#4, 0.7,
levels(AF(data$chick[data$target == tg & data$background == bk])),
col = cls[
which( levels(AF(data$chick)) %in%
levels(AF(data$chick[data$target == tg & data$background == bk]))
)
],
pch = c(22,24)[
as.numeric(
data$sex[data$Colour.difference == max(data$Colour.difference[data$chick %in%
levels(AF(data$chick[data$target == tg & data$background == bk]))]) &
data$chick %in%
levels(AF(data$chick[data$target == tg & data$background == bk]))]
)
],
cex = 1*0.3, bty = 'n', lwd = 1,
lty = as.numeric(
subset(data,
chick %in%
levels(AF(data$chick[data$target == tg & data$background == bk])) &
data$Colour.difference ==
max(data$Colour.difference[
data$chick %in% levels(AF(data$chick[data$target == tg &
data$background == bk]))]
)
)$batch
)
)
if(bk == 'green' & tg == 'diff'){
legend('bottomright',
legend = c('Male', 'Female', paste('Batch',c('A','B','C','D'))),
col = c(2*(2:1), rep(1, 4)),
lty = c(1,1,1:4),
pch = c(24,22,rep(NA,4)),
cex = 0.7
)
}
for(ck in levels(newdata$chick)){
if(sum(data$background == bk & data$chick == ck & data$target == tg)){
sx <- unique(subset(data, chick == ck)$sex)
bc <- unique(subset(data, chick == ck)$batch)
xx <- subset(newdata,
background == bk & chick == ck &
target == tg & sex == sx & batch == bc
)
yy <- nonpsych2.pred[
newdata$background == bk & newdata$chick == ck &
newdata$target == tg & newdata$sex == sx & newdata$batch == bc
]
ORDER <- order(xx$Colour.difference)
lines(xx$Colour.difference[ORDER],
yy[ORDER],
col = cls[which(levels(AF(data$chick)) == ck)],
lty = which(unique(subset(data, chick == ck)$batch) == levels(newdata$batch)),
lwd = 0.5 )
points(data$Colour.difference[
data$background == bk & data$chick == ck &
data$target == tg & data$sex == sx &
data$batch == bc
],
data$pcorr[
data$background == bk & data$chick == ck &
data$target == tg & data$sex == sx & data$batch == bc
],
pch = c(22,24)[
as.numeric(
data$sex[
data$background == bk & data$target == tg & data$chick == ck
]
)
],
bg = 'white',
col = cls[which(levels(AF(data$chick)) == ck)],
lwd = 2)
}
}
}
ndata <- expand.grid(Colour.difference =xseq,
background=unique(data$background),
target=unique(data$target))
prd <- predict(nonpsych2, newdata = newdata, type="response", re.form = NA)
pfun <- function(x){predict(x, newdata = newdata, type = 'response')}Bootstrap model confidence intervals
if(Sys.info()[['sysname']] == 'Windows')
{
require(snow); clt <- makeCluster(parallel::detectCores(), type = 'SOCK')
clusterExport(clt,list('nonpsych2','pfun','newdata'))
boot <- bootMer(nonpsych2, pfun, nsim = 10^2, re.form = NA,
parallel = c("snow"), ncpus = parallel::detectCores(), cl = clt)
stopCluster(clt)
}else
{
boot <- bootMer(nonpsych2, pfun, nsim = 10^2, re.form = NA,
parallel = c("multicore"), ncpus = parallel::detectCores())
}Calculate confidence intervals from bootstrapped model
stqlog <- function(x){ sd(qlogis(x), na.rm = T)}
std.err <- apply(boot$t, 2, stqlog)#standard error for each
CI.lo <- plogis(qlogis(prd) - std.err*1.96)#lower confidence bound (parametric)
CI.hi <- plogis(qlogis(prd) + std.err*1.96)#upper confidence bound (parametric)
bootsdata <- cbind(newdata, CI.lo, CI.hi)
CI.lo. <- with(bootsdata, aggregate(CI.lo, list(Colour.difference = Colour.difference, background = background, target = target), mean))
CI.hi. <- with(bootsdata, aggregate(CI.hi, list(Colour.difference = Colour.difference, background = background, target = target), mean))Also check bootstrapped parameter estimates
ciii <- confint(boot, parallel = c("multicore"),ncpus = parallel::detectCores())
paramdata <- cbind(newdata, ciii)
#mean should be in logit space
lgtmean <- function(x){plogis(mean(qlogis(x)))}
CI.02.5 <- with(paramdata, aggregate(`2.5 %`, list(Colour.difference = Colour.difference, background = background, target = target), lgtmean))#mean))
CI.97.5 <- with(paramdata, aggregate(`97.5 %`, list(Colour.difference = Colour.difference, background = background, target = target), lgtmean))#mean))Find sensible averages for each condition
estdata <- cbind(newdata, prd)[order(newdata$Colour.difference),]
estcurve <- with(estdata, aggregate(prd, list(Colour.difference = Colour.difference, background = background, target = target), lgtmean))
threshquant <- list(same = (0:20/100), diff = (0:40/100))hw <- c(2,2)
par(mfrow = c(hw), mai = .75*c(.8,1,.5,0))
bktg <- data.frame(bk = levels(data$background)[c(1,2,2,1)],
tg = rev(levels(newdata$target))[c(1,2,1,2)])
for(i in 1:dim(bktg)[1]){
bk <- bktg$bk[i]
tg <- bktg$tg[i]
plot(NULL, xlab="Colour difference", ylab="Proportion correct", xlim=c(0, 6), ylim=c(0.4, 1), main = paste('background =',bk, ', target =', tg))
xx <- subset(estcurve, background == bk & target == tg)$Colour.difference
yy <- subset(estcurve, background == bk & target == tg)$x
threshness <-abs(yy - 0.65)
threshx <- mean(xx[threshness %in% quantile(threshness, unlist(threshquant[2 - i %% 2])) ])
print(threshx)
polygon(c(subset(CI.02.5, background == bk & ndata$target == tg)$Colour.difference, rev(subset(CI.97.5, background == bk & ndata$target == tg)$Colour.difference)), c(subset(CI.02.5, background == bk & ndata$target == tg)$x,rev(subset(CI.97.5, background == bk & ndata$target == tg)$x)) , col = 'gray', border = rgb(0,0,0,0))
# lines(xx[ORDER], yy[ORDER], col = 'black', lty = 1, lwd = 5)
lines(xx, yy, col = 'black', lty = 1, lwd = 1)
points(data$Colour.difference[data$background == bk & data$target == tg], data$pcorr[data$background == bk & data$target == tg], pch = c(22,24)[as.numeric(data$sex[data$background == bk & data$target == tg])], bg = 'white', col = rgb(0,0,0,0.3), lwd = 2)
lines(rep(threshx, 2), c(0.36, 0.65), lty = 2)
lines(c(-0.25,threshx), rep(0.65,2), lty = 2)
}## [1] 0.5333333## [1] 0.96## [1] 0.8666667## [1] 0.96legend('bottomright', legend = c('Male', 'Female'),
col = c(1,1), pch = c(24,22), lwd = c(2,2))
std.err <- apply(boot$t, 2, stqlog)#standard error for each
CI.lo.se <- plogis(qlogis(prd) - std.err)#lower se bound (parametric)
CI.hi.se <- plogis(qlogis(prd) + std.err)#upper se bound (parametric)
bootsdata.se <- cbind(newdata, CI.lo.se, CI.hi.se)
CI.lo.se. <- with(bootsdata.se, aggregate(CI.lo.se, list(Colour.difference = Colour.difference, background = background, target = target), lgtmean))
CI.hi.se. <- with(bootsdata.se, aggregate(CI.hi.se, list(Colour.difference = Colour.difference, background = background, target = target), lgtmean))Make a function to evaluate the threshold.
thresholder <- function(xx,yy,lev){
diff. <- sort(yy-lev)
close. <- min(abs(diff.))
if(sign(diff.[abs(diff.) == close.]) == 0){
return( xx[which(abs(diff.) == close.) + c(0)] )
}else{
if(sign(diff.[abs(diff.) == close.]) == 1){
if(which(abs(diff.) == close.) == 1){
return( xx[1])
}else{
ab. <- xx[which(abs(diff.) == close.) + c(-1,0)]
}
}else{
ab. <- xx[which(abs(diff.) == close.) + c(0,1)]
}
ty. <- yy[xx %in% ab.]
xxt. <- seq(min(ab.), max(ab.), length.out = 10^3)
yyt. <- (xxt.-min(ab.))*diff(ty.)/diff(ab.) + min(ty.)
return( mean(xxt.[round(yyt.,2) == lev]) )
}
}par(mfrow = c(hw), mai = .75*c(.8,1,.5,0))
for(i in 1:dim(bktg)[1]){
bk <- bktg$bk[i]
tg <- bktg$tg[i]
plot(NULL, xlab="Colour difference", ylab="Proportion correct",
xlim=c(0, 6), ylim=c(0.4, 1),
main = paste('background =',bk, ', target =', tg))
abline(h = 0.65, lty = 2, lwd = 0.5)
x1 <- subset(estcurve, background == bk & target == tg)$Colour.difference
y1 <- subset(estcurve, background == bk & target == tg)$x
xx <- sort(unique(subset(CI.lo.se., background == bk & target==tg)$Colour.difference ))
yy.lo <- sort(unique( subset(CI.lo.se., background == bk & target == tg)$x ))
yy.hi <- sort(unique( subset(CI.hi.se., background == bk & target == tg)$x ))
lines(xx, yy.lo)
lines(xx, yy.hi)
lines(x1, y1, col = 'darkgreen')
assign(paste0('t.', tg, bk), thresholder(x1,y1,0.65))
lines(rep(get(paste0('t.', tg, bk)),2), c(0.38, 0.65), lty = 2, col = 'green')
assign(paste0('t.se.lo.', tg, bk), thresholder(xx,yy.hi,0.65))
lines(rep(get(paste0('t.se.lo.', tg, bk)),2), c(0.38, 0.65), lty = 2,col ='red')
assign(paste0('t.se.hi.', tg, bk), thresholder(xx,yy.lo,0.65))
lines(rep(get(paste0('t.se.hi.', tg, bk)),2), c(0.38, 0.65), lty = 2,col ='blue')
}
The “binomial” 65% threshold
JND for chickens discriminating orange colours on an orange background
diff(round( c(
t.se.lo.sameorange,
t.se.hi.sameorange ), 2)
)## [1] 0.52JND for chickens discriminating orange colours on a green background
diff(round( c(
t.se.lo.diffgreen,
t.se.hi.diffgreen ), 2)
)## [1] 0.53JND for chickens discriminating green colours on an orange background.
diff(round( c(
t.se.lo.difforange,
t.se.hi.difforange ), 2)
)## [1] 0.8JND for chickens discriminating green colours on an green background.
diff(round( c(
t.se.lo.samegreen,
t.se.hi.samegreen ), 2)
)## [1] 0.48