Data analyses for manuscript: Learners constrain their linguistic generalizations using preemption but not entrenchment: Evidence from artificial language learning and computational modeling with adults and children

Load packages and helper functions

Packages

library(psych)
library(plyr)
library(doBy)
library(cowplot)
library(reshape2)
library(lme4)
library(brms)
library(tidyr)
library(tidyverse)
library(data.table)
library(janitor)
library(brms)
library(yarrr)
library(knitr)


source("lmedrop.R")
source("myCenter.R")
source("lizCenter.R")
source("summarySEwithin.R")
source("summarySE.R")
source("normDataWithin.R")
source("BF.R")
source("Bf_range.R")
source("Bf_powercalc.R")


theme_set(theme_bw())

Helper functions

SummarySE

This function can be found on the website “Cookbook for R”.

http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#Helper functions

It summarizes data, giving count, mean, standard deviation, standard error of the mean, and confidence intervals (default 95%).

data: a data frame.

measurevar: the name of a column that contains the variable to be summariezed

groupvars: a vector containing names of columns that contain grouping variables

na.rm: a boolean that indicates whether to ignore NA’s

conf.interval: the percent range of the confidence interval (default is 95%)

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
    require(plyr)

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

SummarySEwithin

This function can be found on the website “Cookbook for R”.

http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#Helper functions

It summarizes data, handling within-subjects variables by removing inter-subject variability. It will still work if there are no within-S variables. It gives count, un-normed mean, normed mean (with same between-group mean), standard deviation, standard error of the mean, and confidence intervals. If there are within-subject variables, calculate adjusted values using method from Morey (2008).

data: a data frame.

measurevar: the name of a column that contains the variable to be summarized

betweenvars: a vector containing names of columns that are between-subjects variables

withinvars: a vector containing names of columns that are within-subjects variables

idvar: the name of a column that identifies each subject (or matched subjects)

na.rm: a boolean that indicates whether to ignore NA’s

conf.interval: the percent range of the confidence interval (default is 95%)

summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
                            idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {

  # Ensure that the betweenvars and withinvars are factors
  factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
    FUN=is.factor, FUN.VALUE=logical(1))

  if (!all(factorvars)) {
    nonfactorvars <- names(factorvars)[!factorvars]
    message("Automatically converting the following non-factors to factors: ",
            paste(nonfactorvars, collapse = ", "))
    data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
  }

  # Get the means from the un-normed data
  datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
                     na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)

  # Drop all the unused columns (these will be calculated with normed data)
  datac$sd <- NULL
  datac$se <- NULL
  datac$ci <- NULL

  # Norm each subject's data
  ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)

  # This is the name of the new column
  measurevar_n <- paste(measurevar, "_norm", sep="")

  # Collapse the normed data - now we can treat between and within vars the same
  ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
                      na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)

  # Apply correction from Morey (2008) to the standard error and confidence interval
  #  Get the product of the number of conditions of within-S variables
  nWithinGroups    <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
                           FUN.VALUE=numeric(1)))
  correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )

  # Apply the correction factor
  ndatac$sd <- ndatac$sd * correctionFactor
  ndatac$se <- ndatac$se * correctionFactor
  ndatac$ci <- ndatac$ci * correctionFactor

  # Combine the un-normed means with the normed results
  merge(datac, ndatac)
}

normDataWithin

This function is used by the SummarySEWithin fucntion above. It can be found on the website “Cookbook for R”.

http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#Helper functions

From that website: Norms the data within specified groups in a data frame; it normalizes each subject (identified by idvar) so that they have the same mean, within each group specified by betweenvars.

data: a data frame

idvar: the name of a column that identifies each subject (or matched subjects)

measurevar: the name of a column that contains the variable to be summarized

betweenvars: a vector containing names of columns that are between-subjects variables

na.rm: a boolean that indicates whether to ignore NA’s

normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
                           na.rm=FALSE, .drop=TRUE) {
    require(plyr)

    # Measure var on left, idvar + between vars on right of formula.
    data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
     .fun = function(xx, col, na.rm) {
        c(subjMean = mean(xx[,col], na.rm=na.rm))
      },
      measurevar,
      na.rm
    )

    # Put the subject means with original data
    data <- merge(data, data.subjMean)

    # Get the normalized data in a new column
    measureNormedVar <- paste(measurevar, "_norm", sep="")
    data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
                               mean(data[,measurevar], na.rm=na.rm)

    # Remove this subject mean column
    data$subjMean <- NULL

    return(data)
}

myCenter

This function outputs the centered values of an variable, which can be a numeric variable, a factor, or a data frame. It was taken from Florian Jaegers blog https://hlplab.wordpress.com/2009/04/27/centering-several-variables/.

From his blog:

-If the input is a numeric variable, the output is the centered variable.

-If the input is a factor, the output is a numeric variable with centered factor level values. That is, the factor’s levels are converted into numerical values in their inherent order (if not specified otherwise, R defaults to alphanumerical order). More specifically, this centers any binary factor so that the value below 0 will be the 1st level of the original factor, and the value above 0 will be the 2nd level.

-If the input is a data frame or matrix, the output is a new matrix of the same dimension and with the centered values and column names that correspond to the colnames() of the input preceded by “c” (e.g. “Variable1” will be “cVariable1”).

myCenter= function(x) {
  if (is.numeric(x)) { return(x - mean(x, na.rm=T)) }
    if (is.factor(x)) {
        x= as.numeric(x)
        return(x - mean(x, na.rm=T))
    }
    if (is.data.frame(x) || is.matrix(x)) {
        m= matrix(nrow=nrow(x), ncol=ncol(x))
        colnames(m)= paste("c", colnames(x), sep="")
    
        for (i in 1:ncol(x)) {
        
            m[,i]= myCenter(x[,i])
        }
        return(as.data.frame(m))
    }
}

lizCenter

This function provides a wrapper around myCenter allowing you to center a specific list of variables from a dataframe. The input is a dataframe (x) and a list of the names of the variables which you wish to center (listfname). The output is a copy of the dataframe with a column (numeric) added for each of the centered variables with each one labelled with it’s previous name with “.ct” appended. For example, if x is a dataframe with columns “a” and “b” lizCenter(x, list(“a”, “b”)) will return a dataframe with two additional columns, a.ct and b.ct, which are numeric, centered codings of the corresponding variables.

lizCenter= function(x, listfname) 
{
    for (i in 1:length(listfname)) 
    {
        fname = as.character(listfname[i])
        x[paste(fname,".ct", sep="")] = myCenter(x[fname])
    }
        
    return(x)
}

###get_coeffs This function allows us to inspect particular coefficients from the output of an lme model by putting them in table.

x: the output returned when running lmer or glmer (i.e. an object of type lmerMod or glmerMod)

list: a list of the names of the coefficients to be extracted (e.g. c(“variable1”, “variable1:variable2”))

get_coeffs <- function(x,list){(as.data.frame(summary(x)$coefficients)[list,])}

Bf

This function is equivalent to the Dienes (2008) calculator which can be found here: http://www.lifesci.sussex.ac.uk/home/Zoltan_Dienes/inference/Bayes.htm.

The code was provided by Baguely and Kayne (2010) and can be found here: http://www.academia.edu/427288/Review_of_Understanding_psychology_as_a_science_An_introduction_to_scientific_and_statistical_inference

Bf<-function(sd, obtained, uniform, lower=0, upper=1, meanoftheory=0,sdtheory=1, tail=1){
 area <- 0
 if(identical(uniform, 1)){
  theta <- lower
  range <- upper - lower
  incr <- range / 2000
  for (A in -1000:1000){
     theta <- theta + incr
     dist_theta <- 1 / range
     height <- dist_theta * dnorm(obtained, theta, sd)
     area <- area + height * incr
  }
 }else
   {theta <- meanoftheory - 5 * sdtheory
    incr <- sdtheory / 200
    for (A in -1000:1000){
      theta <- theta + incr
      dist_theta <- dnorm(theta, meanoftheory, sdtheory)
      if(identical(tail, 1)){
        if (theta <= 0){
          dist_theta <- 0
        } else {
          dist_theta <- dist_theta * 2
        }
      }
      height <- dist_theta * dnorm(obtained, theta, sd)
      area <- area + height * incr
    }
 }
 LikelihoodTheory <- area
 Likelihoodnull <- dnorm(obtained, 0, sd)
 BayesFactor <- LikelihoodTheory / Likelihoodnull
 ret <- list("LikelihoodTheory" = LikelihoodTheory,"Likelihoodnull" = Likelihoodnull, "BayesFactor" = BayesFactor)
 ret
} 

Bf power calculation

This works with the Bf function above. It requires the same values as that function (i.e. the obtained mean and SE for the current sample, a value for the predicted mean, which is set to be sdtheory (with meanoftheory=0), and the current number of participants N). However, rather than returning a BF for the current sample, it works out what the BF would be for a range of different subject numbers (assuming that the SE scales with sqrt(N)).

Bf_powercalc<-function(sd, obtained, uniform, lower=0, upper=1, meanoftheory=0, sdtheory=1, tail=2, N, min, max)
{
  
  x = c(0)
  y = c(0)

# note: working out what the difference between N and df is (for the contrast between two groups, this is 2; for constraints where there is 4 groups this will be 3, etc.)

  for(newN in min : max)
  {
    B = as.numeric(Bf(sd = sd*sqrt(N/newN), obtained, uniform, lower, upper, meanoftheory, sdtheory, tail)[3])
    x= append(x,newN) 
    y= append(y,B)
    output = cbind(x,y)
    
  } 
  output = output[-1,] 
  return(output) 
}

Bf range

This works with the Bf function above. It requires the obtained mean and SE for the current sample and works out what the BF would be for a range of predicted means (which are set to be sdtheoryrange (with meanoftheory=0)).

Bf_range<-function(sd, obtained, meanoftheory=0, sdtheoryrange, tail=1)
{
  
  x = c(0)
  y = c(0)
  
  for(sdi in sdtheoryrange)
  {
    B = as.numeric(Bf(sd, obtained, meanoftheory=0, uniform = 0, sdtheory=sdi, tail)[3])
    
    x= append(x,sdi)  
    y= append(y,B)
    output = cbind(x,y)
    
  } 
  output = output[-1,] 
  colnames(output) = c("sdtheory", "BF")
  return(output) 
}

Experiment 1

Load data

# Experiment 1 - verb study with adults

#Create the dataframes that we will be working on
combined_production_data.df <- read.csv("exp1_production_data.csv")

combined_judgment_data.df <- read.csv("exp1_judgment_data.csv")
combined_judgment_data.df$restricted_verb <- factor(combined_judgment_data.df$restricted_verb)
combined_judgment_data.df$condition <- factor(combined_judgment_data.df$condition)


#separately for entrenchment and preemption

