Simple detection with random sensitivity effects: Using rstanarm
Stefan M. Herzog
2016-04-13 15:28:18
1 Overview
Adapted from pp. 259–261 in:
Knoblauch, K. & Maloney, L. (2012). Modeling psychophysical data in R. New York, NY: Springer. http://www.springer.com/gp/book/9781461444749
2 Setup
suppressMessages({
library(lme4)
library(rstanarm)
library(MPDiR)
library(dplyr)
library(ggplot2); library(cowplot); library(gridExtra)
})## Warning: package 'lme4' was built under R version 3.2.3## Warning: package 'Matrix' was built under R version 3.2.2## Warning: package 'rstanarm' was built under R version 3.2.3## Warning: package 'Rcpp' was built under R version 3.2.3## Warning: package 'ggplot2' was built under R version 3.2.4## Warning: package 'cowplot' was built under R version 3.2.4## Warning: package 'gridExtra' was built under R version 3.2.43 The data
3.1 The variables
glimpse(Faces2)## Observations: 960
## Variables: 4
## $ Resp (fctr) 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0...
## $ Stim (fctr) P, P, P, P, P, P, P, P, P, P, P, P, P, P, P, A, A, A, A...
## $ Obs (fctr) S1, S1, S1, S1, S1, S1, S1, S1, S1, S1, S1, S1, S1, S1,...
## $ Image (fctr) Im1, Im2, Im3, Im4, Im5, Im6, Im7, Im8, Im9, Im10, Im11...3.2 Example rows
First two lines of data for each participant for images 1 and 16:
Faces2 %>%
filter(Image %in% c("Im1", "Im16")) %>%
group_by(Obs, Stim) %>%
slice(1:2)## Source: local data frame [64 x 4]
## Groups: Obs, Stim [64]
##
## Resp Stim Obs Image
## (fctr) (fctr) (fctr) (fctr)
## 1 0 A S1 Im16
## 2 1 P S1 Im1
## 3 0 A S2 Im16
## 4 1 P S2 Im1
## 5 1 A S3 Im16
## 6 0 P S3 Im1
## 7 0 A S4 Im16
## 8 1 P S4 Im1
## 9 1 A S5 Im16
## 10 1 P S5 Im1
## .. ... ... ... ...4 Classic approach
F2.glmm1 <- glmer(
Resp ~ Stim + (1 | Obs) + (1 | Image),
data = Faces2,
family = binomial(probit))
F2.glmm1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( probit )
## Formula: Resp ~ Stim + (1 | Obs) + (1 | Image)
## Data: Faces2
## AIC BIC logLik deviance df.resid
## 1015.8694 1035.3372 -503.9347 1007.8694 956
## Random effects:
## Groups Name Std.Dev.
## Obs (Intercept) 0.1142
## Image (Intercept) 0.7447
## Number of obs: 960, groups: Obs, 32; Image, 30
## Fixed Effects:
## (Intercept) StimP
## -0.4777 1.25735 Bayesian approach
5.1 Running the model
# simulating draws from the posterior
F2.stan_glmm1 <- stan_glmer(
Resp ~ Stim + (1 | Obs) + (1 | Image),
data = Faces2,
family = binomial(probit),
#
cores = 2,
chains = 2,
seed = 12345,
iter = 500)
# simulating draws from the prior
F2.stan_glmm1_prior <- update(F2.stan_glmm1, prior_PD = TRUE)5.2 Diagnosing MCMC
launch_shinystan(F2.stan_glmm1)5.3 Joint distributions
5.3.1 Joint prior
pairs(F2.stan_glmm1_prior,
pars = c("(Intercept)", "StimP"),
panel = points)
5.3.2 Joint posterior
pairs(F2.stan_glmm1,
pars = c("(Intercept)", "StimP"),
panel = points)
5.4 Model output
F2.stan_glmm1## stan_glmer(formula = Resp ~ Stim + (1 | Obs) + (1 | Image), data = Faces2,
## family = binomial(probit), cores = 2, chains = 2, seed = 12345,
## iter = 500)
##
## Estimates:
## Median MAD_SD
## (Intercept) -0.5 0.2
## StimP 1.2 0.3
##
## Error terms:
## Groups Name Std.Dev.
## Obs (Intercept) 0.12
## Image (Intercept) 0.80
## Num. levels: Obs 32, Image 30
##
## Sample avg. posterior predictive
## distribution of y (X = xbar):
## Median MAD_SD
## mean_PPD 0.6 0.05.5 Plot
plot(F2.stan_glmm1)
5.6 Detailed outputs
5.6.1 Fixed effects coefficients (posterior medians)
fixef(F2.stan_glmm1)## (Intercept) StimP
## -0.4844873 1.24079225.6.2 Posterior uncertainty intervals
posterior_interval(F2.stan_glmm1, pars = c("(Intercept)", "StimP", "b[(Intercept) Obs:S1]", "b[(Intercept) Obs:S2]"))## 5% 95%
## (Intercept) -0.8728055 -0.1230282
## StimP 0.8063392 1.7661444
## b[(Intercept) Obs:S1] -0.1385547 0.2645446
## b[(Intercept) Obs:S2] -0.2291267 0.20305495.6.3 Random effects coefficients (posterior medians)
coef(F2.stan_glmm1)## $Obs
## (Intercept) StimP
## S1 -0.4683157 1.240792
## S2 -0.4906644 1.240792
## S3 -0.4321791 1.240792
## S4 -0.4576067 1.240792
## S5 -0.4820218 1.240792
## S6 -0.5303548 1.240792
## S7 -0.5614861 1.240792
## S8 -0.4672616 1.240792
## S9 -0.4775527 1.240792
## S10 -0.5023560 1.240792
## S11 -0.4809788 1.240792
## S12 -0.5023078 1.240792
## S13 -0.4966944 1.240792
## S14 -0.4928748 1.240792
## S15 -0.4451241 1.240792
## S16 -0.5154148 1.240792
## S17 -0.4332706 1.240792
## S18 -0.3875348 1.240792
## S19 -0.5089230 1.240792
## S20 -0.5436357 1.240792
## S21 -0.4738265 1.240792
## S22 -0.5300315 1.240792
## S23 -0.4762652 1.240792
## S24 -0.4612836 1.240792
## S25 -0.4152815 1.240792
## S26 -0.4751397 1.240792
## S27 -0.4933243 1.240792
## S28 -0.4881294 1.240792
## S29 -0.4601846 1.240792
## S30 -0.5001663 1.240792
## S31 -0.4763255 1.240792
## S32 -0.5380504 1.240792
##
## $Image
## (Intercept) StimP
## Im1 0.01790354 1.240792
## Im2 -0.87803524 1.240792
## Im3 -0.12697195 1.240792
## Im4 -1.33315768 1.240792
## Im5 0.20137797 1.240792
## Im6 0.02844347 1.240792
## Im7 0.20946125 1.240792
## Im8 0.19063208 1.240792
## Im9 0.02454942 1.240792
## Im10 -1.34455395 1.240792
## Im11 -0.58012825 1.240792
## Im12 -1.34748340 1.240792
## Im13 -0.81107526 1.240792
## Im14 -0.25353197 1.240792
## Im15 -1.19489924 1.240792
## Im16 -0.55924805 1.240792
## Im17 -0.86180361 1.240792
## Im18 0.33138676 1.240792
## Im19 -1.95215062 1.240792
## Im20 0.49758406 1.240792
## Im21 -0.56089890 1.240792
## Im22 -0.56664624 1.240792
## Im23 0.67567659 1.240792
## Im24 -0.74923186 1.240792
## Im25 0.27523235 1.240792
## Im26 -1.58980571 1.240792
## Im27 0.76127837 1.240792
## Im28 -1.09086130 1.240792
## Im29 -0.57962939 1.240792
## Im30 -1.09612938 1.240792
##
## attr(,"class")
## [1] "coef.mer"5.7 Extracting posterior samples
as.data.frame(F2.stan_glmm1) %>% tbl_df## Source: local data frame [500 x 64]
##
## (Intercept) StimP b[(Intercept) Obs:S1] b[(Intercept) Obs:S2]
## (dbl) (dbl) (dbl) (dbl)
## 1 -0.5640800 1.3115187 -0.07676334 -0.066036530
## 2 -0.3140877 1.3033781 -0.22867572 0.085261642
## 3 -0.6344988 1.4504917 0.18635972 -0.028821125
## 4 -0.7837881 1.7517712 -0.11102367 0.016254660
## 5 -0.6141691 1.3618942 0.27793383 -0.068009979
## 6 -0.4192527 1.3099360 -0.06511322 0.006768492
## 7 -0.5202288 1.2593128 0.07777688 -0.033934369
## 8 -0.4423881 1.1970167 -0.01036270 -0.029489470
## 9 -0.2156597 0.9296101 0.04179802 -0.104895307
## 10 -0.4671435 1.3886580 -0.10852458 -0.015485983
## .. ... ... ... ...
## Variables not shown: b[(Intercept) Obs:S3] (dbl), b[(Intercept) Obs:S4]
## (dbl), b[(Intercept) Obs:S5] (dbl), b[(Intercept) Obs:S6] (dbl),
## b[(Intercept) Obs:S7] (dbl), b[(Intercept) Obs:S8] (dbl), b[(Intercept)
## Obs:S9] (dbl), b[(Intercept) Obs:S10] (dbl), b[(Intercept) Obs:S11]
## (dbl), b[(Intercept) Obs:S12] (dbl), b[(Intercept) Obs:S13] (dbl),
## b[(Intercept) Obs:S14] (dbl), b[(Intercept) Obs:S15] (dbl),
## b[(Intercept) Obs:S16] (dbl), b[(Intercept) Obs:S17] (dbl),
## b[(Intercept) Obs:S18] (dbl), b[(Intercept) Obs:S19] (dbl),
## b[(Intercept) Obs:S20] (dbl), b[(Intercept) Obs:S21] (dbl),
## b[(Intercept) Obs:S22] (dbl), b[(Intercept) Obs:S23] (dbl),
## b[(Intercept) Obs:S24] (dbl), b[(Intercept) Obs:S25] (dbl),
## b[(Intercept) Obs:S26] (dbl), b[(Intercept) Obs:S27] (dbl),
## b[(Intercept) Obs:S28] (dbl), b[(Intercept) Obs:S29] (dbl),
## b[(Intercept) Obs:S30] (dbl), b[(Intercept) Obs:S31] (dbl),
## b[(Intercept) Obs:S32] (dbl), b[(Intercept) Image:Im1] (dbl),
## b[(Intercept) Image:Im2] (dbl), b[(Intercept) Image:Im3] (dbl),
## b[(Intercept) Image:Im4] (dbl), b[(Intercept) Image:Im5] (dbl),
## b[(Intercept) Image:Im6] (dbl), b[(Intercept) Image:Im7] (dbl),
## b[(Intercept) Image:Im8] (dbl), b[(Intercept) Image:Im9] (dbl),
## b[(Intercept) Image:Im10] (dbl), b[(Intercept) Image:Im11] (dbl),
## b[(Intercept) Image:Im12] (dbl), b[(Intercept) Image:Im13] (dbl),
## b[(Intercept) Image:Im14] (dbl), b[(Intercept) Image:Im15] (dbl),
## b[(Intercept) Image:Im16] (dbl), b[(Intercept) Image:Im17] (dbl),
## b[(Intercept) Image:Im18] (dbl), b[(Intercept) Image:Im19] (dbl),
## b[(Intercept) Image:Im20] (dbl), b[(Intercept) Image:Im21] (dbl),
## b[(Intercept) Image:Im22] (dbl), b[(Intercept) Image:Im23] (dbl),
## b[(Intercept) Image:Im24] (dbl), b[(Intercept) Image:Im25] (dbl),
## b[(Intercept) Image:Im26] (dbl), b[(Intercept) Image:Im27] (dbl),
## b[(Intercept) Image:Im28] (dbl), b[(Intercept) Image:Im29] (dbl),
## b[(Intercept) Image:Im30] (dbl)# or: as.matrix(F2.stan_glmm1)5.8 Simulating posterior predictives
F2.stan_glmm1_pp <- posterior_predict(F2.stan_glmm1, seed = 98765)
dim(F2.stan_glmm1_pp)## [1] 500 960F2.stan_glmm1_pp[1:25, 1:10]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 1 1 0 1 0 1 1 0 1
## [2,] 0 1 1 0 1 1 1 1 1 0
## [3,] 1 1 1 1 1 1 1 1 1 1
## [4,] 0 0 0 1 1 1 0 1 1 0
## [5,] 1 1 1 1 1 1 1 1 1 1
## [6,] 0 0 1 0 1 1 1 1 1 1
## [7,] 1 0 1 1 1 1 1 1 1 1
## [8,] 1 1 0 1 1 1 1 1 0 1
## [9,] 1 1 1 1 1 0 1 1 1 1
## [10,] 1 1 0 1 1 1 1 1 0 0
## [11,] 1 0 1 1 1 1 1 1 1 1
## [12,] 0 1 1 0 1 1 1 1 0 0
## [13,] 1 1 1 0 1 1 1 0 1 0
## [14,] 1 1 1 1 0 1 1 1 1 1
## [15,] 1 0 1 1 1 1 1 1 1 1
## [16,] 1 1 1 1 1 1 1 1 1 1
## [17,] 1 1 1 0 1 1 1 1 1 0
## [18,] 1 0 1 0 1 1 1 1 0 1
## [19,] 1 1 1 0 1 0 1 1 1 0
## [20,] 1 0 1 1 1 1 1 1 1 1
## [21,] 1 0 1 0 1 1 1 1 1 1
## [22,] 1 1 0 1 1 1 1 1 1 1
## [23,] 1 0 1 0 1 0 1 1 0 1
## [24,] 1 0 1 0 1 1 1 1 1 0
## [25,] 1 0 1 0 1 1 1 1 1 05.9 Model comparison
5.9.1 Fit model with logit-link
F2.stan_glmm2 <- update(F2.stan_glmm1,
formula = Resp ~ Stim + (1 | Obs) + (1 | Image),
family = binomial(logit))
F2.stan_glmm2_prior <- update(F2.stan_glmm2, prior_PD = TRUE)F2.stan_glmm2## stan_glmer(formula = Resp ~ Stim + (1 | Obs) + (1 | Image), data = Faces2,
## family = binomial(logit), cores = 2, chains = 2, seed = 12345,
## iter = 500)
##
## Estimates:
## Median MAD_SD
## (Intercept) -0.8 0.4
## StimP 2.1 0.6
##
## Error terms:
## Groups Name Std.Dev.
## Obs (Intercept) 0.2
## Image (Intercept) 1.4
## Num. levels: Obs 32, Image 30
##
## Sample avg. posterior predictive
## distribution of y (X = xbar):
## Median MAD_SD
## mean_PPD 0.6 0.0fixef(F2.stan_glmm2)## (Intercept) StimP
## -0.8072523 2.0802664posterior_interval(F2.stan_glmm2, pars = c("StimP", "b[(Intercept) Obs:S1]", "b[(Intercept) Obs:S2]"))## 5% 95%
## StimP 1.2220056 2.9592993
## b[(Intercept) Obs:S1] -0.2145908 0.3703989
## b[(Intercept) Obs:S2] -0.4156460 0.40243475.9.2 Compare models
5.9.2.1 Parameter estimates
grid.arrange(plot(F2.stan_glmm1) + ggtitle("Probit model"),
plot(F2.stan_glmm2) + ggtitle("Logit model"),
ncol = 2)
5.9.2.2 Using approximate leave-one-out cross-validation (LOO)
F2.stan_glmm1_loo <- loo(F2.stan_glmm1)
F2.stan_glmm2_loo <- loo(F2.stan_glmm2)F2.stan_glmm1_loo## Computed from 500 by 960 log-likelihood matrix
##
## Estimate SE
## elpd_loo -485.4 16.4
## p_loo 32.9 1.9
## looic 970.8 32.9plot(F2.stan_glmm1_loo, label_points = TRUE)
F2.stan_glmm2_loo## Computed from 500 by 960 log-likelihood matrix
##
## Estimate SE
## elpd_loo -485.6 16.5
## p_loo 32.5 1.7
## looic 971.3 33.0plot(F2.stan_glmm2_loo, label_points = TRUE)
compare(F2.stan_glmm1_loo, F2.stan_glmm2_loo)## elpd_diff se
## -0.2 0.56 R session
sessionInfo()## R version 3.2.0 (2015-04-16)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.10.2 (Yosemite)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_0.4.3 MPDiR_0.1-16 rstanarm_2.9.0-3 Rcpp_0.12.3
## [5] lme4_1.1-11 Matrix_1.2-3
##
## loaded via a namespace (and not attached):
## [1] rstan_2.9.0-3 gtools_3.5.0 zoo_1.7-12
## [4] shinyjs_0.4.0 reshape2_1.4.1 splines_3.2.0
## [7] lattice_0.20-33 colorspace_1.2-6 htmltools_0.3
## [10] stats4_3.2.0 loo_0.1.5 yaml_2.1.13
## [13] base64enc_0.1-3 nloptr_1.0.4 DBI_0.3.1
## [16] matrixStats_0.50.1 plyr_1.8.3 stringr_1.0.0
## [19] munsell_0.4.3 gtable_0.2.0 htmlwidgets_0.6
## [22] threejs_0.2.1 codetools_0.2-14 evaluate_0.8
## [25] inline_0.3.14 knitr_1.12.3 httpuv_1.3.3
## [28] parallel_3.2.0 markdown_0.7.7 xts_0.9-7
## [31] xtable_1.8-2 scales_0.4.0 DT_0.1
## [34] shinystan_2.1.0 formatR_1.2.1 mime_0.4
## [37] gridExtra_2.2.1 ggplot2_2.1.0 digest_0.6.9
## [40] stringi_1.0-1 shiny_0.13.1 grid_3.2.0
## [43] tools_3.2.0 magrittr_1.5 shinythemes_1.0.1
## [46] MASS_7.3-45 dygraphs_0.8 assertthat_0.1
## [49] minqa_1.2.4 rmarkdown_0.9.2 R6_2.1.2
## [52] nlme_3.1-126