-
Notifications
You must be signed in to change notification settings - Fork 1
/
model-behavioral-data-final.R
81 lines (71 loc) · 2.9 KB
/
model-behavioral-data-final.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#! /usr/bin/env Rscript
# ================================
# Script 'model-behavioral-data.R'
# ================================
# This script reads in behavioral data from two psychoacoustics experiments
# (vocoder/reverb) and models the conditions with mixed effects regression.
#
# Author: Dan McCloy <[email protected]>
# License: BSD (3-clause)
library(afex)
library(parallel)
# file paths
data_dir <- "data-behavioral"
out_dir <- "models"
log <- file.path(out_dir, "cluster-log.txt")
################################################################################
# LOAD DATA
cat("loading data\n")
load(file.path(out_dir, "dfs.RData"))
# remove dataframes with target-only definition of "truth" from workspace
rm(rev_df2, voc_df2)
################################################################################
# SET UP CLUSTER
invisible(file.remove(log))
cl <- makeForkCluster(nnodes=12, outfile=log)
################################################################################
# FIT MODELS (FAST; LIKELIHOOD RATIO TESTS FOR P-VALUES)
# define models
formulae <- list(
# reverb
rev_full_model=formula(press ~ truth*reverb*gender*attn + (1|subj)),
# vocoder
voc_full_model=formula(press ~ truth*voc_chan*gap_len*attn + (1|subj))
)
datas <- c(list(rev_df), list(voc_df))
stopifnot(length(datas) == length(formulae))
# fit models
cat("fitting models (LRT)\n")
proc.time()
lrt_models <- mapply(function(formula_, name, data_source) {
mod <- mixed(formula_, data=data_source, family=binomial(link="probit"),
method="LRT", check.contrasts=FALSE, cl=cl,
control=glmerControl(optCtrl=list(maxfun=30000)))
cat(paste(name, "finished\n"))
mod
}, formulae, names(formulae), datas)
save(lrt_models, file=file.path(out_dir, "lrt-models.RData"))
stopCluster(cl)
stop()
################################################################################
# FIT MODELS (BOOTSTRAP P-VALUES; NOT USED -- TOO SLOW)
# reverb experiment
cat("fitting models (bootstrap)\n")
proc.time()
form <- formula(press ~ truth*reverb*gender*attn + (1|subj))
rev_mod <- mixed(form, data=rev_df, family=binomial(link="probit"),
method="PB", check.contrasts=FALSE, cl=cl,
args.test=list(nsim=1000, cl=cl, seed=1234, details=2),
control=glmerControl(optCtrl=list(maxfun=30000)))
save(rev_mod, file=file.path(out_dir, "reverb-model.RData"))
# vocoder experiment
form <- formula(press ~ truth*voc_chan*gap_len*attn + (1|subj))
voc_mod <- mixed(form, data=voc_df, family=binomial(link="probit"),
method="PB", check.contrasts=FALSE, cl=cl,
args.test=list(nsim=1000, cl=cl, seed=1234, details=2),
control=glmerControl(optCtrl=list(maxfun=30000)))
save(voc_mod, file=file.path(out_dir, "vocoder-model.RData"))
proc.time()
################################################################################
# STOP CLUSTER
stopCluster(cl)