Author: tiflo

Bayesian replication test

Posted on Updated on

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.

Advertisements

Two post-doc positions available in our lab

Posted on Updated on

We are searching for two outstanding post-doctoral researchers to join our lab. Both positions offer competitive NIH-level post-doctoral salaries for up to 3 years, an annual travel budget, and moving expenses. The lab has a good record at job placement, with three of the four most recent post-docs now holding tenure-track positions in linguistics or psychology.
  • Project 1: Inference and learning during speech perception and adaptation
  • Project 2: Web-based self-administered speech therapy

Although we mention preferred specializations below, applicants from any fields in the cognitive and language sciences are welcome. While candidates will join an active project, candidates are welcome/encouraged to also develop their own independent research program.  In case of doubt, please contact Florian Jaeger at fjaeger@bcs.rochester.edu, rather than to self-select not to apply.

Interested candidates should contact HLP lab manager Olga Nikolayeva (onikolay@u.rochester.edu) along with:

Read the rest of this entry »

What did you read in 2015?

Posted on Updated on

Another year has passed and academic platform bombard us with end-of-year summaries. So, here are the most-read HLP Lab papers of 2015. Congratulations to Dave Kleinschmidt, who according to ResearchGate leads the 2015 HLP Lab pack with his beautiful paper on the ideal adapter framework for speech perception, adaptation, and generalization. The paper was cited 22 times in the first 6 months of being published! Well deserved, I think … as a completely neutral (and non-ideal) observer ;).

Screen Shot 2016-01-20 at 10.07.04 AM
Most read HLP Lab papers on ResearchGate. Speech perception, syntactic alignment in production, and … typology!

Academia.edu mostly agreed, Read the rest of this entry »

The (in)dependence of pronunciation variation on the time course of lexical planning

Posted on Updated on

Language, Cognition, and Neuroscience just published Esteban Buz​’s paper on the relation between the time course of lexical planning  and the detail of articulation (as hypothesized by production ease accounts).

Several recent proposals hold that much if not all of explainable pronunciation variation (variation in the realization of a word) can be reduced to effects on the ease of lexical planning. Such production ease accounts have been proposed, for example, for effects of frequency, predictability, givenness, or phonological overlap to recently produced words on the articulation of a word. According to these account, these effects on articulation are mediated through parallel effects on the time course of lexical planning (e.g., recent research by Jennifer Arnold, Jason Kahn, Duane Watson, and others; see references in paper).

 

This would indeed offer a parsimonious explanation of pronunciation variation. However, the critical test for this claim is a mediation analysis, Read the rest of this entry »

Special issue in “Cross-linguistic Psycholinguistics”

Posted on Updated on

At long last! It’s my great pleasure to announce the publication of the special issue on “Laboratory in the field: advances in cross-linguistic psycholinguistics”, edited by Alice Harris (UMass), Elisabeth Norcliffe (MPI, Nijmegen), and me (Rochester), in Language, Cognition, and Neuroscience. It is an exciting collection of cross-linguistic studies on language production and comprehension and it feels great to see the proofs for the whole shiny issue:

Front cover of special issue
Front cover of special issue

Read the rest of this entry »

CUNY 2015 plenary

Posted on Updated on

As requested by some, here are the slides from my 2015 CUNY Sentence Processing Conference plenary last week:

This slideshow requires JavaScript.

I’m posting them here for discussion purposes only. During the Q&A several interesting points were raised. For example Read the rest of this entry »

HLP Lab at CUNY 2015

Posted on Updated on

We hope to see y’all at CUNY in a few weeks. In the interest of hopefully luring to some of our posters, here’s an overview of the work we’ll be presenting. In particular, we invite our reviewers, who so boldly claimed (but did not provide references for the) triviality of our work ;), to visit our posters and help us mere mortals understand.

  • Articulation and hyper-articulation
  • Unsupervised and supervised learning during speech perception
  • Syntactic priming and implicit learning during sentence comprehension
  • Uncovering the biases underlying language production through artificial language learning

Interested in more details? Read on. And, as always, I welcome feedback. (to prevent spam, first time posters are moderated; after that your posts will always directly show)

Read the rest of this entry »