Accuracy Analyses for Categorization Task
Libraries and Data Files
#load libraries
library(dplyr)
library(sciplot)
library(plotrix)
library(emmeans)
library(mgcv)
library(itsadug)
#Load data
data<- read.csv("AllData_Categorization.csv", header = TRUE)
Data Preprocessing
# Rename useful variables
data$acc <- data$Accuracy
data$sbj <- data$participant
data$block <- data$block.thisN + 1
# Convert variables to factors
data <- mutate_if(data, is.character, as.factor)
data$sbj <- as.factor(data$sbj)
data$group <- as.factor(data$group)
# Exclude only the participant who did not follow instructions
data <- droplevels(data[data$sbj != "20335", ])
# Create trial number
data$trial <- data$block.thisN * 80 + data$trials.thisN + 1
# Check number of participants per group
data %>%
distinct(sbj, group) %>%
count(group)
## group n
## 1 1100 15
## 2 1600 17
## 3 600 18
## 4 RD 16
RD Exposure Duration
rt_by_sbj_RD <- aggregate(RT ~ sbj, data = data[data$group == "RD", ], FUN = mean, na.rm = TRUE)
#descriptives
mean(rt_by_sbj_RD$RT)
## [1] 0.844926
sd(rt_by_sbj_RD$RT)
## [1] 0.210926
t.test(rt_by_sbj_RD$RT)$conf.int
## [1] 0.7325314 0.9573205
## attr(,"conf.level")
## [1] 0.95
Descriptive Statistics
#calculate participant's average accuracy
data_av<-aggregate(data$acc, list(data$sbj, data$group), mean)
colnames(data_av) <-c("sbj","group" ,"acc")
#write data frame, to use it in analysis examining correlation between learning accuracy and perceptual change
write.csv(data_av, "d_cat.csv", row.names = FALSE)
#by-group accuracy
#"600"
mean(data_av[data_av$group=="600",]$acc)
## [1] 0.8303819
sd(data_av[data_av$group=="600",]$acc)
## [1] 0.08761026
#"1100"
mean(data_av[data_av$group=="1100",]$acc)
## [1] 0.8670833
sd(data_av[data_av$group=="1100",]$acc)
## [1] 0.03201882
"#1600"
## [1] "#1600"
mean(data_av[data_av$group=="1600",]$acc)
## [1] 0.8689338
sd(data_av[data_av$group=="1600",]$acc)
## [1] 0.04216537
#"RD"
mean(data_av[data_av$group=="RD",]$acc)
## [1] 0.8714844
sd(data_av[data_av$group=="RD",]$acc)
## [1] 0.03289177
Learning-Curves Plot
##################
# Learning Curves
##################
bins <- 16 #thus, 16 bins of 20 trials (320 trials overall)
data <- data %>%
group_by(sbj) %>%
mutate(trial_bin = ntile(trial, bins)) %>%
ungroup()
data_tb<-aggregate(data$acc, list(data$sbj, data$group, data$trial_bin), mean)
colnames(data_tb) <-c("sbj","group","trial_bin","acc")
data_gr<-aggregate(data_tb$acc, list(data_tb$group, data_tb$trial_bin), mean)
colnames(data_gr)<-c("group","trial_bin", "acc")
temp<- aggregate(data_tb$acc, list(data_tb$group, data_tb$trial_bin), se)
data_gr$se<-temp$x
rm(temp)
data_gr<-data_gr[order(data_gr$group, data_gr$trial_bin),]
# Labels
xlb="Bins of 20 Trials"
ylb="Proportion Correct"
mn="Categorization Accuracy"
off=0.10
pch=c(15,16,17,18)
lty=c(2,3,4,5)
col=c("grey70","grey60","grey50", "grey40" )
x0 <- 1:bins
#"600"
plotCI(x=x0, y=data_gr[data_gr$group=="600",]$acc,
uiw=data_gr[data_gr$group=="600",]$se,
bty="n", ylim=c(0.5,1),xlim=c(1, bins+0.6),
xlab=xlb, ylab=ylb, main=mn, las=1, xaxt="n", yaxt="n",
pch=pch[1], col=col[1], gap=0, cex=1.5)
lines(x=x0, y=data_gr[data_gr$group=="600",]$acc, lty=lty[1], col=col[1], lwd=2)
#"1100"
plotCI(x=x0+off, y=data_gr[data_gr$group=="1100",]$acc,
uiw=data_gr[data_gr$group=="1100",]$se,
add=T, pch=pch[2], col=col[2], gap=0, cex=1.5)
lines(x=x0+off, y=data_gr[data_gr$group=="1100",]$acc, lty=lty[2], col=col[2], lwd=2)
#"1600"
plotCI(x=x0+2*off, y=data_gr[data_gr$group=="1600",]$acc,
uiw=data_gr[data_gr$group=="1600",]$se,
add=T, pch=pch[3], col=col[3], gap=0, cex=1.5)
lines(x=x0+2*off, y=data_gr[data_gr$group=="1600",]$acc, lty=lty[3], col=col[3], lwd=2)
#"RD"
plotCI(x=x0+3*off, y=data_gr[data_gr$group=="RD",]$acc,
uiw=data_gr[data_gr$group=="RD",]$se,
add=T, pch=pch[4], col=col[4], gap=0, cex=1.5)
lines(x=x0+3*off, y=data_gr[data_gr$group=="RD",]$acc, lty=lty[4], col=col[4], lwd=2)
axis(side=1, at=1:bins, labels=1:bins)
axis(side=2, at=c(0.5, 0.6, 0.7, 0.8, 0.9, 1.0),labels=c("0.0", "0.6", "0.7", "0.8", "0.9", "1.0") , las=2)
axis.break(2, 0.55, style="slash")
legend(x=round(bins*0.60), y=0.8, legend=c("600 ms", "1100 ms", "1600 ms", "RD"),
col=col, lty=lty, pch=pch, bty="n", seg.len=4, lwd=2, pt.cex=1.5)
### pixels: 1000 x 650
GAMMs - Inferential Analyses
# Set group order
data$group <- factor(
data$group,
levels = c("600", "1100", "1600", "RD")
)
# Standardize trial number
data$trial_s <- as.numeric(scale(data$trial))
# Create factor for difference smooths
data$group_diff <- ordered(data$group, levels = c("600", "1100", "1600", "RD"))
###############################
# Learning-curve model comparison
###############################
# Fit the shared-curve model using ML.
# ML is used here only to compare this model with a model
# allowing group-specific deviations in the learning curve.
m1_ML <- bam(
acc ~ group +
s(trial_s, k = 10) +
s(trial_s, sbj, bs = "fs", m = 1, k = 10),
data = data,
family = binomial(link = "logit"),
method = "ML",
discrete = FALSE
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
summary(m1_ML)
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.66701 0.08651 19.271 <2e-16 ***
## group1100 0.25303 0.12929 1.957 0.0503 .
## group1600 0.27525 0.12515 2.199 0.0279 *
## groupRD 0.28667 0.12719 2.254 0.0242 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.948 8.69 149.8 <2e-16 ***
## s(trial_s,sbj) 126.064 656.00 487.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0413 Deviance explained = 4.79%
## -ML = 8380.1 Scale est. = 1 n = 21120
###############################
# Group-specific curve model: ML estimation
###############################
# Fit the group-specific curve model using ML.
# This model will be compared with m1_ML to determine whether
# group-specific deviations improve the learning-curve model.
m2_ML <- bam(
acc ~ group +
s(trial_s, k = 10) +
s(trial_s, by = group_diff, k = 10) +
s(trial_s, sbj, bs = "fs", m = 1, k = 10),
data = data,
family = binomial(link = "logit"),
method = "ML",
discrete = FALSE
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
summary(m2_ML)
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, by = group_diff,
## k = 10) + s(trial_s, sbj, bs = "fs", m = 1, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.66642 0.08654 19.256 <2e-16 ***
## group1100 0.25527 0.12940 1.973 0.0485 *
## group1600 0.27784 0.12526 2.218 0.0265 *
## groupRD 0.28542 0.12723 2.243 0.0249 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.949 8.691 127.946 <2e-16 ***
## s(trial_s):group_diff1100 1.005 1.010 0.355 0.560
## s(trial_s):group_diff1600 1.004 1.007 0.407 0.526
## s(trial_s):group_diffRD 1.006 1.010 0.422 0.521
## s(trial_s,sbj) 122.164 656.000 481.029 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0411 Deviance explained = 4.77%
## -ML = 8379.1 Scale est. = 1 n = 21120
AIC(m1_ML, m2_ML)
## df AIC
## m1_ML 140.6691 16667.71
## m2_ML 140.0348 16670.12
delta_AIC <- AIC(m2_ML) - AIC(m1_ML)
delta_AIC
## [1] 2.410365
anova(m1_ML, m2_ML, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
## Model 2: acc ~ group + s(trial_s, k = 10) + s(trial_s, by = group_diff,
## k = 10) + s(trial_s, sbj, bs = "fs", m = 1, k = 10)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 20920 16386
## 2 20922 16390 -2.3279 -3.6789 0.2017
# Model 1 had modestly lower AIC than Model 2 (Delta AIC = 2.41).
# The addition of group-specific smooth deviations was therefore
# not sufficiently supported to prefer the more complex structure.
# Model 1 was retained as the more parsimonious model.
###############################
# Final model: shared learning curve
###############################
m1_fREML <- bam(
acc ~ group +
s(trial_s, k = 10) +
s(trial_s, sbj, bs = "fs", m = 1, k = 10),
data = data,
family = binomial(link = "logit"),
method = "fREML",
discrete = TRUE,
nthreads = 10
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
summary(m1_fREML)
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.66809 0.08913 18.716 <2e-16 ***
## group1100 0.25297 0.13315 1.900 0.0575 .
## group1600 0.27513 0.12889 2.135 0.0328 *
## groupRD 0.28667 0.13098 2.189 0.0286 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.983 8.708 150.3 <2e-16 ***
## s(trial_s,sbj) 126.796 656.000 490.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0413 Deviance explained = 4.8%
## fREML = 8390.5 Scale est. = 1 n = 21120
# Check model
set.seed(1234)
gam.check(m1_fREML,k.rep = 1000)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -1.846125e-07 -2.521784e-05 -4.960516e-06
##
## $hess
## [,1] [,2] [,3]
## [1,] 2.56007093 -0.05002208 -0.03164658
## [2,] -0.05002208 11.58980420 0.08810668
## [3,] -0.03164658 0.08810668 21.92688559
##
## Model rank = 673 / 673
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(trial_s) 9.00 7.98 0.94 0.31
## s(trial_s,sbj) 660.00 126.80 0.94 0.32
###############################
# Primary contrast
###############################
group_emm <- emmeans(m1_fREML, ~ group, at = list(trial_s = 0))
planned_contrast <- contrast( group_emm, method = list("Longer - 600" = c(-1, 1/3, 1/3, 1/3)), adjust = "none")
summary(planned_contrast,infer = c(TRUE, TRUE), type = "response",df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## Longer / 600 1.31 0.138 Inf 1.07 1.61 1 2.588 0.0097
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scale
###############################
# Fixed-duration robustness contrast
###############################
fixed_contrast <- contrast(group_emm, method = list("Fixed longer - 600" = c(-1, 1/2, 1/2, 0)), adjust = "none")
summary(fixed_contrast, infer = c(TRUE, TRUE), type = "response", df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio
## Fixed longer / 600 1.3 0.146 Inf 1.05 1.62 1 2.356
## p.value
## 0.0184
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scale
###############################
# Supplementary omnibus group test
###############################
# Test whether the four groups differ overall.
omnibus_group <- anova(m1_fREML, freq = FALSE)
omnibus_group
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
##
## Parametric Terms:
## df Chi.sq p-value
## group 3 6.784 0.0791
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.983 8.708 150.3 <2e-16
## s(trial_s,sbj) 126.796 656.000 490.8 <2e-16
##############################
# Secondary comparisons with the 600-ms group
###############################
# Compare each longer viewing regime with the 600-ms reference group.
# The three comparisons constitute one family of tests.
dunnett_group <- contrast(group_emm, method = "trt.vs.ctrl", ref = "600")
# Report log-odds differences with multiplicity-adjusted
# confidence intervals and p values.
set.seed(1234)
summary(dunnett_group, infer = c(TRUE, TRUE), adjust = "mvt", df = Inf)
## contrast estimate SE df asymp.LCL asymp.UCL z.ratio p.value
## 1100 - 600 0.253 0.133 Inf -0.0607 0.567 1.900 0.1456
## 1600 - 600 0.275 0.129 Inf -0.0284 0.579 2.135 0.0860
## RD - 600 0.287 0.131 Inf -0.0218 0.595 2.189 0.0756
##
## Degrees-of-freedom method: user-specified
## Results are given on the log odds ratio (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 3 estimates
## P value adjustment: mvt method for 3 tests
# Report the same comparisons as odds ratios.
set.seed(1234)
summary(dunnett_group, infer = c(TRUE, TRUE), type = "response", adjust = "mvt", df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## 1100 / 600 1.29 0.171 Inf 0.941 1.76 1 1.900 0.1456
## 1600 / 600 1.32 0.170 Inf 0.972 1.78 1 2.135 0.0860
## RD / 600 1.33 0.174 Inf 0.978 1.81 1 2.189 0.0756
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 3 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: mvt method for 3 tests
## Tests are performed on the log odds ratio scale
Model-Predictions Plot
fm <- m1_fREML
labx="Trial"
#labx="Trial\n Error Bands: 95% CI"
main="Categorization Accuracy - Model Predictions"
col=c("grey70","grey60","grey50", "grey40" )
lty=c(2,3,4,5)
alpha_shade <- 0.4
# draw band (transparent) + redraw line (opaque)
plot_group <- function(g, i, first = FALSE) {
col_band <- adjustcolor(col[i], alpha.f = alpha_shade)
# 1) band (and faint line)
plot_smooth(
fm, view="trial_s", cond=list(group=g),
rm.ranef=T,
shade=TRUE, se=1.96,
print.summary=FALSE,
ylim=c(0.5,1),
transform=plogis,
xlab=if (first) labx else "",
ylab=if (first) "Proportion Correct" else "",
lwd=2, las=2, rug=FALSE,
col=col_band, lty=lty[i],
hide.label=TRUE,
main=if (first) main else "",
yaxs="i", xaxt=if (first) "n" else "n", yaxt=if (first) "n" else "n",
add = !first
)
# 2) line only (opaque)
plot_smooth(
fm, view="trial_s", cond=list(group=g),
rm.ranef=T,
shade=FALSE, se=0,
print.summary=FALSE,
transform=plogis,
lwd=2, rug=FALSE,
col=col[i], lty=lty[i],
hide.label=TRUE,
add=TRUE
)
}
plot_group("600", 1, first=TRUE)
plot_group("1100", 2)
plot_group("1600", 3)
plot_group("RD", 4)
mu <- attr(scale(data$trial), "scaled:center")
sd <- attr(scale(data$trial), "scaled:scale")
ticks_raw <- c(0,50,100,150,200,250,300,320)
ticks_s <- (ticks_raw - mu) / sd
axis(side=1, at=ticks_s, labels=ticks_raw)
axis(side=2, at=c(0.5, 0.6, 0.7, 0.8, 0.9, 1.0),labels=c("0.0", "0.6", "0.7", "0.8", "0.9", "1.0") , las=2)
axis.break(2, 0.55, style="slash")
legend(x=(200-mu)/sd, y=0.75, legend=c("600 ms", "1100 ms", "1600 ms", "RD"), col=col, lty=lty, bty="n", seg.len=4, lwd=2, pt.cex=1.5)
# graph for paper 1000 x 650
Common Sample - Sensitivity Analysis
###############################
# Common-sample exclusions
###############################
# Participants excluded for not meeting the learning criterion.
accuracy_exclusions <- c("43423", "48888", "63018", "77337", "77623")
# Learning-criterion failure rates across groups
learning_status <- data %>%
distinct(sbj, group) %>%
mutate(nonlearner = as.character(sbj) %in% accuracy_exclusions)
# Rows: groups; columns: learners (FALSE), nonlearners (TRUE)
learning_table <- with(learning_status, table(group, nonlearner))
learning_table
## nonlearner
## group FALSE TRUE
## 600 16 2
## 1100 15 0
## 1600 15 2
## RD 15 1
# Percentage of nonlearners within each group
round(100 * prop.table(learning_table, margin = 1)[, "TRUE"], 1)
## 600 1100 1600 RD
## 11.1 0.0 11.8 6.2
# Fisher's exact test for the 4 × 2 table
fisher.test(learning_table, simulate.p.value = FALSE)
##
## Fisher's Exact Test for Count Data
##
## data: learning_table
## p-value = 0.7453
## alternative hypothesis: two.sided
# Participants excluded because they received the initial
# similarity-rating instructions.
instruction_exclusions <- c("17202", "57393", "94372", "40242", "89443")
data_common <- data %>% filter(!as.character(sbj) %in% c(accuracy_exclusions, instruction_exclusions)) %>% droplevels()
###############################
# Common-sample sensitivity model
###############################
# Fit the selected shared-curve model to the common sample.
# The model structure is held constant relative to the primary analysis.
m1_common_fREML <- bam(acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs", m = 1, k = 10), data = data_common, family = binomial(link = "logit"), method = "fREML", discrete = TRUE, nthreads = 10)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
## 1-d smooths of same variable.
summary(m1_common_fREML)
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.74073 0.06940 25.082 < 2e-16 ***
## group1100 0.16044 0.09915 1.618 0.10561
## group1600 0.30010 0.10020 2.995 0.00274 **
## groupRD 0.25104 0.09979 2.516 0.01188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.945 8.685 131.8 <2e-16 ***
## s(trial_s,sbj) 102.756 556.000 219.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0277 Deviance explained = 3.61%
## fREML = 6921.1 Scale est. = 1 n = 17920
# Check convergence and whether the basis dimensions are adequate.
set.seed(1234)
gam.check(m1_common_fREML, k.rep = 1000)
##
## Method: fREML Optimizer: perf chol
## $grad
## [1] -2.932397e-08 -3.376176e-06 -1.117799e-06
##
## $hess
## [,1] [,2] [,3]
## [1,] 2.46715363 -0.05305422 -0.0175836
## [2,] -0.05305422 11.56434835 0.4206303
## [3,] -0.01758360 0.42063033 9.3625120
##
## Model rank = 573 / 573
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(trial_s) 9.00 7.94 0.96 0.57
## s(trial_s,sbj) 560.00 102.76 0.96 0.58
###############################
# Common-sample contrast
###############################
# Estimate group means at the midpoint of the experiment.
group_emm_common <- emmeans(m1_common_fREML, ~ group, at = list(trial_s = 0))
# Compare the 600-ms group with the average of the three longer viewing regimes.
planned_contrast_common <- contrast(group_emm_common, method = list("Longer - 600" = c(-1, 1/3, 1/3, 1/3)), adjust = "none")
# Report the contrast as an odds ratio with its confidence interval.
summary(planned_contrast_common, infer = c(TRUE, TRUE), type = "response", df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## Longer / 600 1.27 0.102 Inf 1.08 1.49 1 2.937 0.0033
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scale
###############################
# Common-sample fixed-duration robustness contrast
###############################
# Test the same fixed-duration contrast in the common sample.
fixed_contrast_common <- contrast(
group_emm_common,
method = list("Fixed longer - 600" = c(-1, 1/2, 1/2, 0)),
adjust = "none"
)
summary(fixed_contrast_common, infer = c(TRUE, TRUE), type = "response", df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio
## Fixed longer / 600 1.26 0.108 Inf 1.06 1.49 1 2.682
## p.value
## 0.0073
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Intervals are back-transformed from the log odds ratio scale
## Tests are performed on the log odds ratio scale
###############################
# Common-sample omnibus group test
###############################
# Test whether the four groups differ overall in the common sample.
omnibus_group_common <- anova(m1_common_fREML, freq = FALSE)
omnibus_group_common
##
## Family: binomial
## Link function: logit
##
## Formula:
## acc ~ group + s(trial_s, k = 10) + s(trial_s, sbj, bs = "fs",
## m = 1, k = 10)
##
## Parametric Terms:
## df Chi.sq p-value
## group 3 10.52 0.0146
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(trial_s) 7.945 8.685 131.8 <2e-16
## s(trial_s,sbj) 102.756 556.000 219.1 <2e-16
###############################
# Common-sample Dunnett comparisons
###############################
# Compare each longer viewing regime with the 600-ms reference group.
dunnett_group_common <- contrast(group_emm_common, method = "trt.vs.ctrl", ref = "600")
# Report log-odds differences with multiplicity-adjusted
# confidence intervals and p values.
set.seed(1234)
summary(dunnett_group_common, infer = c(TRUE, TRUE), adjust = "mvt", df = Inf)
## contrast estimate SE df asymp.LCL asymp.UCL z.ratio p.value
## 1100 - 600 0.160 0.0991 Inf -0.0727 0.394 1.618 0.2506
## 1600 - 600 0.300 0.1000 Inf 0.0644 0.536 2.995 0.0078
## RD - 600 0.251 0.0998 Inf 0.0163 0.486 2.516 0.0322
##
## Degrees-of-freedom method: user-specified
## Results are given on the log odds ratio (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 3 estimates
## P value adjustment: mvt method for 3 tests
# Report the Dunnett-adjusted comparisons as odds ratios.
set.seed(1234)
summary(dunnett_group_common, infer = c(TRUE, TRUE), type = "response", adjust = "mvt", df = Inf)
## contrast odds.ratio SE df asymp.LCL asymp.UCL null z.ratio p.value
## 1100 / 600 1.17 0.116 Inf 0.93 1.48 1 1.618 0.2506
## 1600 / 600 1.35 0.135 Inf 1.07 1.71 1 2.995 0.0078
## RD / 600 1.29 0.128 Inf 1.02 1.63 1 2.516 0.0322
##
## Degrees-of-freedom method: user-specified
## Confidence level used: 0.95
## Conf-level adjustment: mvt method for 3 estimates
## Intervals are back-transformed from the log odds ratio scale
## P value adjustment: mvt method for 3 tests
## Tests are performed on the log odds ratio scale
Session Information
sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Europe/Athens
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] itsadug_2.5 plotfunctions_1.5 mgcv_1.9-4 nlme_3.1-169
## [5] emmeans_2.0.4 plotrix_3.8-14 sciplot_1.2-0 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.7-5 jsonlite_2.0.0 compiler_4.6.1 tidyselect_1.2.1
## [5] jquerylib_0.1.4 splines_4.6.1 yaml_2.3.12 fastmap_1.2.0
## [9] lattice_0.22-9 R6_2.6.1 generics_0.1.4 knitr_1.51
## [13] tibble_3.3.1 bslib_0.11.0 pillar_1.11.1 rlang_1.3.0
## [17] cachem_1.1.0 xfun_0.60 sass_0.4.10 estimability_2.0.0
## [21] cli_3.6.6 withr_3.0.3 magrittr_2.0.5 digest_0.6.39
## [25] grid_4.6.1 rstudioapi_0.19.0 mvtnorm_1.4-2 lifecycle_1.0.5
## [29] vctrs_0.7.3 evaluate_1.0.5 glue_1.8.1 stats4_4.6.1
## [33] rmarkdown_2.31 tools_4.6.1 pkgconfig_2.0.3 htmltools_0.5.9