Statistics & Methodology
Bayesian replication test
In a very cool paper on different types of replication tests, Verhagen and Wagenmakers proposed their own variant of a Bayesian replication test (see reference below). The test compares the “proponent’s” hypothesis that the new data constitute a replication against the “skeptic’s” hypothesis that the new data constitutes a null effect. The presents simulations comparing the different Bayesian and frequentist measures, and makes a compelling argument for the proposed new measure.
For the case of a t-test, Verhagen and Wagenmakers show how one can approximate the evidential support (Bayes Factor) for the replication hypothesis over the null hypothesis, BF_r0, based on just
- the t-statistics of the original and new data
- the number of data points in the original t-test and t-test over the replication data
It occurred to me that one might be able to use the same measure to assess the replicability of individual effects in linear mixed models (i.e., the evidentiary support for the hypothesis that a specific effect in a multiple linear regression model replicates), in particular in cases where the different predictors in the model are all orthogonal to each other (e.g., in balanced designs). So one purpose of this post is to give everyone something to shout at if they disagree.
The second purpose is to post the code here. The code is really all from Josine Verhagen’s webpage. I’m posting it here for convenience. But I just
- wrote a quick and dirty wrapper function that uses the lmerTest and pbkrtest packages to provide Satterthwaite (fast and often good for large data sets) or Kenward-Roger (more precise, but a lot slower and ever more so, the larger the data) approximations of the degrees of freedom for the t-test, which are used to get the ‘number of data points’ that went into the original and replication t-test.
- update the plotting of the prior and posterior densities to use ggplot2, so that the plot can be more easily modified subsequently.
For documentation and background, please see Verhagen and Wagenmaker’s excellent paper!
whichDFMethod <- function(model) {
l = length(fitted(model))
# This threshold was set more or less arbitrarily, based on some
# initial testing. You might want to change it. Satterthwaite's
# approximation can result in hugely inflated DFs.
if (l < 3000) {
cat(paste0("Using Kenward-Roger's approximation of DFs for t-test over coefficient because the data set is small: ", l, "\n"))
return("Kenward-Roger")
} else { return("Satterthwaite") }
}
# wrapper function that calls ReplicationBayesfactorNCT with t-values and n1, n2 observations
# based on the output of mixed model analyses. The test is conducted for the k=th t-test within
# the mixed model, specified by the "coef" argument. Assumes that the models have been fit with lmerTest,
# so that the approximated DFs can be used as n1 and n2 (number of observations in the original and
# replication experiment). For that reason (and only that reason), the function requires lmerTest.
replicationBF <- function(model1, model2, yhigh = 0, limits.x = NULL, coef = 1, M = 500000) {
require(lmerTest)
require(pbkrtest)
summary1 = summary(model1, ddf=whichDFMethod(model1))
summary2 = summary(model2, ddf=whichDFMethod(model2))
print(summary1)
cat("\n\n")
print(summary2)
cat("\n\nCalculting replication Bayes Factor:\n")
out <- ReplicationBayesfactorNCT(
coef(summary1)[coef,"t value"],
coef(summary2)[coef,"t value"],
coef(summary1)[coef,"df"] + 1, # Since the ReplicationBayesfactorNCT function subtracts 1 for the DFs (for sample == 1)
coef(summary2)[coef,"df"] + 1, # Since the ReplicationBayesfactorNCT function subtracts 1 for the DFs (for sample == 1)
sample = 1,
plot = 1,
post = 1,
yhigh = yhigh,
limits.x = limits.x,
M = M
)
return (out)
}
# Replication Bayes Factor taken from Josine Verhagen's webpage
# (http://josineverhagen.com/?page_id=76).
# See Verhagen and Wagenmakers (2014) for details.
#
# Minor modifications to the plotting made by fjaeger@ur.rochester.edu
ReplicationBayesfactorNCT <- function(
tobs, # t value in first experiment
trep, # t value in replicated experiment
n1, # first experiment: n in group 1 or total n
n2, # second experiment: n in group 1 or total n
m1 = 1, # first experiment: n in group 2 or total n
m2 = 1, # second experiment: n in group 2 or total n
sample = 1, # 1 = one sample t-test (or within), 2 = two sample t-test
wod = dir, # working directory
plot = 0, # 0 = no plot 1 = replication
post = 0, # 0 = no posterior, 1 = estimate posterior
M = 500000,
yhigh = 0,
limits.x = NULL
)
{
require(MCMCpack)
require(ggplot2) # added by TFJ
##################################################################
#STEP 1: compute the prior for delta based on the first experiment
################################### ###############################
D <- tobs
if (sample==1) {
sqrt.n.orig <- sqrt(n1)
df.orig <- n1 -1
}
if (sample==2) {
sqrt.n.orig <- sqrt( 1/(1/n1+1/m1)) # two sample alternative for sqrt(n)
df.orig <- n1 + m1 -2 #degrees of freedom
}
# To find out quickly the lower and upper bound areas, the .025 and .975 quantiles of tobs at a range of values for D are computed.
# To make the algorithm faster, only values in a reasonable area around D are computed.
# Determination of area, larger with large D and small N
range.D <- 4 + abs(D/(2*sqrt.n.orig))
sequence.D <- seq(D,D + range.D,.01) #make sequence with range
# determine which D gives a quantile closest to tobs with an accuracy of .01
options(warn=-1)
approximatelow.D <- sequence.D[which( abs(qt(.025,df.orig,sequence.D)-tobs)==min(abs(qt(.025,df.orig,sequence.D)-tobs)) )]
options(warn=0)
# Then a more accurate interval is computed within this area
# Make sequence within .01 from value found before
sequenceappr.D <- seq((approximatelow.D-.01),(approximatelow.D+.01),.00001)
# determine which D gives a quantile closest to tobs with an accuracy of .00001
low.D <- sequenceappr.D[which( abs(qt(.025,df.orig,sequenceappr.D)-tobs)==min(abs(qt(.025,df.orig,sequenceappr.D)-tobs)) )]
# Compute standard deviation for the corresponding normal distribution.
sdlow.D <- (D-low.D)/qnorm(.025)
# compute prior mean and as for delta
prior.mudelta <- D/sqrt.n.orig
prior.sdelta <- sdlow.D/sqrt.n.orig
##################################################################
#STEP 2: Compute Replication Bayes Factor
##################################################################
# For one sample t-test: within (group2==1) or between (group2==vector)
if (sample==1)
{
df.rep <- n2-1
sqrt.n.rep <- sqrt(n2)
Likelihood.Y.H0 <- dt(trep,df.rep)
sample.prior <- rnorm(M,prior.mudelta,prior.sdelta)
options(warn=-1)
average.Likelihood.H1.delta <- mean(dt(trep,df.rep,sample.prior*sqrt.n.rep))
options(warn=0)
BF <- average.Likelihood.H1.delta/Likelihood.Y.H0
}
# For two sample t-test:
if (sample==2) {
df.rep <- n2 + m2 -2
sqrt.n.rep <- sqrt(1/(1/n2+1/m2))
Likelihood.Y.H0 <- dt(trep,df.rep)
sample.prior <- rnorm(M,prior.mudelta,prior.sdelta)
options(warn=-1)
average.Likelihood.H1.delta <- mean(dt(trep,df.rep,sample.prior*sqrt.n.rep))
options(warn=0)
BF <- average.Likelihood.H1.delta/Likelihood.Y.H0
}
#################################################################
#STEP 3: Posterior distribution
#################################################################
if( post == 1) {
options(warn=-1)
likelihood <- dt(trep,df.rep,sample.prior*sqrt.n.rep)
prior.density <- dnorm(sample.prior,prior.mudelta,prior.sdelta)
likelihood.x.prior <- likelihood * prior.density
LikelihoodXPrior <- function(x) {dnorm(x,prior.mudelta,prior.sdelta) * dt(trep,df.rep,x*sqrt.n.rep) }
fact <- integrate(LikelihoodXPrior,-Inf,Inf)
posterior.density <- likelihood.x.prior/fact$value
PosteriorDensityFunction <- function(x) {(dnorm(x,prior.mudelta,prior.sdelta) * dt(trep,df.rep,x*sqrt.n.rep))/fact$value }
options(warn=0)
mean <- prior.mudelta
sdh <- prior.mudelta + .5*(max(sample.prior)[1] - prior.mudelta)
sdl <- prior.mudelta - .5*(prior.mudelta - min(sample.prior)[1])
dev <- 2
for ( j in 1:10) {
rangem <- seq((mean-dev)[1],(mean+dev)[1],dev/10)
rangesdh <- seq((sdh-dev)[1],(sdh+dev)[1],dev/10)
rangesdl <- seq((sdl-dev)[1],(sdl+dev)[1],dev/10)
perc <- matrix(0,length(rangem),3)
I<-min(length(rangem), length(rangesdh),length(rangesdl) )
for ( i in 1:I) {
options(warn=-1)
vpercm <- integrate(PosteriorDensityFunction, -Inf, rangem[i])
perc[i,1]<- vpercm$value
vpercsh <- integrate(PosteriorDensityFunction, -Inf, rangesdh[i])
perc[i,2]<- vpercsh$value
vpercsl <- integrate(PosteriorDensityFunction, -Inf, rangesdl[i])
perc[i,3]<- vpercsl$value
options(warn=0)
}
mean <- rangem[which(abs(perc[,1]-.5)== min(abs(perc[,1]-.5)))]
sdh <- rangesdh[which(abs(perc[,2]-pnorm(1))== min(abs(perc[,2]-pnorm(1))))]
sdl <- rangesdl[which(abs(perc[,3]-pnorm(-1))== min(abs(perc[,3]-pnorm(-1))))]
dev <- dev/10
}
posterior.mean <- mean
posterior.sd <- mean(c(abs(sdh- mean),abs(sdl - mean)))
}
if (post != 1) {
posterior.mean <- 0
posterior.sd <- 0
}
###########OUT
dat.SD=new.env()
dat.SD$BF <- BF
dat.SD$prior.mean= round(prior.mudelta,2)
dat.SD$prior.sd= round(prior.sdelta,2)
dat.SD$post.mean= round(posterior.mean,2)
dat.SD$post.sd= round(posterior.sd,2)
dat.SD=as.list(dat.SD)
###########################################
#PLOT
###########################################
if (plot==1)
{
if (post != 1) {
options(warn =-1)
rp <- ReplicationPosterior(trep,prior.mudelta,prior.sdelta,n2,m2=1,sample=sample)
options(warn =0)
posterior.mean <- rp[[1]]
posterior.sd <- rp[[2]]
}
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
font.lab = 2, cex.axis = 1.3, bty = "n", las=1)
high <- dnorm(posterior.mean,posterior.mean,posterior.sd)
add <- high/5
if(yhigh==0) {
yhigh <- high + add
}
scale <- 3
if (is.null(limits.x)) {
min.x <- min((posterior.mean - scale*posterior.sd),(prior.mudelta - scale*prior.sdelta))
max.x <- max((posterior.mean + scale*posterior.sd),(prior.mudelta + scale*prior.sdelta))
} else {
min.x <- limits.x[1]
max.x <- limits.x[2]
}
ggplot(
data = data.frame(x = 0),
mapping = aes(
x = x
)
) +
scale_x_continuous(expression(Effect ~ size ~ delta),
limits = c(min.x, max.x)
) +
scale_y_continuous("Density",
limits = c(0, yhigh)
) +
stat_function(
fun = function(x) dnorm(x, prior.mudelta, prior.sdelta),
aes(linetype = "Before replication"),
size = .75
) +
stat_function(
fun = function(x) dnorm(x, posterior.mean, posterior.sd),
aes(linetype = "After replication"),
size = 1
) +
scale_linetype_manual("",
values = c("Before replication" = 2,
"After replication" = 1)
) +
geom_point(
x = 0, y = dnorm(0,prior.mudelta,prior.sdelta),
color = "gray", size = 3
) +
geom_point(
x = 0, y = dnorm(0,posterior.mean,posterior.sd),
color = "gray", size = 3
) +
geom_segment(
x = 0, xend = 0,
y = dnorm(0,prior.mudelta,prior.sdelta), yend= dnorm(0,posterior.mean,posterior.sd),
color = "gray"
) +
annotate(
geom = "text",
x = sum(limits.x) / 2, y = yhigh - 1/10 * yhigh,
label = as.expression(bquote(BF[r0] == .(round(BF, digits=2))))
) +
theme(
legend.position = "bottom",
legend.key.width = unit(1.5, "cm")
)
}
return(dat.SD)
}
Verhagen, J., & Wagenmakers, E. J. (2014). Bayesian tests to quantify the result of a replication attempt. Journal of Experimental Psychology: General, 143(4), 1457.
Some thoughts on Healey et al (2014) failure to find syntactic priming in conversational speech
In a recent PLoS one article, Healey, Purver, and Howes (2014) investigate syntactic priming in conversational speech, both within speakers and across speakers. Healey and colleagues follow Reitter et al (2006) in taking a broad-coverage approach to the corpus-based study of priming. Rather than to focus on one or a few specific structures, Healey and colleagues assess lexical and structural similarity within and across speakers. The paper concludes with the interesting claim that there is no evidence for syntactic priming within speaker and that alignment across speakers is actually less than expected by chance once lexical overlap is controlled for. Given more than 30 years of research on syntactic priming, this is a rather interesting claim. As some folks have Twitter-bugged me (much appreciated!), I wanted to summarize some quick thoughts here. Apologies in advance for the somewhat HLP-lab centric view. If you know of additional studies that seem relevant, please join the discussion and post. Of course, Healey and colleagues are more than welcome to respond and correct me, too.
First, the claim by Healey and colleagues that “previous work has not tested for general syntactic repetition effects in ordinary conversation independently of lexical repetition” (Healey et al 2014, abstract) isn’t quite accurate.
more on old and new lme4
(This is another guest post by Klinton Bicknell.)
This is an update to my previous blog post, in which I observed that post-version-1.0 versions of the lme4 package yielded worse model fits than old pre-version-1.0 versions for typical psycholinguistic datasets, and I gave instructions for installing the legacy lme4.0 package. As I mentioned there, however, lme4 is under active development, the short version of this update post is to say that it seems that the latest versions of the post-version-1.0 lme4 now yield models that are just as good, and often better than lme4.0! This seems to be due to the use of a new optimizer, better convergence checking, and probably other things too. Thus, installing lme4.0 now seems only useful in special situations involving old code that expects the internals of the models to look a certain way. Life is once again easier thanks to the furious work of the lme4 development team!
[update: Since lme4 1.1-7 binaries are now on CRAN, this paragraph is obsolete.] One minor (short-lived) snag is that the current version of lme4 on CRAN (1.1-6) is overzealous in displaying convergence warnings, and displays them inappropriately in many cases where models have in fact converged properly. This will be fixed in 1.1-7 (more info here). To avoid them for now, the easiest thing to do is probably to install the current development version of lme4 1.1-7 from github like so:
library("devtools"); install_github("lme4/lme4")
Read on if you want to hear more details about my comparisons of the versions.
The ‘softer kind’ of tutorial on linear mixed effect regression
I recently was pointed to this nice and very accessible tutorial on linear mixed effects regression and how to run them in R by Bodo Winter (at UC Merced). If you don’t have much or any background in this type of model, I recommend you pair it with a good conceptual introduction to these models like Gelman and Hill 2007 and perhaps some slides from our LSA 2013 tutorial.
There are a few thing I’d like to add to Bodo’s suggestions regarding how to report your results:
- be clear how you coded the variables since this does change the interpretation of the coefficients (the betas that are often reported). E.g. say whether you sum- or treatment-coded your factors, whether you centered or standardized continuous predictors etc. As part of this, also be clear about the direction of the coding. For example, state that you “sum-coded gender as female (1) vs. male (-1)”. Alternatively, report your results in a way that clearly states the directionality (e.g., “Gender=male, beta = XXX”).
- please also report whether collinearity was an issue. E.g., report the highest fixed effect correlations.
Happy reading.
old and new lme4
(This is a guest post by Klinton Bicknell.)
update 2014-06-24: Using lme4.0 probably isn’t necessary anymore. See post here.
The lme4 package‘s major 1.0 release was back in August. I and others have noticed that for typical psycholinguistic datasets, the new >=1.0 versions of lme4 often yield models with substantially poorer fits to the data than the old pre-1.0 versions (sometimes worse by many points of log likelihood), which suggests that the new lme4 isn’t as reliably converging to the actual maximum likelihood (or REML) solution. Since unconverged models yield misleading inferences about model parameters, it’s useful to be able to fit models using the old pre-1.0 lme4.
Happily, the lme4 developers have created a new package (named “lme4.0”), which is a bugfix-only version of the old pre-1.0 lme4. This allows for the installation of both old and new versions of lme4 side-by-side. As of this posting, lme4.0 is not yet on CRAN, but is installable by performing the following steps: Read the rest of this entry »
Another example of recording spoken productions over the web
A few days ago, I posted a summary of some recent work on syntactic alignment with Kodi Weatherholtz and Kathryn Campell-Kibler (both at The Ohio State University), in which we used the WAMI interface to collect speech data for research on language production over Amazon’s Mechanical Turk.
Socially-mediated syntactic alignment
The first step in our OSU-Rochester collaboration on socially-mediated syntactic alignment has been submitted a couple of weeks ago. Kodi Weatherholtz in Linguistics at The Ohio State University took the lead in this project together with Kathryn Campbell-Kibler (same department) and me.
We collected spoken picture descriptions via Amazon’s crowdsourcing platform Mechanical Turk to investigate how social attitude towards an interlocutor and conflict management styles affected syntactic priming. Our paradigm combines Read the rest of this entry »
Ways of plotting map data in R (and python)
Thanks to Scott Jackson, Daniel Ezra Johnson, David Morris, Michael Shvartzman, and Nathanial Smith for the recommendations and pointers to the packages mentioned below.
- R:
- The maps, mapsextra, and maptools packages provide data and tools to plot world, US, and a variety of regional maps (see also mapproj and mapdata). This, combined with ggplot2 is also what we used in Jaeger et al., (2011, 2012) to plot distributions over world maps. Here’s an example from ggplot2 with maps.
Example use of ggplot2 combined with the maps package (similar to the graphs created for Jaeger et al., 2011, 2012).
And a belated welcome to Scott Grimm
Scott Grimm just joined our faculty in Linguistics at Rochester last month. So today he got his belated Rochester welcome:
Scott joins the Center of Language Sciences at Rochester with an unnecessary number of degrees. Read the rest of this entry »
Updated slides on GLM, GLMM, plyr, etc. available
Some of you asked for the slides to the Mixed effect regression class I taught at the 2013 LSA Summer Institute in Ann Arbor, MI. The class covered some Generalized Linear Model, Generalized Linear Mixed Models, extensions beyond the linear model, simulation-based approaches to assessing the validity (or power) of your analysis, data summarization and visualization, and reporting of results. The class included slides from Maureen Gillespie, Dave Kleinschmidt, and Judith Degen (see above link). Dave even came by to Ann Arbor and gave his lecture on the awesome power of plyr (and reshape etc.), which I recommend. You might also just browse through them to get an idea of some new libraries (such as Stargazer for quick and nice looking latex tables). There’s also a small example to work through for time series analysis (for beginners).
Almost all slides were created in knitr and latex (very conveniently integrated into RStudio — I know some purists hate it, but comm’on), so that the code on the slides is the code that generated the output on the slides. Feedback welcome.
Workshop announcement (Tuebingen): Advances in Visual Methods for Linguistics
This workshop on data visualization might be of interest to a lot of you. I wish I could just hop over the pond.
- Date: 24-Sept-2014 – 26-Sept-2014
- Location: Tuebingen, Germany
- Contact Person: Fabian Tomaschek (contact@avml-meeting.com)
- Web Site: http://avml-meeting.com
- Call Deadlines: 21 March / 18 April
The AVML-meeting offers a meeting place for all linguists from all fields who are interested in elaborating and improving their data visualization skills and methods. The meeting consists of a one-day hands-on workshop Read the rest of this entry »
Is my analysis problematic? A simulation-based example
This post is in reply to a recent question on in ling-R-lang by Meredith Tamminga. Meredith was wondering whether an analysis she had in mind for her project was circular, causing the pattern of results predicted by the hypothesis that she was interested in testing. I felt her question (described below in more detail) was an interesting example that might best be answered with some simulations. Reasoning through an analysis can, of course, help a lot in understanding (or better, as in Meredith’s case, anticipating) problems with the interpretation of the results. Not all too infrequently, however, I find that intuition fails or isn’t sufficiently conclusive. In those cases, simulations can be a powerful tool in understanding your analysis. So, I decided to give it a go and use this as an example of how one might approach this type of question.
Going full Bayesian with mixed effects regression models
Thanks to some recently developed tools, it’s becoming very convenient to do full Bayesian inference for generalized linear mixed-effects models. First, Andrew Gelman et al. have developed Stan, a general-purpose sampler (like BUGS/JAGS) with a nice R interface which samples from models with correlated parameters much more efficiently than BUGS/JAGS. Second, Richard McElreath has written glmer2stan, an R package that essentially provides a drop-in replacement for the lmer command that runs Stan on a generalized linear mixed-effects model specified with a lme4-style model formula.
This means that, in many cases, you simply simply replace calls to (g)lmer() with calls to glmer2stan():
library(glmer2stan)
library(lme4)
lmer.fit <- glmer(accuracy ~ (1|item) + (1+condition|subject) + condition,
data=data, family='binomial')
summary(lmer.fit)
library(glmer2stan)
library(rstan)
stan.fit <- glmer2stan(accuracy ~ (1|item) + (1+condition|subject) + condition,
data=data, family='binomial')
stanmer(stan.fit)
There’s the added benefit that you get a sample from the full, joint posterior distribution of the model parameters
Read on for more about the advantage of this approach and how to use it.
Running phonetic (adaptation) experiments online
I’ve developed some JavaScript code that somewhat simplifies running experiments online (over, e.g., Amazon’s Mechanical Turk). There’s a working demo, and you can download or fork the source code to tinker with yourself. The code for the core functionality which controls stimulus display, response collection, etc. is also available in its own repository if you just want to build around that.
If you notice a bug, or have a feature request, open an issue on the issue tracker (preferred), or comment here with questions and ideas. And, of course, if you want to contribute, please go ahead and submit a pull request. Everything’s written in HTML, CSS, and JavaScript (+JQuery) and aims to be as extensible as possible. Happy hacking!
If you find this code useful for your purposes, please refer others to this page. If you’d like to cite something to acknowledge this code or your own code based on this code, the following is the paper in which we first used this paradigm:
- Kleinschmidt, D. F., and Jaeger, T. F. 2012. A continuum of phonetic adaptation: Evaluating an incremental belief-updating model of recalibration and selective adaptation. Proceedings of the 34th Annual Meeting of the Cognitive Science Society (CogSci12), 605-610. Austin, TX: Cognitive Science Society.
A more detailed journal paper is currently under review. If you’re interested, subscribe to this post and get the update when we post the paper here once it’s out (or contact me if you can’t wait).
HLP/Jaeger lab blog 