#entrenchment
entrenchment_production.df <- subset(combined_production_data.df, condition == "entrenchment")
entrenchment_production.df$semantically_correct <- as.numeric(entrenchment_production.df$semantically_correct)
entrenchment_production.df$transitivity_test_scene2 <- factor(entrenchment_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in entrenchment
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

entrenchment_production.df$det1 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction1", 1, 0)
entrenchment_production.df$det2 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction2", 1, 0)
entrenchment_production.df$other <- ifelse(entrenchment_production.df$det_lenient_adapted == "other", 1, 0)
entrenchment_production.df$none <- ifelse(entrenchment_production.df$det_lenient_adapted == "none", 1, 0)


entrenchment_judgment.df <- subset(combined_judgment_data.df, condition == "entrenchment")
entrenchment_judgment.df$semantically_correct <- factor(entrenchment_judgment.df$semantically_correct)
entrenchment_judgment.df$transitivity_test_scene2 <- factor(entrenchment_judgment.df$transitivity_test_scene2)
entrenchment_judgment.df$restricted_verb <- factor(entrenchment_judgment.df$restricted_verb)


#preemption
preemption_production.df <- subset(combined_production_data.df, condition == "preemption")
preemption_production.df$semantically_correct <- as.numeric(preemption_production.df$semantically_correct)  #actually, all NAs here
preemption_production.df$transitivity_test_scene2 <- factor(preemption_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in pre-emption
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

preemption_production.df$det1 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction1", 1, 0)
preemption_production.df$det2 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction2", 1, 0)
preemption_production.df$other <- ifelse(preemption_production.df$det_lenient_adapted == "other", 1, 0)
preemption_production.df$none <- ifelse(preemption_production.df$det_lenient_adapted == "none", 1, 0)


preemption_judgment.df <- subset(combined_judgment_data.df, condition == "preemption")
preemption_judgment.df$semantically_correct <- factor(preemption_judgment.df$semantically_correct)
preemption_judgment.df$transitivity_test_scene2 <- factor(preemption_judgment.df$transitivity_test_scene2)
preemption_judgment.df$restricted_verb <- factor(preemption_judgment.df$restricted_verb)

Preregistered data analyses

Question 1: Have participants picked up on the difference in meaning between the two argument-structure constructions?

Production data

#Figure 2
RQ1_graph_productions.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph1 = aggregate(semantically_correct ~ verb_type_training2 + Participant.Private.ID, RQ1_graph_productions.df, FUN=mean)

aggregated.graph1 <- rename(aggregated.graph1, verb = verb_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = correct  ~ verb,
                  data = aggregated.graph1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "% semantically correct",
                  cex.lab = 1,
                  cex.axis = 1,
                  cex.names = 1,
                  yaxt = "n")

axis(2, at = seq(0, 1, by = 0.25), las=1)
abline(h = 0.50, lty = 2)

#1 alternating verb production

alternating_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating")

#and filter out responses where participants said something other than det1 or det2
alternating_prod.df = subset(alternating_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_alternating_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + Participant.Private.ID, alternating_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_alternating_prod.df$semantically_correct),3)

## [1] 0.945

# average accuracy separately for causative and inchoative scenes
round(tapply(aggregated.means_alternating_prod.df$semantically_correct, aggregated.means_alternating_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.971         0.919

# maximally vague priors for the intercept and the predictors
a = lizCenter(alternating_prod.df, list("transitivity_test_scene2"))  

alternating_model <-brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|Participant.Private.ID), data=a, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct" ))

##                                 Estimate Est.Error      Q2.5     Q97.5
## b_Intercept                    3.1731439 0.3566528  2.554203 3.9548012
## b_transitivity_test_scene2.ct -0.7512829 0.5142606 -1.766796 0.2593774

mcmc_plot(alternating_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##    c(C1, C2)
## 1 0.00000000
## 2 0.07216667

# no difference between construction 1 and construction 2

# Final model
# maximally vague priors for the intercept 
alternating_model_final = brm(formula = semantically_correct~1 + (1|Participant.Private.ID), data=a, family = bernoulli(link = logit),set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.025985 0.3396312 2.443645 3.768338

mcmc_plot(alternating_model_final, variable = "b_Intercept", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

#2 novel verb production

novel_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "novel")

#filter out responses where participants said something other than det1 or det2
novel_prod.df = subset(novel_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_novel_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + Participant.Private.ID, novel_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_novel_prod.df$semantically_correct),3)

## [1] 0.955

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_prod.df$semantically_correct, aggregated.means_novel_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.946         0.964

b = lizCenter(novel_prod.df, list("transitivity_test_scene2"))  

# maximally vague priors for the intercept and the predictors
novel_model <- brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|Participant.Private.ID), data=b, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct"))

##                                Estimate Est.Error      Q2.5    Q97.5
## b_Intercept                   3.4724866 0.3908453  2.782233 4.311712
## b_transitivity_test_scene2.ct 0.0618837 0.5954551 -1.144560 1.201520

mcmc_plot(novel_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.4531667

# no difference between construction 1 and construction 2  
# Final model

# maximally vague priors for the intercept 
novel_model_final <- brm(formula = semantically_correct~1+ (1|Participant.Private.ID), data=b, family = bernoulli(link = logit), set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.237372 0.3708822 2.594387 4.039616

mcmc_plot(novel_model_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

Judgment data

#Figure 3
RQ1_graph_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph2 = aggregate(Response ~ verb_type_training2 + semantically_correct + Participant.Private.ID, RQ1_graph_judgments.df, FUN=mean)
aggregated.graph2$semantically_correct <- recode(aggregated.graph2$semantically_correct, "1" = "yes","0" = "no")

aggregated.graph2 <- rename(aggregated.graph2, verb = verb_type_training2,
                                           correct = semantically_correct)

yarrr::pirateplot(formula = Response ~ correct + verb,
                  data = aggregated.graph2,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

#1 alternating verb judgments

alternating_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating")

# aggregated dataframe for means
aggregated.means_alternating_judgments = aggregate(Response ~ transitivity_test_scene2 + semantically_correct + Participant.Private.ID, alternating_judgments.df, FUN=mean)
aggregated.means_alternating_judgments$semantically_correct<- recode(aggregated.means_alternating_judgments$semantically_correct, "1" = "yes","0" = "no")
aggregated.means_alternating_judgments$transitivity_test_scene2<- recode(aggregated.means_alternating_judgments$transitivity_test_scene2, "construction1" = "transitive causative","construction2" = "intransitive inchoative")

# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_alternating_judgments$Response, aggregated.means_alternating_judgments$semantically_correct, mean),3)

##    no   yes 
## 2.506 4.837

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_alternating_judgments$Response, list(aggregated.means_alternating_judgments$semantically_correct, aggregated.means_alternating_judgments$transitivity_test_scene2), mean),3)

##     transitive causative intransitive inchoative
## no                 2.512                   2.500
## yes                4.779                   4.895

c = lizCenter(alternating_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments <-brm(formula = Response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|Participant.Private.ID), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(alternating_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct", "b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                        Estimate  Est.Error
## b_Intercept                                           3.6833521 0.07869591
## b_transitivity_test_scene2.ct                         0.0499209 0.06826174
## b_semantically_correct.ct                             2.2925296 0.12841265
## b_transitivity_test_scene2.ct:semantically_correct.ct 0.1254788 0.16133039
##                                                              Q2.5     Q97.5
## b_Intercept                                            3.52864089 3.8398844
## b_transitivity_test_scene2.ct                         -0.08445557 0.1831947
## b_semantically_correct.ct                              2.03776426 2.5410904
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.19210586 0.4418731

mcmc_plot(alternating_model_judgments, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] < 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1         0.0000000
## 2         0.2291667
## 3         0.0000000
## 4         0.2155000

# no difference between construction 1 and construction 2

# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments_final <-brm(formula = Response~semantically_correct.ct + (1 + semantically_correct.ct|Participant.Private.ID), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_judgments, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept               3.683352 0.07869591 3.528641 3.839884
## b_semantically_correct.ct 2.292530 0.12841265 2.037764 2.541090

mcmc_plot(alternating_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

#2 novel verb judgments

novel_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.means_novel_judgments = aggregate(Response ~ transitivity_test_scene2 + semantically_correct + Participant.Private.ID, novel_judgments.df, FUN=mean)


# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_novel_judgments$Response, aggregated.means_novel_judgments$semantically_correct, mean),3)

##     0     1 
## 2.250 4.174

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_judgments$Response, list(aggregated.means_novel_judgments$semantically_correct, aggregated.means_novel_judgments$transitivity_test_scene2), mean),3)

##   construction1 construction2
## 0         2.302         2.198
## 1         4.221         4.128

d = lizCenter(novel_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  

# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments <-brm(formula = Response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|Participant.Private.ID), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct", "b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                          Estimate  Est.Error
## b_Intercept                                            3.20303705 0.11356983
## b_transitivity_test_scene2.ct                         -0.09929811 0.09234945
## b_semantically_correct.ct                              1.87535088 0.16285313
## b_transitivity_test_scene2.ct:semantically_correct.ct  0.01177960 0.20139260
##                                                             Q2.5      Q97.5
## b_Intercept                                            2.9798663 3.42679777
## b_transitivity_test_scene2.ct                         -0.2786229 0.08256785
## b_semantically_correct.ct                              1.5565665 2.20330822
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.3815457 0.40806911

mcmc_plot(novel_model_judgments, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1         0.0000000
## 2         0.1409167
## 3         0.0000000
## 4         0.4764167

# no difference between construction 1 and construction 2
# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments_final <-brm(formula = Response~semantically_correct.ct + (1 + semantically_correct.ct|Participant.Private.ID), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_judgments, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept               3.683352 0.07869591 3.528641 3.839884
## b_semantically_correct.ct 2.292530 0.12841265 2.037764 2.541090

mcmc_plot(novel_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Question 2: Does statistical pre-emption constrain verb argument construction generalizations in adults (judgment data)?

#Figure 4

#first, filter our semantically incorrect trials

judgments_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel
judgments_novel.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items
judgment_unattested_constr1.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
judgment_unattested_constr2.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

judgment_unattested_novel.df <- rbind(judgments_novel.df, judgment_unattested_constr1.df, judgment_unattested_constr2.df)

aggregated.means = aggregate(Response ~ condition + restricted_verb + Participant.Private.ID, judgment_unattested_novel.df, FUN=mean)
aggregated.means<- rename(aggregated.means, restricted = restricted_verb)

yarrr::pirateplot(formula = Response ~ restricted + condition,
                  data = aggregated.means,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

judgment_unattested_novel_preemption.df <- subset(judgment_unattested_novel.df, condition == "preemption")   
round(tapply(judgment_unattested_novel_preemption.df$Response, judgment_unattested_novel_preemption.df$restricted_verb, mean),3)

##    no   yes 
## 3.026 2.362

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_preemption_model <- brm(formula = Response~(1 +restricted_verb.ct|Participant.Private.ID)+restricted_verb.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_preemption_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate  Est.Error       Q2.5      Q97.5
## b_Intercept           2.6927167 0.09992365  2.4924601  2.8901188
## b_restricted_verb.ct -0.6459548 0.16772805 -0.9756427 -0.3195383

mcmc_plot(judgments_preemption_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_preemption_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         1

# BF analyses: we use the difference between attested and unattested in the pilot study reported at rpubs.com/AnnaSamara/333562 as the maximum difference we could expect in comparing rating for unattested vs. novel constructions (SD = 3.15/2)
Bf(0.17, 0.65, uniform = 0, meanoftheory = 0, sdtheory = 3.15/2, tail = 1)

## $LikelihoodTheory
## [1] 0.4629748
## 
## $Likelihoodnull
## [1] 0.001570015
## 
## $BayesFactor
## [1] 294.8856

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.17, 0.65, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# RRs for which BF > 3
ev_for_h1 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.06

print(high_threshold)

## [1] 4

Question 3: Does statistical entrenchment constrain verb argument construction generalizations in adults (judgment data)?

#first, filter our semantically incorrect trials
entrenchment_judgment_unattested_novel.df <- subset(entrenchment_judgment.df, semantically_correct == "1")   

#we only want to keep novel

entrenchment_judgment_novel.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

entrenchment_judgment_unattested_constr1.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
entrenchment_judgment_unattested_constr2.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

entrenchment_judgment_unattested_novel.df <- rbind(entrenchment_judgment_novel.df, entrenchment_judgment_unattested_constr1.df, entrenchment_judgment_unattested_constr2.df)
entrenchment_judgment_unattested_novel.df$restricted_verb <- factor(entrenchment_judgment_unattested_novel.df$restricted_verb , levels = c("yes", "no"))

round(tapply(entrenchment_judgment_unattested_novel.df$Response, entrenchment_judgment_unattested_novel.df$restricted_verb, mean),3)

##   yes    no 
## 4.552 4.174

#Center variables of interest using the lizCenter function:
d_unattested_novel_entrenchment = lizCenter(entrenchment_judgment_unattested_novel.df , list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_entrenchment_model <- brm(formula = Response~(1 +restricted_verb.ct|Participant.Private.ID)+restricted_verb.ct, data=d_unattested_novel_entrenchment, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_entrenchment_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate Est.Error       Q2.5       Q97.5
## b_Intercept           4.3694021 0.1132104  4.1481279  4.59418338
## b_restricted_verb.ct -0.3665091 0.1629086 -0.6808387 -0.04604246

mcmc_plot(judgments_entrenchment_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_entrenchment_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

# the effect is in the opposite direction

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1     0.000
## 2     0.987

# drawing a max based on the difference between attested vs. unattested in this experiment 

# BF analyses: we use the difference between attested and unattested in this study (attested > unattested provides supporting evidence for entrenchment) as a maximum we expect when comparing ratings for unattested vs. novel constructions (SD = 0.38/2)
Bf(0.16, -0.36, uniform = 0, meanoftheory = 0, sdtheory = 0.38/2, tail = 1)

## $LikelihoodTheory
## [1] 0.0482932
## 
## $Likelihoodnull
## [1] 0.1983728
## 
## $BayesFactor
## [1] 0.2434466

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.16, -0.36, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# RRs for which BF <1/3
ev_for_h0 <- subset(data.frame(range_test), BF < 1/3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0

print(high_threshold)

## [1] 4

Question 4: Is the effect of statistical pre-emption larger than entrenchment (judgment data)?

#first, filter our semantically incorrect trials
all_judgment_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel

all_judgment_novel.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

all_judgment_unattested_constr1.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
all_judgment_unattested_constr2.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

all_judgment_unattested_novel.df <- rbind(all_judgment_novel.df, all_judgment_unattested_constr1.df, all_judgment_unattested_constr2.df)

round(tapply(all_judgment_unattested_novel.df$Response, list(all_judgment_unattested_novel.df$restricted_verb, all_judgment_unattested_novel.df$condition), mean),3)

##     entrenchment preemption
## no         4.174      3.026
## yes        4.552      2.362

#Center variables of interest using the lizCenter function:
df = lizCenter(all_judgment_unattested_novel.df, list("restricted_verb", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
judgments_pre_vs_ent_model <- brm(formula = Response~(1 +restricted_verb.ct|Participant.Private.ID)+restricted_verb.ct * condition.ct, data=df, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_pre_vs_ent_model, variable = c("b_Intercept", "b_restricted_verb.ct", "b_restricted_verb.ct:condition.ct"))

##                                    Estimate  Est.Error       Q2.5       Q97.5
## b_Intercept                        3.290080 0.07562601  3.1415574  3.43905880
## b_restricted_verb.ct              -0.286182 0.12011078 -0.5218846 -0.05475702
## b_restricted_verb.ct:condition.ct -1.002610 0.22863709 -1.4514705 -0.55007196

mcmc_plot(judgments_pre_vs_ent_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_pre_vs_ent_model))

C1=mean(samps[,"b_restricted_verb.ct"] > 0)
C2=mean(samps[,"b_restricted_verb.ct:condition.ct"] > 0) 
C3=mean(samps[,"b_Intercept"] > 0)

pMCMC=as.data.frame(c(C1,C2,C3))
pMCMC

##   c(C1, C2, C3)
## 1         0.009
## 2         0.000
## 3         1.000

#roughly predicted effect size from adult pilot study was 2.91. Use it as max for unattested vs. novel (SD = 2.91/2)
Bf(0.22, 0.99, uniform = 0, meanoftheory = 0, sdtheory = 2.91/2, tail = 1)

## $LikelihoodTheory
## [1] 0.4323971
## 
## $Likelihoodnull
## [1] 7.265337e-05
## 
## $BayesFactor
## [1] 5951.508

H1RANGE = seq(0,4,by=0.01) # [5-1]-[0] - max effect of preemption minus no effect of entrenchment
range_test <- Bf_range(0.22, 0.99, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.06

print(high_threshold)

## [1] 4

Exploratory data analyses

Effect of statistical pre-emption: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

# Figure 5
judgments_unattested_attested.df <- subset(combined_judgment_data.df, semantically_correct == "1")   
judgments_unattested_attested.df <- subset(judgments_unattested_attested.df, restricted_verb == "yes")   


aggregated.means1 = aggregate(Response ~ condition + attested_unattested + Participant.Private.ID, judgments_unattested_attested.df , FUN=mean)
aggregated.means1<- rename(aggregated.means1, attested = attested_unattested)

aggregated.means1$attested<- recode(aggregated.means1$attested, "1" = "yes","0" = "no")


yarrr::pirateplot(formula = Response ~  attested  + condition,
                  data = aggregated.means1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

# analyses
attested_vs_unattested = subset(preemption_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested$Response, attested_vs_unattested$attested_unattested, mean),3)

##     0     1 
## 2.362 4.952

#Center variables of interest using the lizCenter function:
d_attested_unattested = lizCenter(attested_vs_unattested , list("attested_unattested"))


# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_preemption <- brm(formula =Response~(1 +attested_unattested.ct|Participant.Private.ID)+attested_unattested.ct, data=d_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(attested_unattested_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate Est.Error     Q2.5    Q97.5
## b_Intercept              3.674561 0.0583548 3.560757 3.789909
## b_attested_unattested.ct 2.548472 0.1187866 2.309334 2.780839

mcmc_plot(attested_unattested_preemption, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

# this prior (5.00-1.85 = 3.15) is drawn from previous pilot study with 10 adults in preemption that was preregistered at https://rpubs.com/AnnaSamara/333562
Bf(0.12, 2.55, uniform = 0, meanoftheory = 0, sdtheory = 3.15, tail = 1)

## $LikelihoodTheory
## [1] 0.1824811
## 
## $Likelihoodnull
## [1] 2.92535e-98
## 
## $BayesFactor
## [1] 6.237925e+96

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.12, 2.55, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.01

print(high_threshold)

## [1] 4

Effect of statistical entrenchment: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

attested_vs_unattested_ent = subset(entrenchment_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_ent$Response, attested_vs_unattested_ent$attested_unattested, mean),3)

##     0     1 
## 4.552 4.942

#Center variables of interest using the lizCenter function:
d_attested_unattested_ent = lizCenter(attested_vs_unattested_ent, list("attested_unattested"))


# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment <- brm(formula =Response~(1 +attested_unattested.ct|Participant.Private.ID)+attested_unattested.ct, data=d_attested_unattested_ent, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate  Est.Error      Q2.5     Q97.5
## b_Intercept              4.7495366 0.06132368 4.6302263 4.8692174
## b_attested_unattested.ct 0.3828401 0.12359019 0.1384787 0.6260996

mcmc_plot(attested_unattested_entrenchment, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_entrenchment))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.00175

# we preregistered that the max effect of entrenchment here would be 1 based on adult data suggesting difference never more than 1 with novel verbs, i.e. SD = 0.5
Bf(0.13, 0.38, uniform = 0, meanoftheory = 0, sdtheory = 0.5, tail = 1)

## $LikelihoodTheory
## [1] 1.175537
## 
## $Likelihoodnull
## [1] 0.04281328
## 
## $BayesFactor
## [1] 27.45731

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.13, 0.38, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.06

print(high_threshold)

## [1] 4

Entrenchment vs. preemption: ratings for witnessed vs. unwitnessed forms

attested_vs_unattested_across = subset(combined_judgment_data.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_across$Response, list(attested_vs_unattested_across$condition, attested_vs_unattested_across$attested_unattested), mean),3)

##                  0     1
## entrenchment 4.552 4.942
## preemption   2.362 4.952

#Center variables of interest using the lizCenter function:
df_attested_unattested = lizCenter(attested_vs_unattested_across, list("attested_unattested", "condition"))


# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment_preemption <- brm(formula = Response~(1 +attested_unattested.ct|Participant.Private.ID)+attested_unattested.ct * condition.ct, data=df_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment_preemption, variable = c("b_Intercept","b_condition.ct", "b_attested_unattested.ct","b_attested_unattested.ct:condition.ct"))

##                                        Estimate  Est.Error      Q2.5      Q97.5
## b_Intercept                            4.055831 0.04442185  3.968996  4.1424373
## b_condition.ct                        -1.049018 0.08525871 -1.214662 -0.8800325
## b_attested_unattested.ct               1.781977 0.09074218  1.606896  1.9589410
## b_attested_unattested.ct:condition.ct  2.109871 0.17228987  1.769255  2.4503766

mcmc_plot(attested_unattested_entrenchment_preemption, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_entrenchment_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_condition.ct"] < 0)
C3=mean(samps[,"b_attested_unattested.ct"] < 0) 
C4=mean(samps[,"b_attested_unattested.ct:condition.ct"] < 0) 

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1                 0
## 2                 1
## 3                 0
## 4                 0

#roughly predicted effect size from adult pilot study: 2.91
Bf(0.17, 2.11, uniform = 0, meanoftheory = 0, sdtheory = 2.91, tail = 1)

## $LikelihoodTheory
## [1] 0.210635
## 
## $Likelihoodnull
## [1] 8.289253e-34
## 
## $BayesFactor
## [1] 2.541061e+32

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.17, 2.11, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.02

print(high_threshold)

## [1] 4

Production data: Effect of statistical pre-emption

#Figure 6

data_long <- gather(preemption_production.df, det_type, produced, det1:none, factor_key=TRUE)

p = ggplot(data_long, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme_bw()+
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "periphrastic causative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "periphrastic causative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted verbs against chance 

production_preemption_attested_unattested.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, restricted_verb =="yes")

round(tapply(production_preemption_attested_unattested.df $attested_unattested, production_preemption_attested_unattested.df $verb_type_training2, mean),3)

## construction1 construction2 
##         0.994         0.994

production_preemption_attested_unattested.df$verb_type_training2 <- factor(production_preemption_attested_unattested.df$verb_type_training2)

df_prod = lizCenter(production_preemption_attested_unattested.df , list("verb_type_training2"))  

# maximally vague priors for the predictors and the intercept
prod_attested_unattested = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|Participant.Private.ID), data=df_prod, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested, variable = c("b_Intercept","b_verb_type_training2.ct"))

##                            Estimate Est.Error      Q2.5    Q97.5
## b_Intercept              4.66652282 0.4196348  3.919495 5.565043
## b_verb_type_training2.ct 0.01178403 0.6257371 -1.236330 1.239773

mcmc_plot(prod_attested_unattested, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.4896667

#same analyses without verb_training_type

# maximally vague priors for the intercept
prod_attested_unattested_final = brm(formula = attested_unattested ~1 + (1|Participant.Private.ID), data=df_prod, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 4.511625 0.3851926 3.832252 5.322255

mcmc_plot(prod_attested_unattested_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_final))
C1=mean(samps[,"b_Intercept"] < 0)


# We will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the restricted verbs than for the novel verb


production_preemption_restricted_novel.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_restricted_novel.df<- subset(production_preemption_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_preemption_restricted_novel.df$attested_unattested)

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_preemption_restricted_novel.df$attested_unattested)

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.994         0.994         0.470

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.470 0.994

production_preemption_restricted_novel.df$restricted_verb <- factor(production_preemption_restricted_novel.df$restricted_verb)
production_preemption_restricted_novel1.df = lizCenter(production_preemption_restricted_novel.df, list("restricted_verb"))

# maximally vague priors for the predictors and the intercept
prod_unattested_novel_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|Participant.Private.ID), data=production_preemption_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                      Estimate Est.Error     Q2.5    Q97.5
## b_Intercept          3.253954 0.3456579 2.628447 4.003758
## b_restricted_verb.ct 3.792121 0.6640449 2.410433 5.002237

mcmc_plot(prod_unattested_novel_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_unattested_novel_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Production data: Effect of statistical entrenchment

data_long_e <- gather(entrenchment_production.df, det_type, produced, det1:none, factor_key=TRUE)

data_long_e$transitivity_test_scene2 <-recode(data_long_e$transitivity_test_scene2, "construction1" = "test: transitive causative","construction2" = "test: intransitive inchoative")


p = ggplot(data_long_e, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme(panel.grid.major = element_blank()) +
  facet_grid("transitivity_test_scene2") +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "intransitive inchoative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "intransitive inchoative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#a. Are participants producing more attested than unattested dets?
# here, we want to see how often participants say the unattested e.g. transitive-only det1 for a det2 (intransitive-only) verb in the intransitive condition at test 
# and vice versa 

production_entrenchment_attested_unattested.df  <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, restricted_verb =="yes")

#We want to compare attested vs. unattested trials for transitive verbs in the intransitive inchoative construction at test
production_entrenchment_attested_unattested1.df  <- subset(production_entrenchment_attested_unattested.df, verb_type_training2 == "construction1" & transitivity_test_scene2 == "construction2")

#And intransitive inchoative verbs in the transitive construction at test. Filter out irrelevant trials
production_entrenchment_attested_unattested2.df  <- subset(production_entrenchment_attested_unattested.df, verb_type_training2 == "construction2" & transitivity_test_scene2 == "construction1")


production_entrenchment_attested_unattested.df <- rbind(production_entrenchment_attested_unattested1.df, production_entrenchment_attested_unattested2.df)

#How much of the time are participants producing attested items?
round(mean(production_entrenchment_attested_unattested.df$attested_unattested),3)

## [1] 0.148

# and separately for each verb type
round(tapply(production_entrenchment_attested_unattested.df$attested_unattested, production_entrenchment_attested_unattested.df$verb_type_training2, mean),3)

## construction1 construction2 
##         0.174         0.122

production_entrenchment_attested_unattested.df$verb_type_training2 <- factor(production_entrenchment_attested_unattested.df$verb_type_training2)
df_prod_ent = lizCenter((production_entrenchment_attested_unattested.df), list("verb_type_training2"))  


# maximally vague priors for the predictors and the intercept
prod_attested_unattested_ent = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|Participant.Private.ID), data=df_prod_ent, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent, variable = c("b_Intercept","b_verb_type_training2.ct"))

##                            Estimate Est.Error      Q2.5      Q97.5
## b_Intercept              -2.7140164 0.4912616 -3.705935 -1.7837384
## b_verb_type_training2.ct -0.6347692 0.4664496 -1.556009  0.2897334

mcmc_plot(prod_attested_unattested_ent, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##    c(C1, C2)
## 1 1.00000000
## 2 0.08683333

#same analyses without verb_training_type


# maximally vague priors for the intercept
prod_attested_unattested_ent_final = brm(formula = attested_unattested ~1 + (1|Participant.Private.ID), data=df_prod_ent, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

summary(prod_attested_unattested_ent_final, WAIC=T)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: attested_unattested ~ 1 + (1 | Participant.Private.ID) 
##    Data: df_prod_ent (Number of observations: 344) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2000; thin = 1;
##          total post-warmup draws = 12000
## 
## Group-Level Effects: 
## ~Participant.Private.ID (Number of levels: 43) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     2.82      0.64     1.82     4.31 1.00     3341     4636
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -2.65      0.48    -3.62    -1.74 1.00     3909     4682
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

posterior_summary(prod_attested_unattested_ent_final, variable = c("b_Intercept"))

##              Estimate Est.Error      Q2.5     Q97.5
## b_Intercept -2.650557 0.4759355 -3.616911 -1.739386

mcmc_plot(prod_attested_unattested_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 1

# c. we will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the 2 non-alternating verbs than for the novel verb (presumably the “unwitnessed” form has to be set arbitrarily here)


production_entrenchment_restricted_novel.df <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_restricted_novel.df<- subset(production_entrenchment_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_entrenchment_restricted_novel.df$attested_unattested)
production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_entrenchment_restricted_novel.df$attested_unattested)

# select trials featuring the novel verb in the intransitive inchoative construction
production_entrenchment_restricted_novel1.df <- subset(production_entrenchment_restricted_novel.df, verb_type_training2 == "novel"  & transitivity_test_scene2 == "construction2")


# Select trials featuring transitive verbs in the intransitive inchoative construction at test
production_entrenchment_restricted_novel2.df  <- subset(production_entrenchment_restricted_novel.df, verb_type_training2 == "construction1" & transitivity_test_scene2 == "construction2")

# Select trials featuring intransitive verbs in the transitive construction at test
production_entrenchment_restricted_novel3.df  <- subset(production_entrenchment_restricted_novel.df, verb_type_training2 == "construction2" & transitivity_test_scene2 == "construction1")


production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.174         0.122         0.036

# reverse coding to focus on unattested rather than attested for novel vs. restricted
production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)
production_entrenchment_restricted_novel.df$attested_unattested<- recode(production_entrenchment_restricted_novel.df$attested_unattested, `1` = 0L, `0` = 1L)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.964 0.852

#what this means is that participants produce *unattested forms* less for the restricted than they do for the novel

production_entrenchment_restricted_novel.df$restricted_verb <- factor(production_entrenchment_restricted_novel.df$restricted_verb)
production_entrenchment_restricted_novel1.df = lizCenter(production_entrenchment_restricted_novel.df, list("restricted_verb"))


# maximally vague priors for the predictors and the intercept
prod_unattested_novel_ent_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|Participant.Private.ID), data=production_entrenchment_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_ent_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                        Estimate Est.Error      Q2.5     Q97.5
## b_Intercept           3.1474443 0.3975774  2.406269 3.9667440
## b_restricted_verb.ct -0.7137785 0.6248564 -1.972300 0.4830665

mcmc_plot(prod_unattested_novel_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_unattested_novel_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.12375

Experiment 2

Load data

#Create the dataframes that we will be working on
combined_production_data.df <- read.csv("exp2_production_data.csv")

combined_judgment_data.df <- read.csv("exp2_judgment_data.csv")
combined_judgment_data.df$restricted_verb <- factor(combined_judgment_data.df$restricted_verb)
combined_judgment_data.df$condition <- factor(combined_judgment_data.df$condition)


#separately for entrenchment and preemption

#entrenchment
entrenchment_production.df <- subset(combined_production_data.df, condition == "entrenchment")
entrenchment_production.df$semantically_correct <- as.numeric(entrenchment_production.df$semantically_correct)
entrenchment_production.df$transitivity_test_scene2 <- factor(entrenchment_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in entrenchment
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

entrenchment_production.df$det1 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction1", 1, 0)
entrenchment_production.df$det2 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction2", 1, 0)
entrenchment_production.df$other <- ifelse(entrenchment_production.df$det_lenient_adapted == "other", 1, 0)
entrenchment_production.df$none <- ifelse(entrenchment_production.df$det_lenient_adapted == "none", 1, 0)


entrenchment_judgment.df <- subset(combined_judgment_data.df, condition == "entrenchment")
entrenchment_judgment.df$semantically_correct <- factor(entrenchment_judgment.df$semantically_correct)
entrenchment_judgment.df$transitivity_test_scene2 <- factor(entrenchment_judgment.df$transitivity_test_scene2)
entrenchment_judgment.df$restricted_verb <- factor(entrenchment_judgment.df$restricted_verb)


#preemption
preemption_production.df <- subset(combined_production_data.df, condition == "preemption")
preemption_production.df$semantically_correct <- as.numeric(preemption_production.df$semantically_correct)  #actually, all NAs here
preemption_production.df$transitivity_test_scene2 <- factor(preemption_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in pre-emption
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

preemption_production.df$det1 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction1", 1, 0)
preemption_production.df$det2 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction2", 1, 0)
preemption_production.df$other <- ifelse(preemption_production.df$det_lenient_adapted == "other", 1, 0)
preemption_production.df$none <- ifelse(preemption_production.df$det_lenient_adapted == "none", 1, 0)


preemption_judgment.df <- subset(combined_judgment_data.df, condition == "preemption")
preemption_judgment.df$semantically_correct <- factor(preemption_judgment.df$semantically_correct)
preemption_judgment.df$transitivity_test_scene2 <- factor(preemption_judgment.df$transitivity_test_scene2)
preemption_judgment.df$restricted_verb <- factor(preemption_judgment.df$restricted_verb)

Preregistered data analyses

Question 1: Have participants picked up on the difference in meaning between the two argument-structure constructions?

Production data

#Figure 7

RQ1_graph_productions.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph1 = aggregate(semantically_correct ~ verb_type_training2 + participant_private_id, RQ1_graph_productions.df, FUN=mean)

aggregated.graph1 <- rename(aggregated.graph1, verb = verb_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = correct  ~ verb,
                  data = aggregated.graph1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "% semantically correct",
                  cex.lab = 1,
                  cex.axis = 1,
                  cex.names = 1,
                  yaxt = "n")

axis(2, at = seq(0, 1, by = 0.25), las=1)
abline(h = 0.50, lty = 2)

#1 alternating verb production

alternating_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating")

#and filter out responses where participants said something other than det1 or det2
alternating_prod.df = subset(alternating_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_alternating_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + participant_private_id, alternating_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_alternating_prod.df$semantically_correct),3)

## [1] 0.975

# average accuracy separately for causative and inchoative scenes
round(tapply(aggregated.means_alternating_prod.df$semantically_correct, aggregated.means_alternating_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.963         0.988

a = lizCenter(alternating_prod.df, list("transitivity_test_scene2"))   

# maximally vague priors for the intercept and the predictors
alternating_model <- brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|participant_private_id), data=a, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct" ))

##                                Estimate Est.Error       Q2.5    Q97.5
## b_Intercept                   3.6240014 0.3711631  2.9720934 4.419299
## b_transitivity_test_scene2.ct 0.5169631 0.5830616 -0.6092265 1.667905

mcmc_plot(alternating_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.81225

# no difference between construction 1 and construction 2

# Final model
# maximally vague priors for the intercept 
alternating_model_final = brm(formula = semantically_correct~1 + (1|participant_private_id), data=a, family = bernoulli(link = logit),set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error    Q2.5    Q97.5
## b_Intercept 3.451023 0.3381242 2.84819 4.177323

mcmc_plot(alternating_model_final, variable = "b_Intercept", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

#2 novel verb production

novel_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "novel")

#and filter out responses where participants said something other than det1 or det2
novel_prod.df = subset(novel_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_novel_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + participant_private_id, novel_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_novel_prod.df$semantically_correct),3)

## [1] 0.96

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_prod.df$semantically_correct, aggregated.means_novel_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.975         0.946

b = lizCenter(novel_prod.df, list("transitivity_test_scene2"))  


# maximally vague priors for the intercept and the predictors
novel_model <- brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|participant_private_id), data=b, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct"))

##                                 Estimate Est.Error      Q2.5     Q97.5
## b_Intercept                    3.6346413 0.4290942  2.860954 4.5508293
## b_transitivity_test_scene2.ct -0.3121872 0.6075546 -1.504537 0.9067366

mcmc_plot(novel_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.7031667

# no difference between construction 1 and construction 2  
# Final model

# maximally vague priors for the intercept 
novel_model_final <- brm(formula = semantically_correct~1+ (1 + 1|participant_private_id), data=b, family = bernoulli(link = logit), set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.486862  0.412819 2.753102 4.372883

mcmc_plot(novel_model_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

Judgment data

#Figure 8
RQ1_graph_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph2 = aggregate(response ~ verb_type_training2+ semantically_correct + participant_private_id, RQ1_graph_judgments.df, FUN=mean)
aggregated.graph2$semantically_correct <- recode(aggregated.graph2$semantically_correct, "1" = "yes","0" = "no")

aggregated.graph2 <- rename(aggregated.graph2, verb = verb_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = response ~ correct + verb,
                  data = aggregated.graph2,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

alternating_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating")

#1 alternating verb judgments

# aggregated dataframe for means
aggregated.means_alternating_judgments = aggregate(response ~ transitivity_test_scene2 + semantically_correct + participant_private_id, alternating_judgments.df, FUN=mean)

# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_alternating_judgments$response, aggregated.means_alternating_judgments$semantically_correct, mean),3)

##     0     1 
## 2.294 4.775

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_alternating_judgments$response, list(aggregated.means_alternating_judgments$semantically_correct, aggregated.means_alternating_judgments$transitivity_test_scene2), mean),3)

##   construction1 construction2
## 0         2.350         2.237
## 1         4.812         4.737

c = lizCenter(alternating_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments <-brm(formula = response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|participant_private_id), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct","b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                          Estimate  Est.Error
## b_Intercept                                            3.54612840 0.07060380
## b_transitivity_test_scene2.ct                         -0.09279547 0.07879941
## b_semantically_correct.ct                              2.43230358 0.14190796
## b_transitivity_test_scene2.ct:semantically_correct.ct  0.03688559 0.15197141
##                                                             Q2.5      Q97.5
## b_Intercept                                            3.4085529 3.68650734
## b_transitivity_test_scene2.ct                         -0.2472334 0.06333402
## b_semantically_correct.ct                              2.1468487 2.70419033
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.2618707 0.33391845

samps = as.matrix(as.mcmc(alternating_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1         0.0000000
## 2         0.1205000
## 3         0.0000000
## 4         0.4018333

# no difference between construction 1 and construction 2

# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments_final <-brm(formula = response~semantically_correct.ct + (1 + semantically_correct.ct|participant_private_id), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(alternating_model_judgments_final, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept               3.548481 0.07085625 3.409274 3.687733
## b_semantically_correct.ct 2.431171 0.14062245 2.153712 2.699466

mcmc_plot(alternating_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

#2 novel verb judgments

novel_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.means_novel_judgments = aggregate(response ~ transitivity_test_scene2 + semantically_correct + participant_private_id, novel_judgments.df, FUN=mean)

# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_novel_judgments$response, aggregated.means_novel_judgments$semantically_correct, mean),3)

##     0     1 
## 2.163 3.944

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_judgments$response, list(aggregated.means_novel_judgments$semantically_correct, aggregated.means_novel_judgments$transitivity_test_scene2), mean),3)

##   construction1 construction2
## 0         2.188         2.138
## 1         3.938         3.950

d = lizCenter(novel_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  


# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments <-brm(formula = response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|participant_private_id), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct","b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                          Estimate  Est.Error
## b_Intercept                                            3.04288435 0.11071576
## b_transitivity_test_scene2.ct                         -0.02139120 0.09471911
## b_semantically_correct.ct                              1.72392627 0.18029796
## b_transitivity_test_scene2.ct:semantically_correct.ct  0.05942706 0.19077868
##                                                             Q2.5     Q97.5
## b_Intercept                                            2.8259439 3.2635877
## b_transitivity_test_scene2.ct                         -0.2062836 0.1623206
## b_semantically_correct.ct                              1.3697743 2.0783839
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.3101031 0.4366050

mcmc_plot(novel_model_judgments, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1         0.0000000
## 2         0.4081667
## 3         0.0000000
## 4         0.3770000

# no difference between construction 1 and construction 2
# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments_final <-brm(formula = response~semantically_correct.ct + (1 + semantically_correct.ct|participant_private_id), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(novel_model_judgments_final, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate Est.Error     Q2.5    Q97.5
## b_Intercept               3.040811 0.1117169 2.817727 3.254773
## b_semantically_correct.ct 1.721705 0.1838920 1.353737 2.083251

mcmc_plot(novel_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Question 2: Does statistical pre-emption constrain verb argument construction generalizations in adults (judgment data)?

#Figure 9

#first, filter our semantically incorrect trials

judgments_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel

judgments_novel.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

judgment_unattested_constr1.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
judgment_unattested_constr2.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

judgment_unattested_novel.df <- rbind(judgments_novel.df, judgment_unattested_constr1.df, judgment_unattested_constr2.df)

aggregated.means = aggregate(response ~ condition + restricted_verb + participant_private_id, judgment_unattested_novel.df, FUN=mean)
aggregated.means<- rename(aggregated.means, restricted = restricted_verb)

yarrr::pirateplot(formula = response ~ restricted + condition,
                  data = aggregated.means,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

judgment_unattested_novel_preemption.df <- subset(judgment_unattested_novel.df, condition == "preemption")   
round(tapply(judgment_unattested_novel_preemption.df$response, judgment_unattested_novel_preemption.df$restricted_verb, mean),3)

##    no   yes 
## 2.644 2.259

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_preemption_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(judgments_preemption_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate  Est.Error       Q2.5      Q97.5
## b_Intercept           2.4515551 0.09074484  2.2699686  2.6281726
## b_restricted_verb.ct -0.3704174 0.18317117 -0.7348612 -0.0107444

mcmc_plot(judgments_preemption_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_preemption_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.9775833

# BF analyses: we use the difference between attested and unattested in Experiment1 (SD = 0.65) as an estimate of the difference we expect in comparing rating for unattested vs. novel constructions 
Bf(0.18, 0.37, uniform = 0, meanoftheory = 0, sdtheory = 0.65, tail = 1)

## $LikelihoodTheory
## [1] 0.9929748
## 
## $Likelihoodnull
## [1] 0.267993
## 
## $BayesFactor
## [1] 3.705227

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.18, 0.37, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# RRs for which BF > 3
ev_for_h1 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.14

print(high_threshold)

## [1] 0.87

Question 3: Does statistical entrenchment constrain verb argument construction generalizations in adults (judgment data)?

#first, filter our semantically incorrect trials

entrenchment_judgment_unattested_novel.df <- subset(entrenchment_judgment.df, semantically_correct == "1")   

#we only want to keep novel

entrenchment_judgment_novel.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

entrenchment_judgment_unattested_constr1.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
entrenchment_judgment_unattested_constr2.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

entrenchment_judgment_unattested_novel.df <- rbind(entrenchment_judgment_novel.df, entrenchment_judgment_unattested_constr1.df, entrenchment_judgment_unattested_constr2.df)
entrenchment_judgment_unattested_novel.df$restricted_verb <- factor(entrenchment_judgment_unattested_novel.df$restricted_verb , levels = c("yes", "no"))

round(tapply(entrenchment_judgment_unattested_novel.df$response, entrenchment_judgment_unattested_novel.df$restricted_verb, mean),3)

##   yes    no 
## 4.356 3.944

#Center variables of interest using the lizCenter function:
d_unattested_novel_entrenchment = lizCenter(entrenchment_judgment_unattested_novel.df, list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_entrenchment_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct, data=d_unattested_novel_entrenchment, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_entrenchment_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate Est.Error       Q2.5       Q97.5
## b_Intercept           4.1534153 0.1271332  3.9006063  4.40281833
## b_restricted_verb.ct -0.4001877 0.1704516 -0.7337064 -0.06430223

mcmc_plot(judgments_entrenchment_model, variable = "b_Intercept", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_entrenchment_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.9893333

# drawing a max based on the difference between attested vs. unattested in this experiment (this was sig. evidence for entrenchment)
Bf(0.17, -0.40, uniform = 0, meanoftheory = 0, sdtheory = 0.38/2, tail = 1)

## $LikelihoodTheory
## [1] 0.03663446
## 
## $Likelihoodnull
## [1] 0.1473201
## 
## $BayesFactor
## [1] 0.2486726

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.17, -0.40, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF < 3
ev_for_h1 <- subset(data.frame(range_test), BF < 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0

print(high_threshold)

## [1] 4

Question 4: Is the effect of statistical pre-emption larger than entrenchment (judgment data)?

#first, filter our semantically incorrect trials

all_judgment_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel

all_judgment_novel.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

all_judgment_unattested_constr1.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
all_judgment_unattested_constr2.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

all_judgment_unattested_novel.df <- rbind(all_judgment_novel.df, all_judgment_unattested_constr1.df, all_judgment_unattested_constr2.df)

round(tapply(all_judgment_unattested_novel.df$response, list(all_judgment_unattested_novel.df$restricted_verb, all_judgment_unattested_novel.df$condition), mean),3)

##     entrenchment preemption
## no         3.944      2.644
## yes        4.356      2.259

# preemption worked and opposite effect for entrenchment

#Center variables of interest using the lizCenter function:
df = lizCenter(all_judgment_unattested_novel.df, list("restricted_verb", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
judgments_pre_vs_ent_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct * condition.ct, data=df, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_pre_vs_ent_model, variable = c("b_Intercept", "b_restricted_verb.ct","b_condition.ct", "b_restricted_verb.ct:condition.ct"))

##                                     Estimate  Est.Error       Q2.5      Q97.5
## b_Intercept                        3.0187066 0.07843126  2.8630026  3.1717817
## b_restricted_verb.ct              -0.1098356 0.13089887 -0.3688386  0.1533700
## b_condition.ct                    -1.6678278 0.14827473 -1.9547380 -1.3734727
## b_restricted_verb.ct:condition.ct -0.7573872 0.24531036 -1.2402619 -0.2727975

mcmc_plot(judgments_pre_vs_ent_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_pre_vs_ent_model))

C1=mean(samps[,"b_Intercept"] < 0) 
C2=mean(samps[,"b_restricted_verb.ct"] > 0)
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_restricted_verb.ct:condition.ct"] > 0) 


pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1            0.0000
## 2            0.1965
## 3            0.0000
## 4            0.0015

#roughly predicted effect size from previous study was 1.0. Use it as an estimate of the effect we expect here
Bf(0.24, 0.77, uniform = 0, meanoftheory = 0, sdtheory = 1.00, tail = 1)

## $LikelihoodTheory
## [1] 0.5856484
## 
## $Likelihoodnull
## [1] 0.009671967
## 
## $BayesFactor
## [1] 60.55112

H1RANGE = seq(0,4,by=0.01) # [5-1]-[0] - max effect of preemption minus no effect of entrenchment
range_test <- Bf_range(0.24, 0.77, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.09

print(high_threshold)

## [1] 4

Exploratory data analyses

Effect of statistical pre-emption: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

# Figure 10
judgments_unattested_attested.df <- subset(combined_judgment_data.df, semantically_correct == "1")   
judgments_unattested_attested.df <- subset(judgments_unattested_attested.df, restricted_verb == "yes")   


aggregated.means1 = aggregate(response ~ condition + attested_unattested + participant_private_id, judgments_unattested_attested.df , FUN=mean)
aggregated.means1<- rename(aggregated.means1, attested = attested_unattested)

aggregated.means1$attested<- recode(aggregated.means1$attested, "1" = "yes","0" = "no")


yarrr::pirateplot(formula = response ~  attested  + condition,
                  data = aggregated.means1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

# analyses
attested_vs_unattested = subset(preemption_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested$response, attested_vs_unattested$attested_unattested, mean),3)

##     0     1 
## 2.259 4.881

#Center variables of interest using the lizCenter function:
d_attested_unattested = lizCenter(attested_vs_unattested , list("attested_unattested"))


# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_preemption <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(attested_unattested_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept              3.582812 0.05517553 3.476109 3.693045
## b_attested_unattested.ct 2.586153 0.11775131 2.353747 2.816768

mcmc_plot(attested_unattested_preemption, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

# this prioris drawn from Experiment 1
Bf(0.12, 2.59, uniform = 0, meanoftheory = 0, sdtheory = 2.55, tail = 1)

## $LikelihoodTheory
## [1] 0.1868105
## 
## $Likelihoodnull
## [1] 2.321655e-101
## 
## $BayesFactor
## [1] 8.046437e+99

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.12, 2.59, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.01

print(high_threshold)

## [1] 4

Effect of statistical entrenchment: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

attested_vs_unattested_ent = subset(entrenchment_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_ent$response, attested_vs_unattested_ent$attested_unattested, mean),3)

##     0     1 
## 4.356 4.925

#Center variables of interest using the lizCenter function:
d_attested_unattested_ent = lizCenter(attested_vs_unattested_ent, list("attested_unattested"))

attested_unattested_entrenchment <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested_ent, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate  Est.Error      Q2.5     Q97.5
## b_Intercept              4.6440445 0.07739698 4.4930217 4.7973083
## b_attested_unattested.ct 0.5592226 0.14355854 0.2752635 0.8401481

mcmc_plot(attested_unattested_entrenchment, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_entrenchment))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##      c(C1, C2)
## 1 0.0000000000
## 2 0.0004166667

# expect a difference of 0.38 from previous work
Bf(0.14, 0.56, uniform = 0, meanoftheory = 0, sdtheory = 0.38, tail = 1)

## $LikelihoodTheory
## [1] 0.7572747
## 
## $Likelihoodnull
## [1] 0.0009559302
## 
## $BayesFactor
## [1] 792.1862

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.14, 0.56, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.04

print(high_threshold)

## [1] 4

Entrenchment vs. preemption: ratings for witnessed vs. unwitnessed forms

attested_vs_unattested_across = subset(combined_judgment_data.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_across$response, list(attested_vs_unattested_across$condition, attested_vs_unattested_across$attested_unattested), mean),3)

##                  0     1
## entrenchment 4.356 4.925
## preemption   2.259 4.881

#Center variables of interest using the lizCenter function:
df_attested_unattested = lizCenter(attested_vs_unattested_across, list("attested_unattested", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment_preemption <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct * condition.ct, data=df_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment_preemption, variable = c("b_Intercept", "b_attested_unattested.ct",  "b_condition.ct", "b_attested_unattested.ct:condition.ct"))

##                                        Estimate  Est.Error      Q2.5      Q97.5
## b_Intercept                            3.938438 0.04626309  3.846235  4.0300119
## b_attested_unattested.ct               1.911223 0.09244683  1.728814  2.0931752
## b_condition.ct                        -1.034068 0.08943181 -1.205561 -0.8555843
## b_attested_unattested.ct:condition.ct  1.968536 0.17904181  1.621374  2.3122350

mcmc_plot(attested_unattested_entrenchment_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_attested_unattested.ct:condition.ct"] < 0) 

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1                 0
## 2                 0
## 3                 0
## 4                 0

#roughly predicted effect size from previous study 2.11

Bf(0.18, 1.97, uniform = 0, meanoftheory = 0, sdtheory = 2.11, tail = 1)

## $LikelihoodTheory
## [1] 0.2444349
## 
## $Likelihoodnull
## [1] 2.165476e-26
## 
## $BayesFactor
## [1] 1.128781e+25

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.18, 1.97, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.02

print(high_threshold)

## [1] 4

Production data: Effect of statistical pre-emption

#Figure 11

data_long <- gather(preemption_production.df, det_type, produced, det1:none, factor_key=TRUE)

p = ggplot(data_long, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme_bw()+
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "periphrastic causative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "periphrastic causative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_verb"))


#a. Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted verbs against chance 

production_preemption_attested_unattested.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, restricted_verb =="yes")

round(tapply(production_preemption_attested_unattested.df$attested_unattested, production_preemption_attested_unattested.df$verb_type_training2, mean),3)

## construction1 construction2 
##         0.997         0.990

production_preemption_attested_unattested.df$verb_type_training2 <- factor(production_preemption_attested_unattested.df$verb_type_training2)

df_prod = lizCenter(production_preemption_attested_unattested.df , list("verb_type_training2"))  

# maximally vague priors for the predictors and the intercept
prod_attested_unattested = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|participant_private_id), data=df_prod, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested, variable = c("b_Intercept","b_verb_type_training2.ct"))

##                           Estimate Est.Error      Q2.5     Q97.5
## b_Intercept               4.586672 0.3980398  3.875482 5.4411610
## b_verb_type_training2.ct -0.300786 0.6179178 -1.532312 0.9139024

mcmc_plot(prod_attested_unattested, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.3151667

#same analyses without verb_training_type

# maximally vague priors for the intercept
prod_attested_unattested_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 4.436624 0.3634798 3.786213 5.206571

mcmc_plot(prod_attested_unattested_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

# We will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the restricted verbs than for the novel verb

production_preemption_restricted_novel.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_restricted_novel.df<- subset(production_preemption_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_preemption_restricted_novel.df$attested_unattested)

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_preemption_restricted_novel.df$attested_unattested)

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.997         0.990         0.478

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.478 0.993

production_preemption_restricted_novel.df$restricted_verb <- factor(production_preemption_restricted_novel.df$restricted_verb)
production_preemption_restricted_novel1.df = lizCenter(production_preemption_restricted_novel.df, list("restricted_verb"))

# maximally vague priors for the predictors and the intercept
prod_unattested_novel_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|participant_private_id), data=production_preemption_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                      Estimate Est.Error     Q2.5    Q97.5
## b_Intercept          3.062048 0.2654188 2.575785 3.624023
## b_restricted_verb.ct 3.876751 0.5355444 2.774008 4.879432

mcmc_plot(prod_unattested_novel_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_unattested_novel_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Production data: Effect of statistical entrenchment

#Figure 12

# here, we want to see how often participants say det1 (e.g. transitive-only) for a det2 (intransitive-only) verb in the intransitive condition at test. 
# in other words, do they produce the unattested in a semantically correct trial?  

data_long_e <- gather(entrenchment_production.df, det_type, produced, det1:none, factor_key=TRUE)

data_long_e$transitivity_test_scene2 <-recode(data_long_e$transitivity_test_scene2, "construction1" = "test: transitive causative","construction2" = "test: intransitive inchoative")


p = ggplot(data_long_e, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme(panel.grid.major = element_blank()) +
  facet_grid("transitivity_test_scene2") +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "intransitive inchoative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "intransitive inchoative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#a. Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted verbs against chance 

production_entrenchment_attested_unattested.df  <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, restricted_verb =="yes")

#How much of the time are participants producing attested items?
round(mean(production_entrenchment_attested_unattested.df$attested_unattested),3)

## [1] 0.56

# and separately for each verb type
round(tapply(production_entrenchment_attested_unattested.df$attested_unattested, production_entrenchment_attested_unattested.df$verb_type_training2, mean),3)

## construction1 construction2 
##         0.571         0.549

production_entrenchment_attested_unattested.df$verb_type_training2 <- factor(production_entrenchment_attested_unattested.df$verb_type_training2)
df_prod_ent = lizCenter((production_entrenchment_attested_unattested.df), list("verb_type_training2"))  


# maximally vague priors for the predictors and the intercept
prod_attested_unattested_ent = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent, variable = c("b_Intercept", "b_verb_type_training2.ct"))

##                             Estimate Est.Error        Q2.5     Q97.5
## b_Intercept               0.24175335 0.0820546  0.08133591 0.4025836
## b_verb_type_training2.ct -0.08641562 0.1585316 -0.39282603 0.2232650

mcmc_plot(prod_attested_unattested_ent, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_ent))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00125
## 2   0.28725

# maximally vague priors for the intercept
prod_attested_unattested_ent_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent_final, variable = c("b_Intercept"))

##              Estimate  Est.Error       Q2.5     Q97.5
## b_Intercept 0.2400776 0.08259075 0.07946107 0.4006658

mcmc_plot(prod_attested_unattested_ent_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0.001916667

# c. we will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the 2 non-alternating verbs than for the novel verb (presumably the “unwitnessed” form has to be set arbitrarily here)


production_entrenchment_restricted_novel.df <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_restricted_novel.df<- subset(production_entrenchment_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_entrenchment_restricted_novel.df$attested_unattested)

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_entrenchment_restricted_novel.df$attested_unattested)

round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.571         0.549         0.516

# proportion of attested items for each verb type - we will now compare constr1/2 vs. novel


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.516 0.560

#what this means is that participants produce *unattested forms* less for the restricted (0.433) than they do for the novel (0.511)

production_entrenchment_restricted_novel.df$restricted_verb <- factor(production_entrenchment_restricted_novel.df$restricted_verb)
production_entrenchment_restricted_novel1.df = lizCenter(production_entrenchment_restricted_novel.df, list("restricted_verb"))


# maximally vague priors for the predictors and the intercept
prod_unattested_novel_ent_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|participant_private_id), data=production_entrenchment_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(prod_unattested_novel_ent_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                       Estimate  Est.Error        Q2.5     Q97.5
## b_Intercept          0.1822426 0.06794533  0.05096797 0.3168139
## b_restricted_verb.ct 0.1745369 0.13887380 -0.10107238 0.4479459

mcmc_plot(prod_unattested_novel_ent_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_unattested_novel_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0030000
## 2 0.1033333

Experiment 3

Load data

#Create the dataframes that we will be working on
combined_production_data.df <- read.csv("exp3_production_data.csv")

combined_judgment_data.df <- read.csv("exp3_judgment_data.csv")
combined_judgment_data.df$restricted_verb <- factor(combined_judgment_data.df$restricted_verb)
combined_judgment_data.df$condition <- factor(combined_judgment_data.df$condition)


#separately for entrenchment and preemption

#entrenchment
entrenchment_production.df <- subset(combined_production_data.df, condition == "entrenchment")
entrenchment_production.df$semantically_correct <- as.numeric(entrenchment_production.df$semantically_correct)
entrenchment_production.df$transitivity_test_scene2 <- factor(entrenchment_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in entrenchment
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

entrenchment_production.df$det1 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction1", 1, 0)
entrenchment_production.df$det2 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction2", 1, 0)
entrenchment_production.df$other <- ifelse(entrenchment_production.df$det_lenient_adapted == "other", 1, 0)
entrenchment_production.df$none <- ifelse(entrenchment_production.df$det_lenient_adapted == "none", 1, 0)


entrenchment_judgment.df <- subset(combined_judgment_data.df, condition == "entrenchment")
entrenchment_judgment.df$semantically_correct <- factor(entrenchment_judgment.df$semantically_correct)
entrenchment_judgment.df$transitivity_test_scene2 <- factor(entrenchment_judgment.df$transitivity_test_scene2)
entrenchment_judgment.df$restricted_verb <- factor(entrenchment_judgment.df$restricted_verb)


#preemption
preemption_production.df <- subset(combined_production_data.df, condition == "preemption")
preemption_production.df$semantically_correct <- as.numeric(preemption_production.df$semantically_correct)  #actually, all NAs here
preemption_production.df$transitivity_test_scene2 <- factor(preemption_production.df$transitivity_test_scene2)

# Create columns that we will need to run production analyses in pre-emption
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

preemption_production.df$det1 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction1", 1, 0)
preemption_production.df$det2 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction2", 1, 0)
preemption_production.df$other <- ifelse(preemption_production.df$det_lenient_adapted == "other", 1, 0)
preemption_production.df$none <- ifelse(preemption_production.df$det_lenient_adapted == "none", 1, 0)


preemption_judgment.df <- subset(combined_judgment_data.df, condition == "preemption")
preemption_judgment.df$semantically_correct <- factor(preemption_judgment.df$semantically_correct)
preemption_judgment.df$transitivity_test_scene2 <- factor(preemption_judgment.df$transitivity_test_scene2)
preemption_judgment.df$restricted_verb <- factor(preemption_judgment.df$restricted_verb)

Preregistered data analyses

Question 1: Have participants picked up on the difference in meaning between the two argument-structure constructions?

Production data

#Figure 13

RQ1_graph_productions.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph1 = aggregate(semantically_correct ~ verb_type_training2 + participant_private_id, RQ1_graph_productions.df, FUN=mean)

aggregated.graph1 <- rename(aggregated.graph1, verb = verb_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = correct  ~ verb,
                  data = aggregated.graph1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "% semantically correct",
                  cex.lab = 1,
                  cex.axis = 1,
                  cex.names = 1,
                  yaxt = "n")

axis(2, at = seq(0, 1, by = 0.25), las=1)
abline(h = 0.50, lty = 2)

#1 alternating verb production

alternating_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "alternating")

#and filter out responses where participants said something other than det1 or det2
alternating_prod.df = subset(alternating_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_alternating_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + participant_private_id, alternating_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_alternating_prod.df$semantically_correct),3)

## [1] 0.972

# average accuracy separately for causative and inchoative scenes
round(tapply(aggregated.means_alternating_prod.df$semantically_correct, aggregated.means_alternating_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.988         0.956

a = lizCenter(alternating_prod.df, list("transitivity_test_scene2"))   

# maximally vague priors for the intercept and the predictors
alternating_model <- brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|participant_private_id), data=a, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct" ))

##                                 Estimate Est.Error      Q2.5     Q97.5
## b_Intercept                    3.4843861 0.3395370  2.869377 4.1889350
## b_transitivity_test_scene2.ct -0.6580628 0.5485043 -1.752750 0.3920934

mcmc_plot(alternating_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.1121667

# no difference between construction 1 and construction 2

# Final model
# maximally vague priors for the intercept 
alternating_model_final = brm(formula = semantically_correct~1 + (1|participant_private_id), data=a, family = bernoulli(link = logit),set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error  Q2.5    Q97.5
## b_Intercept 3.352933 0.3172057 2.786 4.028725

mcmc_plot(alternating_model_final, variable = "b_Intercept", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

#2 novel verb production

novel_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & verb_type_training2 == "novel")

#and filter out responses where participants said something other than det1 or det2
novel_prod.df = subset(novel_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_novel_prod.df = aggregate(semantically_correct ~ transitivity_test_scene2 + participant_private_id, novel_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_novel_prod.df$semantically_correct),3)

## [1] 0.959

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_prod.df$semantically_correct, aggregated.means_novel_prod.df$transitivity_test_scene2, mean),3)

## construction1 construction2 
##         0.969         0.950

b = lizCenter(novel_prod.df, list("transitivity_test_scene2"))  


# maximally vague priors for the intercept and the predictors
novel_model <- brm(formula = semantically_correct~transitivity_test_scene2.ct + (1 + transitivity_test_scene2.ct|participant_private_id), data=b, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model, variable = c("b_Intercept", "b_transitivity_test_scene2.ct"))

##                                 Estimate Est.Error      Q2.5     Q97.5
## b_Intercept                    3.4395932 0.3795265  2.769943 4.2391092
## b_transitivity_test_scene2.ct -0.5921284 0.5890993 -1.797123 0.5219745

mcmc_plot(novel_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.8466667

# no difference between construction 1 and construction 2  
# Final model

# maximally vague priors for the intercept 
novel_model_final <- brm(formula = semantically_correct~1+ (1 + 1|participant_private_id), data=b, family = bernoulli(link = logit), set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.268721  0.372564 2.625329 4.079778

mcmc_plot(novel_model_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

Judgment data

#Figure 14
RQ1_graph_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating" |verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph2 = aggregate(response ~ verb_type_training2+ semantically_correct + participant_private_id, RQ1_graph_judgments.df, FUN=mean)
aggregated.graph2$semantically_correct <- recode(aggregated.graph2$semantically_correct, "1" = "yes","0" = "no")

aggregated.graph2 <- rename(aggregated.graph2, verb = verb_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = response ~ correct + verb,
                  data = aggregated.graph2,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

alternating_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "alternating")

#1 alternating verb judgments

# aggregated dataframe for means
aggregated.means_alternating_judgments = aggregate(response ~ transitivity_test_scene2 + semantically_correct + participant_private_id, alternating_judgments.df, FUN=mean)

# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_alternating_judgments$response, aggregated.means_alternating_judgments$semantically_correct, mean),3)

##     0     1 
## 2.350 4.862

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_alternating_judgments$response, list(aggregated.means_alternating_judgments$semantically_correct, aggregated.means_alternating_judgments$transitivity_test_scene2), mean),3)

##   construction1 construction2
## 0         2.288         2.413
## 1         4.888         4.838

c = lizCenter(alternating_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments <-brm(formula = response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|participant_private_id), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct","b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                          Estimate  Est.Error
## b_Intercept                                            3.61999537 0.06585723
## b_transitivity_test_scene2.ct                          0.03562877 0.07396541
## b_semantically_correct.ct                              2.47052282 0.13187948
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.17330322 0.14033210
##                                                             Q2.5     Q97.5
## b_Intercept                                            3.4921246 3.7491860
## b_transitivity_test_scene2.ct                         -0.1087127 0.1821042
## b_semantically_correct.ct                              2.2076735 2.7251085
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.4474400 0.1025763

samps = as.matrix(as.mcmc(alternating_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1         0.0000000
## 2         0.6845000
## 3         0.0000000
## 4         0.8964167

# no difference between construction 1 and construction 2

# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here)
alternating_model_judgments_final <-brm(formula = response~semantically_correct.ct + (1 + semantically_correct.ct|participant_private_id), data=c, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(alternating_model_judgments_final, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept               3.622345 0.06670889 3.491441 3.753979
## b_semantically_correct.ct 2.467649 0.13371877 2.200672 2.726265

mcmc_plot(alternating_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(alternating_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

#2 novel verb judgments

novel_judgments.df = subset(entrenchment_judgment.df, condition == "entrenchment" & verb_type_training2 == "novel")

# aggregated dataframe for means
aggregated.means_novel_judgments = aggregate(response ~ transitivity_test_scene2 + semantically_correct + participant_private_id, novel_judgments.df, FUN=mean)

# average accuracy for semantically correct vs. incorrect trials across causative and noncausative trial types
round(tapply(aggregated.means_novel_judgments$response, aggregated.means_novel_judgments$semantically_correct, mean),3)

##     0     1 
## 2.169 3.825

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_judgments$response, list(aggregated.means_novel_judgments$semantically_correct, aggregated.means_novel_judgments$transitivity_test_scene2), mean),3)

##   construction1 construction2
## 0         2.112         2.225
## 1         3.775         3.875

d = lizCenter(novel_judgments.df, list("transitivity_test_scene2", "semantically_correct"))  


# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments <-brm(formula = response~transitivity_test_scene2.ct * semantically_correct.ct + (1 + transitivity_test_scene2.ct*semantically_correct.ct|participant_private_id), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_judgments, variable = c("b_Intercept", "b_transitivity_test_scene2.ct","b_semantically_correct.ct", "b_transitivity_test_scene2.ct:semantically_correct.ct"))

##                                                           Estimate  Est.Error
## b_Intercept                                            2.975265540 0.14272054
## b_transitivity_test_scene2.ct                          0.106603607 0.07614699
## b_semantically_correct.ct                              1.591040410 0.20337653
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.006157815 0.16080378
##                                                              Q2.5     Q97.5
## b_Intercept                                            2.69626235 3.2564555
## b_transitivity_test_scene2.ct                         -0.04459483 0.2551572
## b_semantically_correct.ct                              1.18152180 1.9772712
## b_transitivity_test_scene2.ct:semantically_correct.ct -0.32542141 0.3082339

mcmc_plot(novel_model_judgments, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_transitivity_test_scene2.ct"] > 0)
C3=mean(samps[,"b_semantically_correct.ct"] < 0)
C4=mean(samps[,"b_transitivity_test_scene2.ct:semantically_correct.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1           0.00000
## 2           0.91925
## 3           0.00000
## 4           0.51550

# no difference between construction 1 and construction 2
# Final model

# maximally vague priors for the predictors (we don't interpret the intercept here) 
novel_model_judgments_final <-brm(formula = response~semantically_correct.ct + (1 + semantically_correct.ct|participant_private_id), data=d, family = gaussian(), set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(novel_model_judgments_final, variable = c("b_Intercept", "b_semantically_correct.ct"))

##                           Estimate Est.Error     Q2.5    Q97.5
## b_Intercept               2.975554 0.1375423 2.699915 3.244376
## b_semantically_correct.ct 1.592595 0.1955210 1.209105 1.969879

mcmc_plot(novel_model_judgments_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_judgments_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_semantically_correct.ct"] < 0)


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Question 2: Does statistical pre-emption constrain verb argument construction generalizations in adults (judgment data)?

#Figure 15

#first, filter our semantically incorrect trials

judgments_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel

judgments_novel.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

judgment_unattested_constr1.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
judgment_unattested_constr2.df <- subset(judgments_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

judgment_unattested_novel.df <- rbind(judgments_novel.df, judgment_unattested_constr1.df, judgment_unattested_constr2.df)

aggregated.means = aggregate(response ~ condition + restricted_verb + participant_private_id, judgment_unattested_novel.df, FUN=mean)
aggregated.means<- rename(aggregated.means, restricted = restricted_verb)

yarrr::pirateplot(formula = response ~ restricted + condition,
                  data = aggregated.means,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

judgment_unattested_novel_preemption.df <- subset(judgment_unattested_novel.df, condition == "preemption")   
round(tapply(judgment_unattested_novel_preemption.df$response, judgment_unattested_novel_preemption.df$restricted_verb, mean),3)

##    no   yes 
## 2.722 2.423

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_preemption_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(judgments_preemption_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate  Est.Error       Q2.5      Q97.5
## b_Intercept           2.5727318 0.07591153  2.4225884 2.72058222
## b_restricted_verb.ct -0.2825362 0.17159042 -0.6194975 0.05789206

mcmc_plot(judgments_preemption_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_preemption_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.95025

# BF analyses: we use the difference between attested and unattested in Experiment1 (SD = 0.65) as an estimate of the difference we expect in comparing rating for unattested vs. novel constructions 
Bf(0.18, 0.29, uniform = 0, meanoftheory = 0, sdtheory = 0.65, tail = 1)

## $LikelihoodTheory
## [1] 1.012346
## 
## $Likelihoodnull
## [1] 0.6053312
## 
## $BayesFactor
## [1] 1.672383

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.18, 0.29, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# RRs for which BF > 3
#ev_for_h1 <- subset(data.frame(range_test), BF > 3)
#low_threshold <- min(ev_for_h1$sdtheory)
#high_threshold <- max(ev_for_h1$sdtheory)
#print(low_threshold)
#print(high_threshold)

Question 3: Does statistical entrenchment constrain verb argument construction generalizations in adults (judgment data)?

#first, filter our semantically incorrect trials

entrenchment_judgment_unattested_novel.df <- subset(entrenchment_judgment.df, semantically_correct == "1")   

#we only want to keep novel

entrenchment_judgment_novel.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

entrenchment_judgment_unattested_constr1.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
entrenchment_judgment_unattested_constr2.df <- subset(entrenchment_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

entrenchment_judgment_unattested_novel.df <- rbind(entrenchment_judgment_novel.df, entrenchment_judgment_unattested_constr1.df, entrenchment_judgment_unattested_constr2.df)
entrenchment_judgment_unattested_novel.df$restricted_verb <- factor(entrenchment_judgment_unattested_novel.df$restricted_verb , levels = c("yes", "no"))

round(tapply(entrenchment_judgment_unattested_novel.df$response, entrenchment_judgment_unattested_novel.df$restricted_verb, mean),3)

##   yes    no 
## 4.425 3.825

#Center variables of interest using the lizCenter function:
d_unattested_novel_entrenchment = lizCenter(entrenchment_judgment_unattested_novel.df, list("restricted_verb"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_entrenchment_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct, data=d_unattested_novel_entrenchment, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_entrenchment_model, variable = c("b_Intercept", "b_restricted_verb.ct"))

##                        Estimate Est.Error       Q2.5      Q97.5
## b_Intercept           4.1339581 0.1460519  3.8464818  4.4246023
## b_restricted_verb.ct -0.5747833 0.1891309 -0.9555493 -0.2039924

mcmc_plot(judgments_entrenchment_model, variable = "b_Intercept", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_entrenchment_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.9983333

# drawing a max based on the difference between attested vs. unattested in this experiment (this was sig. evidence for entrenchment)
Bf(0.19, -0.58, uniform = 0, meanoftheory = 0, sdtheory = 0.38/2, tail = 1)

## $LikelihoodTheory
## [1] 0.004503071
## 
## $Likelihoodnull
## [1] 0.01989102
## 
## $BayesFactor
## [1] 0.2263872

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.19, -0.58, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF < 3
ev_for_h1 <- subset(data.frame(range_test), BF < 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0

print(high_threshold)

## [1] 4

Question 4: Is the effect of statistical pre-emption larger than entrenchment (judgment data)?

#first, filter our semantically incorrect trials

all_judgment_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "1")   

#we only want to keep novel

all_judgment_novel.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "novel")   

#and restricted items

all_judgment_unattested_constr1.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction1" & attested_unattested == "0")   
all_judgment_unattested_constr2.df <- subset(all_judgment_unattested_novel.df, verb_type_training2 == "construction2" & attested_unattested == "0")   

all_judgment_unattested_novel.df <- rbind(all_judgment_novel.df, all_judgment_unattested_constr1.df, all_judgment_unattested_constr2.df)

round(tapply(all_judgment_unattested_novel.df$response, list(all_judgment_unattested_novel.df$restricted_verb, all_judgment_unattested_novel.df$condition), mean),3)

##     entrenchment preemption
## no         3.825      2.722
## yes        4.425      2.423

# preemption worked and opposite effect for entrenchment

#Center variables of interest using the lizCenter function:
df = lizCenter(all_judgment_unattested_novel.df, list("restricted_verb", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
judgments_pre_vs_ent_model <- brm(formula = response~(1 +restricted_verb.ct|participant_private_id)+restricted_verb.ct * condition.ct, data=df, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_pre_vs_ent_model, variable = c("b_Intercept", "b_restricted_verb.ct","b_condition.ct", "b_restricted_verb.ct:condition.ct"))

##                                     Estimate  Est.Error       Q2.5      Q97.5
## b_Intercept                        3.0271605 0.07853672  2.8742014  3.1814381
## b_restricted_verb.ct              -0.0384540 0.13266064 -0.2976345  0.2224343
## b_condition.ct                    -1.5223832 0.14999047 -1.8098310 -1.2241765
## b_restricted_verb.ct:condition.ct -0.8682916 0.25253906 -1.3632946 -0.3705756

mcmc_plot(judgments_pre_vs_ent_model, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(judgments_pre_vs_ent_model))

C1=mean(samps[,"b_Intercept"] < 0) 
C2=mean(samps[,"b_restricted_verb.ct"] > 0)
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_restricted_verb.ct:condition.ct"] > 0) 


pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1      0.0000000000
## 2      0.3829166667
## 3      0.0000000000
## 4      0.0001666667

#roughly predicted effect size from previous study was 1.0. Use it as an estimate of the effect we expect here
Bf(0.25, 0.85, uniform = 0, meanoftheory = 0, sdtheory = 1.00, tail = 1)

## $LikelihoodTheory
## [1] 0.5506764
## 
## $Likelihoodnull
## [1] 0.004928877
## 
## $BayesFactor
## [1] 111.7245

H1RANGE = seq(0,4,by=0.01) # [5-1]-[0] - max effect of preemption minus no effect of entrenchment
range_test <- Bf_range(0.25, 0.85, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.09

print(high_threshold)

## [1] 4

Exploratory data analyses

Effect of statistical pre-emption: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

# Figure 16
judgments_unattested_attested.df <- subset(combined_judgment_data.df, semantically_correct == "1")   
judgments_unattested_attested.df <- subset(judgments_unattested_attested.df, restricted_verb == "yes")   


aggregated.means1 = aggregate(response ~ condition + attested_unattested + participant_private_id, judgments_unattested_attested.df , FUN=mean)
aggregated.means1<- rename(aggregated.means1, attested = attested_unattested)

aggregated.means1$attested<- recode(aggregated.means1$attested, "1" = "yes","0" = "no")


yarrr::pirateplot(formula = response ~  attested  + condition,
                  data = aggregated.means1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

# analyses
attested_vs_unattested = subset(preemption_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested$response, attested_vs_unattested$attested_unattested, mean),3)

##     0     1 
## 2.423 4.908

#Center variables of interest using the lizCenter function:
d_attested_unattested = lizCenter(attested_vs_unattested , list("attested_unattested"))


# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_preemption <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(attested_unattested_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept              3.679930 0.06111709 3.561695 3.801650
## b_attested_unattested.ct 2.453843 0.11126771 2.234294 2.669445

mcmc_plot(attested_unattested_preemption, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

# this prioris drawn from Experiment 1
Bf(0.11, 2.45, uniform = 0, meanoftheory = 0, sdtheory = 2.55, tail = 1)

## $LikelihoodTheory
## [1] 0.1972052
## 
## $Likelihoodnull
## [1] 6.891828e-108
## 
## $BayesFactor
## [1] 2.861436e+106

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.11, 2.45, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.01

print(high_threshold)

## [1] 4

Effect of statistical entrenchment: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

attested_vs_unattested_ent = subset(entrenchment_judgment.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_ent$response, attested_vs_unattested_ent$attested_unattested, mean),3)

##     0     1 
## 4.425 4.950

#Center variables of interest using the lizCenter function:
d_attested_unattested_ent = lizCenter(attested_vs_unattested_ent, list("attested_unattested"))

attested_unattested_entrenchment <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested_ent, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate  Est.Error      Q2.5     Q97.5
## b_Intercept              4.6920092 0.06668343 4.5577081 4.8226758
## b_attested_unattested.ct 0.5123686 0.13193648 0.2563227 0.7733299

mcmc_plot(attested_unattested_entrenchment, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(attested_unattested_entrenchment))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##      c(C1, C2)
## 1 0.0000000000
## 2 0.0001666667

# expect a difference of 0.38 from previous work
Bf(0.13, 0.51, uniform = 0, meanoftheory = 0, sdtheory = 0.38, tail = 1)

## $LikelihoodTheory
## [1] 0.887002
## 
## $Likelihoodnull
## [1] 0.001396224
## 
## $BayesFactor
## [1] 635.2864

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.13, 0.51, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.04

print(high_threshold)

## [1] 4

Entrenchment vs. preemption: ratings for witnessed vs. unwitnessed forms

attested_vs_unattested_across = subset(combined_judgment_data.df, restricted_verb == "yes" & semantically_correct == "1")

round(tapply(attested_vs_unattested_across$response, list(attested_vs_unattested_across$condition, attested_vs_unattested_across$attested_unattested), mean),3)

##                  0     1
## entrenchment 4.425 4.950
## preemption   2.423 4.908

#Center variables of interest using the lizCenter function:
df_attested_unattested = lizCenter(attested_vs_unattested_across, list("attested_unattested", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment_preemption <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct * condition.ct, data=df_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment_preemption, variable = c("b_Intercept", "b_attested_unattested.ct",  "b_condition.ct", "b_attested_unattested.ct:condition.ct"))

##                                        Estimate  Est.Error      Q2.5      Q97.5
## b_Intercept                            3.970834 0.04530757  3.883186  4.0611365
## b_attested_unattested.ct               1.898064 0.08475568  1.729359  2.0620102
## b_condition.ct                        -0.990629 0.08904309 -1.165525 -0.8136171
## b_attested_unattested.ct:condition.ct  1.890947 0.16686121  1.557075  2.2169299

mcmc_plot(attested_unattested_entrenchment_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_attested_unattested.ct:condition.ct"] < 0) 

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1                 0
## 2                 0
## 3                 0
## 4                 0

#roughly predicted effect size from previous study 2.11

Bf(0.17, 1.89, uniform = 0, meanoftheory = 0, sdtheory = 2.11, tail = 1)

## $LikelihoodTheory
## [1] 0.2530173
## 
## $Likelihoodnull
## [1] 3.393215e-27
## 
## $BayesFactor
## [1] 7.456566e+25

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.17, 1.89, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.02

print(high_threshold)

## [1] 4

Production data: Effect of statistical pre-emption

#Figure 17

data_long <- gather(preemption_production.df, det_type, produced, det1:none, factor_key=TRUE)

p = ggplot(data_long, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme_bw()+
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "periphrastic causative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "periphrastic causative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_verb"))


#a. Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted verbs against chance 

production_preemption_attested_unattested.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, restricted_verb =="yes")

round(tapply(production_preemption_attested_unattested.df$attested_unattested, production_preemption_attested_unattested.df$verb_type_training2, mean),3)

## construction1 construction2 
##         0.992         0.997

production_preemption_attested_unattested.df$verb_type_training2 <- factor(production_preemption_attested_unattested.df$verb_type_training2)

df_prod = lizCenter(production_preemption_attested_unattested.df , list("verb_type_training2"))  

# maximally vague priors for the predictors and the intercept
prod_attested_unattested = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|participant_private_id), data=df_prod, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested, variable = c("b_Intercept","b_verb_type_training2.ct"))

##                           Estimate Est.Error      Q2.5    Q97.5
## b_Intercept              4.7657488 0.3724566  4.102279 5.560532
## b_verb_type_training2.ct 0.3238239 0.6041177 -0.828422 1.507577

mcmc_plot(prod_attested_unattested, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.70025

#same analyses without verb_training_type

# maximally vague priors for the intercept
prod_attested_unattested_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 4.667878 0.3682112 4.006309 5.467753

mcmc_plot(prod_attested_unattested_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

# We will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the restricted verbs than for the novel verb

production_preemption_restricted_novel.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_restricted_novel.df<- subset(production_preemption_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_preemption_restricted_novel.df$attested_unattested)

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$verb_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_preemption_restricted_novel.df$attested_unattested)

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.992         0.997         0.565

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.565 0.995

production_preemption_restricted_novel.df$restricted_verb <- factor(production_preemption_restricted_novel.df$restricted_verb)
production_preemption_restricted_novel1.df = lizCenter(production_preemption_restricted_novel.df, list("restricted_verb"))

# maximally vague priors for the predictors and the intercept
prod_unattested_novel_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|participant_private_id), data=production_preemption_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                      Estimate Est.Error     Q2.5    Q97.5
## b_Intercept          3.452220 0.2774725 2.940953 4.029338
## b_restricted_verb.ct 3.496777 0.5378190 2.402721 4.516781

mcmc_plot(prod_unattested_novel_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_unattested_novel_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

Production data: Effect of statistical entrenchment

#Figure 12

# here, we want to see how often participants say det1 (e.g. transitive-only) for a det2 (intransitive-only) verb in the intransitive condition at test. 
# in other words, do they produce the unattested in a semantically correct trial?  

data_long_e <- gather(entrenchment_production.df, det_type, produced, det1:none, factor_key=TRUE)

data_long_e$transitivity_test_scene2 <-recode(data_long_e$transitivity_test_scene2, "construction1" = "test: transitive causative","construction2" = "test: intransitive inchoative")


p = ggplot(data_long_e, aes(x = verb_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme(panel.grid.major = element_blank()) +
  facet_grid("transitivity_test_scene2") +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("transitive causative", "intransitive inchoative", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "transitive causative", "construction2" = "intransitive inchoative", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#a. Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted verbs against chance 

production_entrenchment_attested_unattested.df  <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, restricted_verb =="yes")

#How much of the time are participants producing attested items?
round(mean(production_entrenchment_attested_unattested.df$attested_unattested),3)

## [1] 0.572

# and separately for each verb type
round(tapply(production_entrenchment_attested_unattested.df$attested_unattested, production_entrenchment_attested_unattested.df$verb_type_training2, mean),3)

## construction1 construction2 
##         0.588         0.556

production_entrenchment_attested_unattested.df$verb_type_training2 <- factor(production_entrenchment_attested_unattested.df$verb_type_training2)
df_prod_ent = lizCenter((production_entrenchment_attested_unattested.df), list("verb_type_training2"))  


# maximally vague priors for the predictors and the intercept
prod_attested_unattested_ent = brm(formula = attested_unattested ~verb_type_training2.ct + (1 + verb_type_training2.ct|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent, variable = c("b_Intercept", "b_verb_type_training2.ct"))

##                            Estimate  Est.Error       Q2.5     Q97.5
## b_Intercept               0.2949456 0.09496386  0.1108591 0.4827951
## b_verb_type_training2.ct -0.1281749 0.16161117 -0.4451122 0.1847387

mcmc_plot(prod_attested_unattested_ent, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_ent))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_verb_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##     c(C1, C2)
## 1 0.001666667
## 2 0.219500000

# maximally vague priors for the intercept
prod_attested_unattested_ent_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent_final, variable = c("b_Intercept"))

##              Estimate  Est.Error      Q2.5     Q97.5
## b_Intercept 0.2948758 0.09324019 0.1152931 0.4782878

mcmc_plot(prod_attested_unattested_ent_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0.001

# c. we will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the 2 non-alternating verbs than for the novel verb (presumably the “unwitnessed” form has to be set arbitrarily here)


production_entrenchment_restricted_novel.df <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_restricted_novel.df<- subset(production_entrenchment_restricted_novel.df, verb_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_entrenchment_restricted_novel.df$attested_unattested)

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$verb_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_entrenchment_restricted_novel.df$attested_unattested)

round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$verb_type_training2, mean),3)

## construction1 construction2         novel 
##         0.588         0.556         0.511

# proportion of attested items for each verb type - we will now compare constr1/2 vs. novel


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$restricted_verb, mean),3)

##    no   yes 
## 0.511 0.572

#what this means is that participants produce *unattested forms* less for the restricted than they do for the novel

production_entrenchment_restricted_novel.df$restricted_verb <- factor(production_entrenchment_restricted_novel.df$restricted_verb)
production_entrenchment_restricted_novel1.df = lizCenter(production_entrenchment_restricted_novel.df, list("restricted_verb"))


# maximally vague priors for the predictors and the intercept
prod_unattested_novel_ent_final = brm(formula = attested_unattested ~restricted_verb.ct + (1 + restricted_verb.ct|participant_private_id), data=production_entrenchment_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(prod_unattested_novel_ent_final, variable = c("b_Intercept","b_restricted_verb.ct"))

##                       Estimate  Est.Error        Q2.5     Q97.5
## b_Intercept          0.2103495 0.07072125  0.07361361 0.3470228
## b_restricted_verb.ct 0.2432940 0.14806638 -0.04346983 0.5390849

mcmc_plot(prod_unattested_novel_ent_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_unattested_novel_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_verb.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##     c(C1, C2)
## 1 0.001833333
## 2 0.047416667

Experiment 4

Load data

#Create the dataframes that we will be working on
combined_production_data.df <- read.csv("exp4_production_data.csv")

combined_judgment_data.df <- read.csv("exp4_judgment_data.csv")
combined_judgment_data.df$restricted_noun <- factor(combined_judgment_data.df$restricted_noun)
combined_judgment_data.df$condition <- factor(combined_judgment_data.df$condition)


#separately for entrenchment and preemption

#entrenchment
entrenchment_production.df <- subset(combined_production_data.df, condition == "entrenchment")
entrenchment_production.df$semantically_correct <- as.numeric(entrenchment_production.df$semantically_correct)
entrenchment_production.df$noun_type_test <- factor(entrenchment_production.df$noun_type_test)

# Create columns that we will need to run production analyses in entrenchment
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

entrenchment_production.df$det1 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction1", 1, 0)
entrenchment_production.df$det2 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction2", 1, 0)
entrenchment_production.df$other <- ifelse(entrenchment_production.df$det_lenient_adapted == "other", 1, 0)
entrenchment_production.df$none <- ifelse(entrenchment_production.df$det_lenient_adapted == "none", 1, 0)


entrenchment_judgment.df <- subset(combined_judgment_data.df, condition == "entrenchment")
entrenchment_judgment.df$semantically_correct <- factor(entrenchment_judgment.df$semantically_correct)
entrenchment_judgment.df$noun_type_test <- factor(entrenchment_judgment.df$noun_type_test)
entrenchment_judgment.df$restricted_noun <- factor(entrenchment_judgment.df$restricted_noun)


#preemption
preemption_production.df <- subset(combined_production_data.df, condition == "preemption")
preemption_production.df$semantically_correct <- as.numeric(preemption_production.df$semantically_correct)  #actually, all NAs here
preemption_production.df$noun_type_test <- factor(preemption_production.df$noun_type_test)

# Create columns that we will need to run production analyses in pre-emption
# We want some columns coding which of particle 1, particle 2, 'other' and 'none' was produced

preemption_production.df$det1 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction1", 1, 0)
preemption_production.df$det2 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction2", 1, 0)
preemption_production.df$other <- ifelse(preemption_production.df$det_lenient_adapted == "other", 1, 0)
preemption_production.df$none <- ifelse(preemption_production.df$det_lenient_adapted == "none", 1, 0)


preemption_judgment.df <- subset(combined_judgment_data.df, condition == "preemption")
preemption_judgment.df$semantically_correct <- factor(preemption_judgment.df$semantically_correct)
preemption_judgment.df$noun_type_test <- factor(preemption_judgment.df$noun_type_test)
preemption_judgment.df$restricted_noun <- factor(preemption_judgment.df$restricted_noun)

Preregistered data analyses

Question 1: Have participants picked up on the difference in meaning between singular/plural marking?

Production data

#Figure X
RQ1_graph_productions.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "alternating" |noun_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph1 = aggregate(semantically_correct ~ noun_type_training2 + participant_private_id, RQ1_graph_productions.df, FUN=mean)

aggregated.graph1 <- rename(aggregated.graph1, noun = noun_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = correct  ~ noun,
                  data = aggregated.graph1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "% semantically correct",
                  cex.lab = 1,
                  cex.axis = 1,
                  cex.names = 1,
                  yaxt = "n")

axis(2, at = seq(0, 1, by = 0.25), las=1)
abline(h = 0.50, lty = 2)

#1 alternating noun production

alternating_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "alternating")
#and filter out responses where participants said something other than det1 or det2
alternating_prod.df = subset(alternating_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_alternating_prod.df = aggregate(semantically_correct ~ noun_type_test + participant_private_id, alternating_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_alternating_prod.df$semantically_correct),3)

## [1] 0.977

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_alternating_prod.df$semantically_correct, aggregated.means_alternating_prod.df$ noun_type_test, mean),3)

## construction1 construction2 
##         0.994         0.959

a = lizCenter(alternating_prod.df, list("noun_type_test"))  


# maximally vague priors for the intercept and the predictors
alternating_model <-brm(formula = semantically_correct~noun_type_test.ct + (1 + noun_type_test.ct|participant_private_id), data=a, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model, variable = c("b_Intercept", "b_noun_type_test.ct"))

##                       Estimate Est.Error      Q2.5     Q97.5
## b_Intercept          4.0096271 0.4599040  3.196430 4.9878499
## b_noun_type_test.ct -0.4515319 0.6749788 -1.749938 0.9003246

mcmc_plot(alternating_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(alternating_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_test.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.7560833

# no difference between construction 1 and construction 2  
# Final model

# maximally vague priors for the intercept 
alternating_model_final = brm(formula = semantically_correct~1 + (1|participant_private_id), data=a, family = bernoulli(link = logit),set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error    Q2.5    Q97.5
## b_Intercept 3.708707 0.4318145 2.96409 4.647277

mcmc_plot(alternating_model_final, variable = "b_Intercept", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(alternating_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

#2 novel noun production

novel_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "novel")
#and filter out responses where participants said something other than det1 or det2
novel_prod.df = subset(novel_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_novel_prod.df = aggregate(semantically_correct ~ noun_type_test + participant_private_id, novel_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_novel_prod.df$semantically_correct),3)

## [1] 0.939

# average accuracy separately for causative and noncausative scenes
round(tapply(aggregated.means_novel_prod.df$semantically_correct, aggregated.means_novel_prod.df$noun_type_test, mean),3)

## construction1 construction2 
##         0.954         0.924

b = lizCenter(novel_prod.df, list("noun_type_test"))  

# maximally vague priors for the intercept and the predictors
novel_model <- brm(formula = semantically_correct~noun_type_test.ct + (1 + noun_type_test.ct|participant_private_id), data=b, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model, variable = c("b_Intercept", "b_noun_type_test.ct"))

##                       Estimate Est.Error      Q2.5     Q97.5
## b_Intercept          3.6299602 0.5136117  2.650576 4.6771417
## b_noun_type_test.ct -0.3836952 0.6075892 -1.538038 0.8298317

mcmc_plot(novel_model, variable = "^b_", regex = TRUE)


samps = as.matrix(as.mcmc(novel_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_test.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.7394167

# no difference between construction 1 and construction 2  
# Final model   

# maximally vague priors for the intercept 
novel_model_final <- brm(formula = semantically_correct~1+ (1|participant_private_id), data=b, family = bernoulli(link = logit), set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.559027 0.4955513 2.630708 4.571067

mcmc_plot(novel_model_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

Judgment data

There are no semantically incorrect trials in Experiment 4, thus, these analyses are not possible.

Question 2: Does statistical preemption constrain morphological generalizations in adults (judgment data)?

#Figure X

#first, filter our semantically incorrect trials

judgments_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "yes")   

#we only want to keep novel

judgments_novel.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "novel")   

#and restricted items

judgment_unattested_constr1.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "construction1" & attested_unattested == "0")   
judgment_unattested_constr2.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

judgment_unattested_novel.df <- rbind(judgments_novel.df, judgment_unattested_constr1.df, judgment_unattested_constr2.df)

aggregated.means = aggregate(response ~ condition + restricted_noun + participant_private_id, judgment_unattested_novel.df, FUN=mean)
aggregated.means<- rename(aggregated.means, restricted = restricted_noun)

yarrr::pirateplot(formula = response ~ restricted + condition,
                  data = aggregated.means,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

judgment_unattested_novel_preemption.df <- subset(judgment_unattested_novel.df, condition == "preemption")   
round(tapply(judgment_unattested_novel_preemption.df$response, judgment_unattested_novel_preemption.df$restricted_noun, mean),3)

##    no   yes 
## 3.475 2.717

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_noun"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_preemption_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_preemption_model, variable = c("b_Intercept","b_restricted_noun.ct"))

##                        Estimate Est.Error      Q2.5      Q97.5
## b_Intercept           3.1047592 0.1200532  2.871136  3.3451980
## b_restricted_noun.ct -0.7237057 0.2145152 -1.145436 -0.3040576

mcmc_plot(judgments_preemption_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_preemption_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1    0.0000
## 2    0.9995

# BF analyses: we use the difference between attested and novel in Experiment 1 (SD = 0.65) as an estimate of the difference we expect here

Bf(0.22, 0.72, uniform = 0, meanoftheory = 0, sdtheory = 0.65, tail = 1)

## $LikelihoodTheory
## [1] 0.6698739
## 
## $Likelihoodnull
## [1] 0.008564045
## 
## $BayesFactor
## [1] 78.21934

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.22, 0.72, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# RRs for which BF > 3
ev_for_h1 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.08

print(high_threshold)

## [1] 4

Question 3: Does statistical entrenchment constrain morphological generalizations in adults (judgment data)?

#no semantically incorrect trials here

#we only want to keep novel

entrenchment_judgment_novel.df <- subset(entrenchment_judgment.df, noun_type_training2 == "novel")   

#and restricted items

entrenchment_judgment_unattested_constr1.df <- subset(entrenchment_judgment.df, noun_type_training2 == "construction1" & attested_unattested == "0")

entrenchment_judgment_unattested_constr2.df <- subset(entrenchment_judgment.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

#bind new dataframe
entrenchment_judgment_unattested_novel.df <- rbind(entrenchment_judgment_novel.df, entrenchment_judgment_unattested_constr1.df, entrenchment_judgment_unattested_constr2.df)
entrenchment_judgment_unattested_novel.df$restricted_noun <- factor(entrenchment_judgment_unattested_novel.df$restricted_noun , levels = c("yes", "no"))

round(tapply(entrenchment_judgment_unattested_novel.df$response, entrenchment_judgment_unattested_novel.df$restricted_noun, mean),3)

##   yes    no 
## 4.537 4.646

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(entrenchment_judgment_unattested_novel.df, list("restricted_noun"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_entrenchment_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(judgments_entrenchment_model, variable = c("b_Intercept","b_restricted_noun.ct"))

##                       Estimate  Est.Error       Q2.5     Q97.5
## b_Intercept          4.5919269 0.09243736  4.4121429 4.7760506
## b_restricted_noun.ct 0.1070878 0.14117688 -0.1732026 0.3819051

mcmc_plot(judgments_entrenchment_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_entrenchment_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1   0.00000
## 2   0.22475

# use unattested vs. novel in original study as an estimate of difference in unattested vs. novel
Bf(0.14, 0.11, uniform = 0, meanoftheory = 0, sdtheory = 0.38/2, tail = 1)

## $LikelihoodTheory
## [1] 2.237763
## 
## $Likelihoodnull
## [1] 2.092796
## 
## $BayesFactor
## [1] 1.06927

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.14, 0.11, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)


# find values for which BF is inconclusive 
ev_for_h1 <- subset(data.frame(range_test), BF < 3 & BF > 1/3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.01

print(high_threshold)

## [1] 0.88

# find out how many more participants we would need for conclusive evidence for entrenchment (BF > 3)
invisible(Bf_powercalc(0.14,  0.11, uniform=0, meanoftheory=0, sdtheory=0.38/2, tail=1, N=40, min=30, max=400))

#N = 238

Question 4: Is the effect of statistical pre-emption larger than entrenchment (judgment data)?

#all are semantically correct trials

#we only want to keep novel

all_judgment_novel.df <- subset(combined_judgment_data.df, noun_type_training2 == "novel")   

#and restricted items

all_judgment_unattested_constr1.df <- subset(combined_judgment_data.df, noun_type_training2 == "construction1" & attested_unattested == "0")   
all_judgment_unattested_constr2.df <- subset(combined_judgment_data.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

all_judgment_unattested_novel.df <- rbind(all_judgment_novel.df, all_judgment_unattested_constr1.df, all_judgment_unattested_constr2.df)
all_judgment_unattested_novel.df$restricted_noun <- factor(all_judgment_unattested_novel.df$restricted_noun , levels = c("yes", "no"))

round(tapply(all_judgment_unattested_novel.df$response, list(all_judgment_unattested_novel.df$restricted_noun, all_judgment_unattested_novel.df$condition), mean),3)

##     entrenchment preemption
## yes        4.537      2.717
## no         4.646      3.475

round(tapply(all_judgment_unattested_novel.df$response, all_judgment_unattested_novel.df$condition, mean),3)

## entrenchment   preemption 
##        4.592        3.096

#Center variables of interest using the lizCenter function:
df = lizCenter(all_judgment_unattested_novel.df, list("restricted_noun", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
judgments_pre_vs_ent_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct * condition.ct, data=df, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_pre_vs_ent_model, variable = c("b_Intercept", "b_restricted_noun.ct","b_condition.ct", "b_restricted_noun.ct:condition.ct"))

##                                     Estimate  Est.Error        Q2.5      Q97.5
## b_Intercept                        3.8472420 0.07696599  3.69505159  3.9985573
## b_restricted_noun.ct               0.4260149 0.13109895  0.17015045  0.6872192
## b_condition.ct                    -1.4517177 0.14792109 -1.74229562 -1.1628732
## b_restricted_noun.ct:condition.ct  0.5878361 0.25376077  0.09281596  1.0863075

mcmc_plot(judgments_pre_vs_ent_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_pre_vs_ent_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)
C3=mean(samps[,"b_condition.ct"] > 0)
C4=mean(samps[,"b_restricted_noun.ct:condition.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1      0.0000000000
## 2      0.0008333333
## 3      0.0000000000
## 4      0.0101666667

#roughly predicted effect size from previous study was 1.0. Use it as an estimate of the effect we expect here
Bf(0.25, 0.57, uniform = 0, meanoftheory = 0, sdtheory = 1.00, tail = 1)

## $LikelihoodTheory
## [1] 0.6551184
## 
## $Likelihoodnull
## [1] 0.1186183
## 
## $BayesFactor
## [1] 5.52291

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.25, 0.57, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.15

print(high_threshold)

## [1] 2.12

Exploratory data analyses

Exploratory data analyses

Effect of statistical pre-emption: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

# Figure x
judgments_unattested_attested.df <- subset(combined_judgment_data.df, semantically_correct == "yes")   
judgments_unattested_attested.df <- subset(judgments_unattested_attested.df, restricted_noun == "yes")   


aggregated.means1 = aggregate(response ~ condition + attested_unattested + participant_private_id, judgments_unattested_attested.df , FUN=mean)
aggregated.means1<- rename(aggregated.means1, attested = attested_unattested)

aggregated.means1$attested<- recode(aggregated.means1$attested, "1" = "yes","0" = "no")


yarrr::pirateplot(formula = response ~  attested  + condition,
                  data = aggregated.means1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

# analyses
attested_vs_unattested = subset(preemption_judgment.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested$response, attested_vs_unattested$attested_unattested, mean),3)

##     0     1 
## 2.717 4.375

#Center variables of interest using the lizCenter function:
d_attested_unattested = lizCenter(attested_vs_unattested , list("attested_unattested"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_preemption <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate  Est.Error      Q2.5    Q97.5
## b_Intercept              3.555419 0.07734695 3.4042914 3.707816
## b_attested_unattested.ct 1.488844 0.32121247 0.8546146 2.111764

mcmc_plot(attested_unattested_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

# prior from previous study with adults: 2.55
Bf(0.33, 1.49, uniform = 0, meanoftheory = 0, sdtheory = 2.55, tail = 1)

## $LikelihoodTheory
## [1] 0.2623458
## 
## $Likelihoodnull
## [1] 4.523803e-05
## 
## $BayesFactor
## [1] 5799.231

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.33, 1.49, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.09

print(high_threshold)

## [1] 4

Effect of statistical entrenchment: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

attested_vs_unattested_ent = subset(entrenchment_judgment.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested_ent$response, attested_vs_unattested_ent$attested_unattested, mean),3)

##     0     1 
## 4.537 4.850

#Center variables of interest using the lizCenter function:
d_attested_unattested_ent = lizCenter(attested_vs_unattested_ent, list("attested_unattested"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested_ent, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate  Est.Error       Q2.5     Q97.5
## b_Intercept              4.6962774 0.08827469 4.52136338 4.8680174
## b_attested_unattested.ct 0.3063458 0.14363335 0.02140561 0.5858677

mcmc_plot(attested_unattested_entrenchment, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##    c(C1, C2)
## 1 0.00000000
## 2 0.01808333

# expect a difference of 0.38 from previous work
Bf(0.14, 0.31, uniform = 0, meanoftheory = 0, sdtheory = 0.38, tail = 1)

## $LikelihoodTheory
## [1] 1.442615
## 
## $Likelihoodnull
## [1] 0.2455251
## 
## $BayesFactor
## [1] 5.875631

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.14, 0.31, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.09

print(high_threshold)

## [1] 1.01

Entrenchment vs. preemption: ratings for witnessed vs. unwitnessed forms

attested_vs_unattested_across = subset(combined_judgment_data.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested_across$response, list(attested_vs_unattested_across$condition, attested_vs_unattested_across$attested_unattested), mean),3)

##                  0     1
## entrenchment 4.537 4.850
## preemption   2.717 4.375

#Center variables of interest using the lizCenter function:
df_attested_unattested = lizCenter(attested_vs_unattested_across, list("attested_unattested", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment_preemption <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct * condition.ct, data=df_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate  Est.Error      Q2.5    Q97.5
## b_Intercept              4.1241709 0.05724223 4.0115279 4.236412
## b_attested_unattested.ct 0.9514644 0.18102707 0.5988433 1.300932

mcmc_plot(attested_unattested_entrenchment_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_attested_unattested.ct:condition.ct"] < 0) 

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1      0.0000000000
## 2      0.0000000000
## 3      0.0000000000
## 4      0.0003333333

#roughly predicted effect size from previous study 2.11
Bf(0.34, 1.18, uniform = 0, meanoftheory = 0, sdtheory = 2.11, tail = 1)

## $LikelihoodTheory
## [1] 0.3204584
## 
## $Likelihoodnull
## [1] 0.002843783
## 
## $BayesFactor
## [1] 112.6874

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.34, 1.18, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.12

print(high_threshold)

## [1] 4

Production data: Effect of statistical pre-emption

data_long <- gather(preemption_production.df, det_type, produced, det1:none, factor_key=TRUE)

p = ggplot(data_long, aes(x = noun_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme_bw()+
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("construction1 ", "construction2", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "construction 1", "construction2" = "construction2", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted nouns against chance 

production_preemption_attested_unattested.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, restricted_noun =="yes")

round(tapply(production_preemption_attested_unattested.df $attested_unattested, production_preemption_attested_unattested.df $noun_type_training2, mean),3)

## construction1 construction2 
##         0.915         0.823

production_preemption_attested_unattested.df$noun_type_training2 <- factor(production_preemption_attested_unattested.df$noun_type_training2)

df_prod = lizCenter(production_preemption_attested_unattested.df , list("noun_type_training2"))  

# maximally vague priors for the predictors and the intercept
prod_attested_unattested = brm(formula = attested_unattested ~noun_type_training2.ct + (1 + noun_type_training2.ct|participant_private_id), data=df_prod, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested, variable = c("b_Intercept","b_noun_type_training2.ct"))

##                            Estimate Est.Error      Q2.5     Q97.5
## b_Intercept               3.0753259 0.6690992  1.709709 4.3204241
## b_noun_type_training2.ct -0.7684097 0.6016029 -1.901445 0.4503228

mcmc_plot(prod_attested_unattested, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##      c(C1, C2)
## 1 8.333333e-05
## 2 1.049167e-01

# maximally vague priors for the intercept
prod_attested_unattested_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 2.977183  0.626751 1.693474 4.173563

mcmc_plot(prod_attested_unattested_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

# We will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the restricted verbs than for the novel verb

production_preemption_restricted_novel.df <- subset(preemption_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_restricted_novel.df<- subset(production_preemption_restricted_novel.df, noun_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$noun_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_preemption_restricted_novel.df$attested_unattested)

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$noun_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_preemption_restricted_novel.df$attested_unattested)

round(tapply(production_preemption_restricted_novel.df$attested_unattested, production_preemption_restricted_novel.df$noun_type_training2, mean),3)

## construction1 construction2         novel 
##         0.915         0.823         0.497

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$restricted_noun, mean),3)

##    no   yes 
## 0.497 0.869

production_preemption_restricted_novel.df$restricted_noun <- factor(production_preemption_restricted_novel.df$restricted_noun)
production_preemption_restricted_novel1.df = lizCenter(production_preemption_restricted_novel.df, list("restricted_noun"))

# maximally vague priors for the predictors and the intercept
prod_unattested_novel_final = brm(formula = attested_unattested ~restricted_noun.ct + (1 + restricted_noun.ct|participant_private_id), data=production_preemption_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(prod_unattested_novel_final, variable = c("b_Intercept","b_restricted_noun.ct"))

##                      Estimate Est.Error      Q2.5    Q97.5
## b_Intercept          2.054763 0.4324699 1.1813734 2.890399
## b_restricted_noun.ct 1.951592 0.7373833 0.4448778 3.342794

mcmc_plot(prod_unattested_novel_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_unattested_novel_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1    0.0000
## 2    0.0065

Production data: Effect of statistical entrenchment

data_long_e <- gather(entrenchment_production.df, det_type, produced, det1:none, factor_key=TRUE)

data_long_e$noun_type_test <-recode(data_long_e$noun_type_test, "construction1" = "test: singular","construction2" = "test: plural")


p = ggplot(data_long_e, aes(x = noun_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme(panel.grid.major = element_blank()) +
  facet_grid("noun_type_test") +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("singular", "plural", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "singular", "construction2" = "plural", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("noun type at training")
p

#a. Are participants producing more attested than unattested dets?
# here, we want to see how often participants say the unattested e.g. transitive-only det1 for a det2 (intransitive-only) noun in the intransitive condition at test 
# and vice versa 

production_entrenchment_attested_unattested.df  <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, restricted_noun =="yes")

#We want to compare attested vs. unattested trials for transitive nouns in the intransitive inchoative construction at test
production_entrenchment_attested_unattested1.df  <- subset(production_entrenchment_attested_unattested.df, noun_type_training2 == "construction1" & noun_type_test == "construction2")

#And intransitive inchoative nouns in the transitive construction at test. Filter out irrelevant trials
production_entrenchment_attested_unattested2.df  <- subset(production_entrenchment_attested_unattested.df, noun_type_training2 == "construction2" & noun_type_test == "construction1")


production_entrenchment_attested_unattested.df <- rbind(production_entrenchment_attested_unattested1.df, production_entrenchment_attested_unattested2.df)

#How much of the time are participants producing attested items?
round(mean(production_entrenchment_attested_unattested.df$attested_unattested),3)

## [1] 0.123

# and separately for each noun type
round(tapply(production_entrenchment_attested_unattested.df$attested_unattested, production_entrenchment_attested_unattested.df$noun_type_training2, mean),3)

## construction1 construction2 
##         0.137         0.109

production_entrenchment_attested_unattested.df$noun_type_training2 <- factor(production_entrenchment_attested_unattested.df$noun_type_training2)
df_prod_ent = lizCenter((production_entrenchment_attested_unattested.df), list("noun_type_training2"))  


# maximally vague priors for the predictors and the intercept
prod_attested_unattested_ent = brm(formula = attested_unattested ~noun_type_training2.ct + (1 + noun_type_training2.ct|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent, variable = c("b_Intercept","b_noun_type_training2.ct"))

##                            Estimate Est.Error      Q2.5     Q97.5
## b_Intercept              -3.1011716 0.8132802 -4.548023 -1.287716
## b_noun_type_training2.ct -0.3592539 0.7048490 -1.723719  1.044658

mcmc_plot(prod_attested_unattested_ent, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_training2.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.9980833
## 2 0.3002500

#same analyses without noun_training_type


# maximally vague priors for the intercept
prod_attested_unattested_ent_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent_final, variable = c("b_Intercept"))

##              Estimate Est.Error      Q2.5     Q97.5
## b_Intercept -3.138102  0.659476 -4.374398 -1.757928

mcmc_plot(prod_attested_unattested_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0.9999167

# c. we will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the 2 non-alternating nouns than for the novel noun (presumably the “unwitnessed” form has to be set arbitrarily here)


production_entrenchment_restricted_novel.df <- subset(entrenchment_production.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_restricted_novel.df<- subset(production_entrenchment_restricted_novel.df, noun_type_training2 != "alternating")

# all forms are unwitnessed for the novel noun so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$noun_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_entrenchment_restricted_novel.df$attested_unattested)
production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$noun_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_entrenchment_restricted_novel.df$attested_unattested)

# select trials featuring the novel noun in plural scences
production_entrenchment_restricted_novel1.df <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "novel"  & noun_type_test == "construction2")


# Select trials featuring singular only nouns in plural scenes
production_entrenchment_restricted_novel2.df  <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "construction1" & noun_type_test == "construction2")

# Select trials featuring plural only nouns in singular scenes
production_entrenchment_restricted_novel3.df  <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "construction2" & noun_type_test == "construction1")


production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$noun_type_training2, mean),3)

## construction1 construction2         novel 
##         0.137         0.109         0.077

# reverse coding to focus on unattested rather than attested for novel vs. restricted
production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)
production_entrenchment_restricted_novel.df$attested_unattested<- recode(production_entrenchment_restricted_novel.df$attested_unattested, `1` = 0L, `0` = 1L)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$restricted_noun, mean),3)

##    no   yes 
## 0.923 0.877

#what this means is that participants produce *unattested forms* less for the restricted than they do for the novel

production_entrenchment_restricted_novel.df$restricted_noun <- factor(production_entrenchment_restricted_novel.df$restricted_noun)
production_entrenchment_restricted_novel1.df = lizCenter(production_entrenchment_restricted_novel.df, list("restricted_noun"))


# maximally vague priors for the predictors and the intercept
prod_unattested_novel_ent_final = brm(formula = attested_unattested ~restricted_noun.ct + (1 + restricted_noun.ct|participant_private_id), data=production_entrenchment_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(prod_unattested_novel_ent_final, variable = c("b_Intercept","b_restricted_noun.ct"))

##                        Estimate Est.Error      Q2.5     Q97.5
## b_Intercept           3.3952558 0.6658364  1.983958 4.6391261
## b_restricted_noun.ct -0.3182357 0.6604142 -1.644561 0.9749368

mcmc_plot(prod_unattested_novel_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_unattested_novel_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##      c(C1, C2)
## 1 8.333333e-05
## 2 3.108333e-01

Experiment 5

#Create the dataframes that we will be working on
combined_production_data.df <- read.csv("exp5_production_data.csv")

combined_judgment_data.df <- read.csv("exp5_judgment_data.csv")
combined_judgment_data.df$restricted_noun <- factor(combined_judgment_data.df$restricted_noun)
combined_judgment_data.df$condition <- factor(combined_judgment_data.df$condition)


#separately for entrenchment and preemption

#entrenchment
entrenchment_production.df <- subset(combined_production_data.df, condition == "entrenchment")
entrenchment_production.df$semantically_correct <- as.numeric(entrenchment_production.df$semantically_correct)
entrenchment_production.df$noun_type_test <- factor(entrenchment_production.df$noun_type_test)

# We want to first create some columns coding whether det1, det2, other and none was produced. we check proportions for each type of noun

entrenchment_production.df$det1 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction1", 1, 0)
entrenchment_production.df$det2 <- ifelse(entrenchment_production.df$det_lenient_adapted == "det_construction2", 1, 0)
entrenchment_production.df$other <- ifelse(entrenchment_production.df$det_lenient_adapted == "other", 1, 0)
entrenchment_production.df$none <- ifelse(entrenchment_production.df$det_lenient_adapted == "none", 1, 0)

entrenchment_judgment.df <- subset(combined_judgment_data.df, condition == "entrenchment")
entrenchment_judgment.df$semantically_correct <- factor(entrenchment_judgment.df$semantically_correct)
entrenchment_judgment.df$restricted_noun <- factor(entrenchment_judgment.df$restricted_noun)
entrenchment_judgment.df$noun_type_test <- factor(entrenchment_judgment.df$noun_type_test)

#preemption
preemption_production.df <- subset(combined_production_data.df, condition == "preemption")
preemption_production.df$semantically_correct <- as.numeric(preemption_production.df$semantically_correct)
preemption_production.df$noun_type_test <- factor(preemption_production.df$noun_type_test)

# We want to create some columns coding whether det1, det2, other and none was produced

preemption_production.df$det1 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction1", 1, 0)
preemption_production.df$det2 <- ifelse(preemption_production.df$det_lenient_adapted == "det_construction2", 1, 0)
preemption_production.df$other <- ifelse(preemption_production.df$det_lenient_adapted == "other", 1, 0)
preemption_production.df$none <- ifelse(preemption_production.df$det_lenient_adapted == "none", 1, 0)


preemption_judgment.df <- subset(combined_judgment_data.df, condition == "preemption")
preemption_judgment.df$semantically_correct <- factor(preemption_judgment.df$semantically_correct)
preemption_judgment.df$restricted_noun <- factor(preemption_judgment.df$restricted_noun)
preemption_judgment.df$noun_type_test <- factor(preemption_judgment.df$noun_type_test)

Preregistered data analyses

Question 1: Have participants picked up on the difference in meaning between the two argument-structure constructions?

Production data

#Figure X
RQ1_graph_productions.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "alternating" |noun_type_training2 == "novel")

# aggregated dataframe for means
aggregated.graph1 = aggregate(semantically_correct ~ noun_type_training2 + participant_private_id, RQ1_graph_productions.df, FUN=mean)

aggregated.graph1 <- rename(aggregated.graph1, noun = noun_type_training2,
                            correct = semantically_correct)

yarrr::pirateplot(formula = correct  ~ noun,
                  data = aggregated.graph1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "% semantically correct",
                  cex.lab = 1,
                  cex.axis = 1,
                  cex.names = 1,
                  yaxt = "n")

axis(2, at = seq(0, 1, by = 0.25), las=1)
abline(h = 0.50, lty = 2)

#1 alternating noun production

alternating_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "alternating")
#and filter out responses where participants said something other than det1 or det2
alternating_prod.df = subset(alternating_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_alternating_prod.df = aggregate(semantically_correct ~ noun_type_test + participant_private_id, alternating_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_alternating_prod.df$semantically_correct),3)

## [1] 0.958

# average accuracy separately for construction 1 (singular) and construction 2 (plural) scenes
round(tapply(aggregated.means_alternating_prod.df$semantically_correct, aggregated.means_alternating_prod.df$ noun_type_test, mean),3)

## construction1 construction2 
##         0.976         0.940

a = lizCenter(alternating_prod.df, list("noun_type_test"))  


# maximally vague priors for the intercept and the predictors
alternating_model <-brm(formula = semantically_correct~noun_type_test.ct + (1 + noun_type_test.ct|participant_private_id), data=a, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model, variable = c("b_Intercept", "b_noun_type_test.ct"))

##                       Estimate Est.Error      Q2.5     Q97.5
## b_Intercept          3.0892552 0.4308844  2.310566 3.9850731
## b_noun_type_test.ct -0.4485327 0.6173624 -1.658085 0.7572271

mcmc_plot(alternating_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(alternating_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_test.ct"] > 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1     0.000
## 2     0.235

# no difference between construction 1 and construction 2  
# Final model

# maximally vague priors for the intercept 
alternating_model_final = brm(formula = semantically_correct~1 + (1|participant_private_id), data=a, family = bernoulli(link = logit),set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(alternating_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5   Q97.5
## b_Intercept 2.971149 0.4027446 2.259124 3.84618

mcmc_plot(alternating_model_final, variable = "b_Intercept", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(alternating_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

#2 novel noun production

novel_prod.df = subset(entrenchment_production.df, condition == "entrenchment" & noun_type_training2 == "novel")
#and filter out responses where participants said something other than det1 or det2
novel_prod.df = subset(novel_prod.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")

# aggregated dataframe for means
aggregated.means_novel_prod.df = aggregate(semantically_correct ~ noun_type_test + participant_private_id, novel_prod.df, FUN=mean)

# average accuracy across trial types
round(mean(aggregated.means_novel_prod.df$semantically_correct),3)

## [1] 0.94

# average accuracy separately for construction 1 (singular) and construction 2 (plural) scenes
round(tapply(aggregated.means_novel_prod.df$semantically_correct, aggregated.means_novel_prod.df$noun_type_test, mean),3)

## construction1 construction2 
##          0.94          0.94

b = lizCenter(novel_prod.df, list("noun_type_test"))  

# maximally vague priors for the intercept and the predictors
novel_model <- brm(formula = semantically_correct~noun_type_test.ct + (1 + noun_type_test.ct|participant_private_id), data=b, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model, variable = c("b_Intercept", "b_noun_type_test.ct"))

##                        Estimate Est.Error      Q2.5    Q97.5
## b_Intercept          2.89078393 0.4470614  2.066223 3.824239
## b_noun_type_test.ct -0.04543544 0.5951371 -1.229992 1.116706

mcmc_plot(novel_model, variable = "^b_", regex = TRUE)


samps = as.matrix(as.mcmc(novel_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_test.ct"] < 0)
pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.5301667

# no difference between construction 1 and construction 2  
# Final model   

# maximally vague priors for the intercept 
novel_model_final <- brm(formula = semantically_correct~1+ (1|participant_private_id), data=b, family = bernoulli(link = logit), set_prior("normal(0,1)", class="Intercept"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(novel_model_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 2.792343 0.4229372 2.028788 3.700823

mcmc_plot(novel_model_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(novel_model_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0

Judgment data

There are no semantically incorrect trials in Experiment 5, thus, these analyses are not possible.

Question 2: Does statistical preemption constrain morphological generalizations in children (judgment data)?

#Figure X

#first, filter our semantically incorrect trials

judgments_unattested_novel.df <- subset(combined_judgment_data.df, semantically_correct == "yes")   

#we only want to keep novel

judgments_novel.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "novel")   

#and restricted items

judgment_unattested_constr1.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "construction1" & attested_unattested == "0")   
judgment_unattested_constr2.df <- subset(judgments_unattested_novel.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

judgment_unattested_novel.df <- rbind(judgments_novel.df, judgment_unattested_constr1.df, judgment_unattested_constr2.df)

aggregated.means = aggregate(response ~ condition + restricted_noun + participant_private_id, judgment_unattested_novel.df, FUN=mean)
aggregated.means<- rename(aggregated.means, restricted = restricted_noun)

yarrr::pirateplot(formula = response ~ restricted + condition,
                  data = aggregated.means,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

judgment_unattested_novel_preemption.df <- subset(judgment_unattested_novel.df, condition == "preemption")   
round(tapply(judgment_unattested_novel_preemption.df$response, judgment_unattested_novel_preemption.df$restricted_noun, mean),3)

##    no   yes 
## 3.141 2.363

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(judgment_unattested_novel_preemption.df, list("restricted_noun"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_preemption_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_preemption_model, variable = c("b_Intercept","b_restricted_noun.ct"))

##                        Estimate Est.Error      Q2.5      Q97.5
## b_Intercept           2.7567991 0.1400212  2.482302  3.0337493
## b_restricted_noun.ct -0.7476843 0.2055766 -1.155490 -0.3355538

mcmc_plot(judgments_preemption_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_preemption_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.9999167

# BF analyses: we use the difference between attested and novel in Experiment 1 (SD = 0.65) as an estimate of the difference we expect here

Bf(0.21, 0.75, uniform = 0, meanoftheory = 0, sdtheory = 0.65/2, tail = 1)

## $LikelihoodTheory
## [1] 0.3147016
## 
## $Likelihoodnull
## [1] 0.003228164
## 
## $BayesFactor
## [1] 97.48626

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.21, 0.75, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h1 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.07

print(high_threshold)

## [1] 4

Question 3: Does statistical entrenchment constrain morphological generalizations in children (judgment data)?

#no semantically incorrect trials here

#we only want to keep novel

entrenchment_judgment_novel.df <- subset(entrenchment_judgment.df, noun_type_training2 == "novel")   

#and restricted items

entrenchment_judgment_unattested_constr1.df <- subset(entrenchment_judgment.df, noun_type_training2 == "construction1" & attested_unattested == "0")

entrenchment_judgment_unattested_constr2.df <- subset(entrenchment_judgment.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

#bind new dataframe
entrenchment_judgment_unattested_novel.df <- rbind(entrenchment_judgment_novel.df, entrenchment_judgment_unattested_constr1.df, entrenchment_judgment_unattested_constr2.df)
entrenchment_judgment_unattested_novel.df$restricted_noun <- factor(entrenchment_judgment_unattested_novel.df$restricted_noun , levels = c("yes", "no"))

round(tapply(entrenchment_judgment_unattested_novel.df$response, entrenchment_judgment_unattested_novel.df$restricted_noun, mean),3)

##   yes    no 
## 4.267 4.242

#Center variables of interest using the lizCenter function:
d_unattested_novel = lizCenter(entrenchment_judgment_unattested_novel.df, list("restricted_noun"))

# maximally vague priors for the predictors (we don't interpret the intercept here) 
judgments_entrenchment_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct, data=d_unattested_novel, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))


posterior_summary(judgments_entrenchment_model, variable = c("b_Intercept","b_restricted_noun.ct"))

##                        Estimate Est.Error       Q2.5     Q97.5
## b_Intercept           4.2605402 0.1479714  3.9670195 4.5486003
## b_restricted_noun.ct -0.0204897 0.2136363 -0.4488669 0.3985278

mcmc_plot(judgments_entrenchment_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_entrenchment_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0000000
## 2 0.5370833

# this one is based on the final N in the adult study (**attested vs. unattested** used as a max)
Bf(0.21, -0.02, uniform = 0, meanoftheory = 0, sdtheory = 0.38/2, tail = 1)

## $LikelihoodTheory
## [1] 1.337387
## 
## $Likelihoodnull
## [1] 1.891129
## 
## $BayesFactor
## [1] 0.7071894

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.21, -0.02, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)


# find values for which BF is inconclusive 
ev_for_h1 <- subset(data.frame(range_test), BF < 3 & BF > 1/3)
low_threshold <- min(ev_for_h1$sdtheory)
high_threshold <- max(ev_for_h1$sdtheory)
print(low_threshold)

## [1] 0.01

print(high_threshold)

## [1] 0.54

# find out how many more participants we would need for conclusive evidence for entrenchment (BF > 3)
invisible(Bf_powercalc(0.21, -0.02, uniform=0, meanoftheory=0, sdtheory=0.38/2, tail=1, N=40, min=30, max=400))

#N = 327

Question 4: Is the effect of statistical pre-emption larger than entrenchment (judgment data)?

#all are semantically correct trials

#we only want to keep novel

all_judgment_novel.df <- subset(combined_judgment_data.df, noun_type_training2 == "novel")   

#and restricted items

all_judgment_unattested_constr1.df <- subset(combined_judgment_data.df, noun_type_training2 == "construction1" & attested_unattested == "0")   
all_judgment_unattested_constr2.df <- subset(combined_judgment_data.df, noun_type_training2 == "construction2" & attested_unattested == "0")   

all_judgment_unattested_novel.df <- rbind(all_judgment_novel.df, all_judgment_unattested_constr1.df, all_judgment_unattested_constr2.df)
all_judgment_unattested_novel.df$restricted_noun <- factor(all_judgment_unattested_novel.df$restricted_noun , levels = c("yes", "no"))

round(tapply(all_judgment_unattested_novel.df$response, list(all_judgment_unattested_novel.df$restricted_noun, all_judgment_unattested_novel.df$condition), mean),3)

##     entrenchment preemption
## yes        4.267      2.363
## no         4.242      3.141

round(tapply(all_judgment_unattested_novel.df$response, all_judgment_unattested_novel.df$condition, mean),3)

## entrenchment   preemption 
##        4.254        2.752

#Center variables of interest using the lizCenter function:
df = lizCenter(all_judgment_unattested_novel.df, list("restricted_noun", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
judgments_pre_vs_ent_model <- brm(formula = response~(1 +restricted_noun.ct|participant_private_id)+restricted_noun.ct * condition.ct, data=df, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(judgments_pre_vs_ent_model, variable = c("b_Intercept", "b_restricted_noun.ct","b_condition.ct", "b_restricted_noun.ct:condition.ct"))

##                                     Estimate  Est.Error        Q2.5      Q97.5
## b_Intercept                        3.4593249 0.09964572  3.26811895  3.6588504
## b_restricted_noun.ct               0.3893279 0.14860623  0.09554679  0.6764106
## b_condition.ct                    -1.4400544 0.19709344 -1.82936761 -1.0446270
## b_restricted_noun.ct:condition.ct  0.7370877 0.28891731  0.17763038  1.2945767

mcmc_plot(judgments_pre_vs_ent_model, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(judgments_pre_vs_ent_model))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)
C3=mean(samps[,"b_condition.ct"] > 0)
C4=mean(samps[,"b_restricted_noun.ct:condition.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1       0.000000000
## 2       0.005833333
## 3       0.000000000
## 4       0.005083333

#roughly predicted effect size from previous study was 1.0. Use it as an estimate of the max effect we expect here
Bf(0.29, 0.73, uniform = 0, meanoftheory = 0, sdtheory = 1.00/2, tail = 1)

## $LikelihoodTheory
## [1] 0.6125222
## 
## $Likelihoodnull
## [1] 0.05788388
## 
## $BayesFactor
## [1] 10.58191

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.29, 0.73, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.15

print(high_threshold)

## [1] 4

Exploratory data analyses

Exploratory data analyses

Effect of statistical pre-emption: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

# Figure x
judgments_unattested_attested.df <- subset(combined_judgment_data.df, semantically_correct == "yes")   
judgments_unattested_attested.df <- subset(judgments_unattested_attested.df, restricted_noun == "yes")   


aggregated.means1 = aggregate(response ~ condition + attested_unattested + participant_private_id, judgments_unattested_attested.df , FUN=mean)
aggregated.means1<- rename(aggregated.means1, attested = attested_unattested)

aggregated.means1$attested<- recode(aggregated.means1$attested, "1" = "yes","0" = "no")


yarrr::pirateplot(formula = response ~  attested  + condition,
                  data = aggregated.means1,
                  main = "",
                  theme=2,
                  point.o = .3,
                  gl.col = 'white',
                  ylab = "Rating",
                  cex.lab = 0.8,
                  cex.axis = 1,
                  cex.names = 0.8,
                  yaxt = "n")

axis(2, at = seq(1, 9, by = 1), las=1)

# analyses
attested_vs_unattested = subset(preemption_judgment.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested$response, attested_vs_unattested$attested_unattested, mean),3)

##     0     1 
## 2.363 4.526

#Center variables of interest using the lizCenter function:
d_attested_unattested = lizCenter(attested_vs_unattested , list("attested_unattested"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_preemption <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept              3.461733 0.08495619 3.295578 3.629453
## b_attested_unattested.ct 2.045099 0.23683486 1.569918 2.502679

mcmc_plot(attested_unattested_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1         0
## 2         0

# prior from previous study with adults: 2.55 as a max
Bf(0.24, 2.04, uniform = 0, meanoftheory = 0, sdtheory = 2.55/2 , tail = 1)

## $LikelihoodTheory
## [1] 0.1786466
## 
## $Likelihoodnull
## [1] 3.402598e-16
## 
## $BayesFactor
## [1] 5.250301e+14

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.24, 2.04, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.04

print(high_threshold)

## [1] 4

Effect of statistical entrenchment: Comparison of adults’ judgment ratings (acceptability) for witnessed versus unwitnessed forms

attested_vs_unattested_ent = subset(entrenchment_judgment.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested_ent$response, attested_vs_unattested_ent$attested_unattested, mean),3)

##     0     1 
## 4.267 4.717

#Center variables of interest using the lizCenter function:
d_attested_unattested_ent = lizCenter(attested_vs_unattested_ent, list("attested_unattested"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment <- brm(formula =response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct, data=d_attested_unattested_ent, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                           Estimate Est.Error      Q2.5     Q97.5
## b_Intercept              4.5008207  0.123330 4.2620827 4.7408144
## b_attested_unattested.ct 0.4352817  0.168091 0.1050292 0.7670532

mcmc_plot(attested_unattested_entrenchment, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 


pMCMC=as.data.frame(c(C1,C2))
pMCMC

##     c(C1, C2)
## 1 0.000000000
## 2 0.004916667

# expect a difference of 0.38 from previous work
Bf(0.17, 0.43, uniform = 0, meanoftheory = 0, sdtheory = 0.346/2, tail = 1)

## $LikelihoodTheory
## [1] 0.6592187
## 
## $Likelihoodnull
## [1] 0.0957568
## 
## $BayesFactor
## [1] 6.884301

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.17, 0.43, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.09

print(high_threshold)

## [1] 2.72

Entrenchment vs. preemption: ratings for witnessed vs. unwitnessed forms

attested_vs_unattested_across = subset(combined_judgment_data.df, restricted_noun == "yes")

round(tapply(attested_vs_unattested_across$response, list(attested_vs_unattested_across$condition, attested_vs_unattested_across$attested_unattested), mean),3)

##                  0     1
## entrenchment 4.267 4.717
## preemption   2.363 4.526

#Center variables of interest using the lizCenter function:
df_attested_unattested = lizCenter(attested_vs_unattested_across, list("attested_unattested", "condition"))

# maximally vague priors for the predictors (we don't interpret the intercept here)
attested_unattested_entrenchment_preemption <- brm(formula = response~(1 +attested_unattested.ct|participant_private_id)+attested_unattested.ct * condition.ct, data=df_attested_unattested, family=gaussian(),set_prior("normal(0,1)", class="b"),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(attested_unattested_entrenchment_preemption, variable = c("b_Intercept", "b_attested_unattested.ct"))

##                          Estimate  Est.Error     Q2.5    Q97.5
## b_Intercept              3.947383 0.07395113 3.802026 4.091341
## b_attested_unattested.ct 1.323288 0.14946990 1.030987 1.616206

mcmc_plot(attested_unattested_entrenchment_preemption, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(attested_unattested_entrenchment_preemption))

C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_attested_unattested.ct"] < 0) 
C3=mean(samps[,"b_condition.ct"] > 0) 
C4=mean(samps[,"b_attested_unattested.ct:condition.ct"] < 0) 

pMCMC=as.data.frame(c(C1,C2,C3,C4))
pMCMC

##   c(C1, C2, C3, C4)
## 1                 0
## 2                 0
## 3                 0
## 4                 0

#max predicted effect size from previous study 2.11
Bf(0.29, 1.56, uniform = 0, meanoftheory = 0, sdtheory = 2.12/2, tail = 1)

## $LikelihoodTheory
## [1] 0.2650902
## 
## $Likelihoodnull
## [1] 7.160226e-07
## 
## $BayesFactor
## [1] 370225.9

H1RANGE = seq(0,4,by=0.01)
range_test <- Bf_range(0.29, 1.56, meanoftheory=0, sdtheoryrange= H1RANGE, tail=1)

# find values for which BF > 3
ev_for_h0 <- subset(data.frame(range_test), BF > 3)
low_threshold <- min(ev_for_h0$sdtheory)
high_threshold <- max(ev_for_h0$sdtheory)
print(low_threshold)

## [1] 0.06

print(high_threshold)

## [1] 4

Production data: Effect of statistical pre-emption

# filter out missing data

data_long <- subset(preemption_production.df, experimenter == "GS")


data_long <- gather(data_long, det_type, produced, det1:none, factor_key=TRUE)

p = ggplot(data_long, aes(x = noun_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme_bw()+
  theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("construction1 ", "construction2", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "construction 1", "construction2" = "construction2", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("verb type at training")
p

#Are participants producing more attested than unattested dets? we will now compare proportion of attested dets (that's the intercept) for the restricted nouns against chance 

production_preemption_attested_unattested.df <- subset(preemption_production.df, experimenter == "GS")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_attested_unattested.df <- subset(production_preemption_attested_unattested.df, restricted_noun =="yes")

round(tapply(production_preemption_attested_unattested.df $attested_unattested, production_preemption_attested_unattested.df $noun_type_training2, mean),3)

## construction1 construction2 
##         0.940         0.944

#construction1 construction2 
#0.940         0.944


production_preemption_attested_unattested.df$noun_type_training2 <- factor(production_preemption_attested_unattested.df$noun_type_training2)

df_prod = lizCenter(production_preemption_attested_unattested.df , list("noun_type_training2"))  

# maximally vague priors for the predictors and the intercept
prod_attested_unattested = brm(formula = attested_unattested ~noun_type_training2.ct + (1 + noun_type_training2.ct|participant_private_id), data=df_prod, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested, variable = c("b_Intercept","b_noun_type_training2.ct"))

##                           Estimate Est.Error       Q2.5    Q97.5
## b_Intercept              3.1455491 0.9749168  0.8337138 4.706238
## b_noun_type_training2.ct 0.1470705 0.5889979 -1.0013175 1.311157

mcmc_plot(prod_attested_unattested, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_training2.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.0060000
## 2 0.4024167

#same analyses without noun_training_type

# maximally vague priors for the intercept
prod_attested_unattested_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_final, variable = c("b_Intercept"))

##             Estimate Est.Error     Q2.5    Q97.5
## b_Intercept 3.284324 0.8919063 1.097174 4.743904

mcmc_plot(prod_attested_unattested_final, variable = "^b_", regex = TRUE)

samps = as.matrix(as.mcmc(prod_attested_unattested_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0.001833333

# We will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the restricted verbs than for the novel verb


production_preemption_restricted_novel.df <- subset(preemption_production.df, experimenter == "GS")
production_preemption_restricted_novel.df <- subset(production_preemption_restricted_novel.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_preemption_restricted_novel.df<- subset(production_preemption_restricted_novel.df, noun_type_training2 != "alternating")

# all forms are unwitnessed for the novel verb so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$noun_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_preemption_restricted_novel.df$attested_unattested)

production_preemption_restricted_novel.df$attested_unattested <- ifelse(production_preemption_restricted_novel.df$noun_type_training2 == "novel" & production_preemption_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_preemption_restricted_novel.df$attested_unattested)

round(tapply(production_preemption_restricted_novel.df$attested_unattested, production_preemption_restricted_novel.df$noun_type_training2, mean),3)

## construction1 construction2         novel 
##         0.940         0.944         0.463

round(tapply(production_preemption_restricted_novel.df$attested_unattested , production_preemption_restricted_novel.df$restricted_noun, mean),3)

##    no   yes 
## 0.463 0.942

production_preemption_restricted_novel.df$restricted_noun <- factor(production_preemption_restricted_novel.df$restricted_noun)
production_preemption_restricted_novel1.df = lizCenter(production_preemption_restricted_novel.df, list("restricted_noun"))

# maximally vague priors for the predictors and the intercept
prod_unattested_novel_final = brm(formula = attested_unattested ~restricted_noun.ct + (1 + restricted_noun.ct|participant_private_id), data=production_preemption_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_final, variable = c("b_Intercept","b_restricted_noun.ct"))

##                      Estimate Est.Error        Q2.5    Q97.5
## b_Intercept          2.423371 0.5029382 1.346299953 3.349363
## b_restricted_noun.ct 1.882204 0.9328292 0.009959651 3.617911

mcmc_plot(prod_unattested_novel_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_unattested_novel_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##      c(C1, C2)
## 1 0.0003333333
## 2 0.0241666667

Production data: Effect of statistical entrenchment

# filter out missing data

data_long_e <- subset(entrenchment_production.df, experimenter == "GS")

data_long_e <- gather(data_long_e, det_type, produced, det1:none, factor_key=TRUE)

data_long_e$noun_type_test <-recode(data_long_e$noun_type_test, "construction1" = "test: singular","construction2" = "test: plural")


p = ggplot(data_long_e, aes(x = noun_type_training2, y = produced, fill = det_type)) +
  geom_bar(stat = "identity", position = "fill") +
  theme(panel.grid.major = element_blank()) +
  facet_grid("noun_type_test") +
  theme(panel.grid.minor = element_blank()) +
  scale_fill_manual(values=c("grey", "grey15", "azure3","azure4"), name="particle", labels=c("singular", "plural", "other", "none")) +
  scale_x_discrete(labels=c("alternating" = "alternating", "construction1" = "singular", "construction2" = "plural", "novel" = "novel")) +
  ylab("proportion produced") +
  xlab("noun type at training")
p

#a. Are participants producing more attested than unattested dets?
# here, we want to see how often participants say the unattested e.g. transitive-only det1 for a det2 (intransitive-only) noun in the intransitive condition at test 
# and vice versa 
production_entrenchment_attested_unattested.df <- subset(entrenchment_production.df, experimenter == "GS")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_attested_unattested.df  <- subset(production_entrenchment_attested_unattested.df, restricted_noun =="yes")

#We want to compare attested vs. unattested trials for transitive nouns in the intransitive inchoative construction at test
production_entrenchment_attested_unattested1.df  <- subset(production_entrenchment_attested_unattested.df, noun_type_training2 == "construction1" & noun_type_test == "construction2")

#And intransitive inchoative nouns in the transitive construction at test. Filter out irrelevant trials
production_entrenchment_attested_unattested2.df  <- subset(production_entrenchment_attested_unattested.df, noun_type_training2 == "construction2" & noun_type_test == "construction1")

production_entrenchment_attested_unattested.df <- rbind(production_entrenchment_attested_unattested1.df, production_entrenchment_attested_unattested2.df)

#How much of the time are participants producing attested items?
round(mean(production_entrenchment_attested_unattested.df$attested_unattested),3)

## [1] 0.167

# and separately for each noun type
round(tapply(production_entrenchment_attested_unattested.df$attested_unattested, production_entrenchment_attested_unattested.df$noun_type_training2, mean),3)

## construction1 construction2 
##         0.167         0.167

production_entrenchment_attested_unattested.df$noun_type_training2 <- factor(production_entrenchment_attested_unattested.df$noun_type_training2)
df_prod_ent = lizCenter((production_entrenchment_attested_unattested.df), list("noun_type_training2"))  


# maximally vague priors for the predictors and the intercept
prod_attested_unattested_ent = brm(formula = attested_unattested ~noun_type_training2.ct + (1 + noun_type_training2.ct|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)),cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent, variable = c("b_Intercept","b_noun_type_training2.ct"))

##                             Estimate Est.Error       Q2.5     Q97.5
## b_Intercept              -1.90249539 0.4239548 -2.7842160 -1.098850
## b_noun_type_training2.ct  0.06800616 0.4748525 -0.8465454  1.026435

mcmc_plot(prod_attested_unattested_ent, variable = "^b_", regex = TRUE)

dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_noun_type_training2.ct"] < 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##   c(C1, C2)
## 1 0.9999167
## 2 0.4418333

#same analyses without noun_training_type


# maximally vague priors for the intercept
prod_attested_unattested_ent_final = brm(formula = attested_unattested ~1 + (1|participant_private_id), data=df_prod_ent, family = bernoulli(link = logit), set_prior("normal(0, 1)", class = "Intercept"), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_attested_unattested_ent_final, variable = c("b_Intercept"))

##              Estimate Est.Error      Q2.5    Q97.5
## b_Intercept -1.868686 0.4132989 -2.726052 -1.09757

mcmc_plot(prod_attested_unattested_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_attested_unattested_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C1

## [1] 0.9999167

# c. we will now compare unattested for restricted vs. novel
# Do participants produce the unwitnessed form less for the 2 non-alternating nouns than for the novel noun (presumably the “unwitnessed” form has to be set arbitrarily here)

production_entrenchment_restricted_novel.df <- subset(entrenchment_production.df, experimenter == "GS")
production_entrenchment_restricted_novel.df <- subset(production_entrenchment_restricted_novel.df, det_lenient_adapted == "det_construction1" | det_lenient_adapted == "det_construction2")
production_entrenchment_restricted_novel.df<- subset(production_entrenchment_restricted_novel.df, noun_type_training2 != "alternating")

# all forms are unwitnessed for the novel noun so we are going to randomly set all det1s as attested and all dets2 as unattested 

production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$noun_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction1", 1, production_entrenchment_restricted_novel.df$attested_unattested)
production_entrenchment_restricted_novel.df$attested_unattested <- ifelse(production_entrenchment_restricted_novel.df$noun_type_training2 == "novel" & production_entrenchment_restricted_novel.df$det_lenient_adapted == "det_construction2", 0, production_entrenchment_restricted_novel.df$attested_unattested)

# select trials featuring the novel noun in the intransitive inchoative construction
production_entrenchment_restricted_novel1.df <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "novel"  & noun_type_test == "construction2")


# Select trials featuring transitive nouns in the intransitive inchoative construction at test
production_entrenchment_restricted_novel2.df  <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "construction1" & noun_type_test == "construction2")

# Select trials featuring intransitive nouns in the transitive construction at test
production_entrenchment_restricted_novel3.df  <- subset(production_entrenchment_restricted_novel.df, noun_type_training2 == "construction2" & noun_type_test == "construction1")


production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$noun_type_training2, mean),3)

## construction1 construction2         novel 
##         0.167         0.167         0.060

# reverse coding to focus on unattested rather than attested for novel vs. restricted
production_entrenchment_restricted_novel.df <- rbind(production_entrenchment_restricted_novel1.df, production_entrenchment_restricted_novel2.df, production_entrenchment_restricted_novel3.df)
production_entrenchment_restricted_novel.df$attested_unattested<- recode(production_entrenchment_restricted_novel.df$attested_unattested, `1` = 0L, `0` = 1L)


round(tapply(production_entrenchment_restricted_novel.df$attested_unattested , production_entrenchment_restricted_novel.df$restricted_noun, mean),3)

##    no   yes 
## 0.940 0.833

#what this means is that participants produce *unattested forms* less for the restricted than they do for the novel

production_entrenchment_restricted_novel.df$restricted_noun <- factor(production_entrenchment_restricted_novel.df$restricted_noun)
production_entrenchment_restricted_novel1.df = lizCenter(production_entrenchment_restricted_novel.df, list("restricted_noun"))


# maximally vague priors for the predictors and the intercept
prod_unattested_novel_ent_final = brm(formula = attested_unattested ~restricted_noun.ct + (1 + restricted_noun.ct|participant_private_id), data=production_entrenchment_restricted_novel1.df, family = bernoulli(link = logit), prior = c(prior(normal(0, 1), class = Intercept), prior(normal(0, 1), class = b)), cores=4, warmup = 2000, iter=5000, chains=4, control=list(adapt_delta = 0.99))

posterior_summary(prod_unattested_novel_ent_final, variable = c("b_Intercept","b_restricted_noun.ct"))

##                        Estimate Est.Error      Q2.5       Q97.5
## b_Intercept           2.2581055 0.4041438  1.482745  3.07371980
## b_restricted_noun.ct -0.9937755 0.5139680 -2.027903 -0.01354035

mcmc_plot(prod_unattested_novel_ent_final, variable = "^b_", regex = TRUE)
dev.off()

## null device 
##           1

samps = as.matrix(as.mcmc(prod_unattested_novel_ent_final))
C1=mean(samps[,"b_Intercept"] < 0)
C2=mean(samps[,"b_restricted_noun.ct"] > 0)

pMCMC=as.data.frame(c(C1,C2))
pMCMC

##    c(C1, C2)
## 1 0.00000000
## 2 0.02366667