Exp6Mandarin: Interim Analyses

see http://rpubs.com/ewonnacott/288885 for plan

Description of data collection up to this point

The data below were collected in June/July 2017. As discussed in our plan http://rpubs.com/ewonnacott/288885, in our initial sample, we anticipated being able to obtain around 20 participants in this first testing period, intending to inspect the data at that point. In fact, we obtained a final sample of N=26 in the Pictures only condition, and N=25 in the pictures + diacritics condition (further participants were tested but their data were excluded due to either (i) they were found to be speakers of a tonal language (ii) experimenter error during the testing process).

Notes on where/how analyses deviate from original plan

The data were analysed according to the plan in (http://rpubs.com/ewonnacott/288885) with the following exceptions:

1. There are places where we were planning to inform H1 using a value from pilot data with adults, and treating this as an estimate of the value of the theory (thus setting sd to of H1 to be equal to this value). However we noted in the plan that this value was actually higher than we would expect (and so would bias null and be very conservative). However we felt that this was better than setting them as the maximum of a uniform, since unlike the half normal this distribution does not favour smaller values. However in workshop at UCL (10th July 2017), Z Dienes recommended that an alternative when you have a maximum is to use a half normal but set the sd to be HALF the value of the maximum (so that 2*sd is the maximum but small values are still considered more likely). Cases where we now use a half value are:

1. Training prediction 3
   
2. AFC Pictures and AFC Diacritics, prediction 1 (though note here it is unclear whether the value should be regarded as an estimate or a maximum, so we do both)
3. Three interval oddity prediction 3
4. Three interval oddity additional questions: (though note here the maximum comes from the scale and actual scores at pre-test)

2. In the plan, there were places where we did not have a clear prediction and so considered two reversed predictions. In the plan, we proposed to use a two tailed test, since we didn’t have a clear direction. However following the workshop with Dienes (UCL, July 2010) we have become convinced that when there are two clear theories which could apply, it is more informative to test each separately with a different (one tailed) test, and we apply this approach below.

3. We had previously not planned to include random effects for item, due to concerns about power. However reading a recent paper (Brysbaert & Stevens: Power analysis and effect size in mixed effects models: A tutorial, pre-print at https://www.researchgate.net/project/When-do-we-have-enough-power-in-word-recognition-research/update/596c8612b53d2f270e5f71ff ) has made us revisit our understanding of power in mixed models, and we have decided to include “item” as a random effect. We include random intercepts throughout. For by item slopes, we test whether they contribute to the model following the approach recommended in this paper: (Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H., & Bates, D. (2017). Balancing Type I error and power in linear mixed models. Journal of Memory and Language, 94, 305-315). Given that we now include random effects for items, we do not include “version” however we do now include talker (coded speaker) for the training data (since item and talker were both previously subsumed in version) as a fixed effect and since it is known that individual talkers may be more/less difficult (NB: fixed rather than random factor since 4 levels is too few to code as a random effect). We use a coding that ensures that intercept and other fixed effects are reported such that they are averaged across all four talkers.

We continue to use the same by participants slope structure as in the plan (i.e. maximal slopes for factors of interest, not control factors)

4. There was an error in our coding of the control variable “tone contrast” in each of our models in the plan. Our aim was to use a coding which ensures that other fixed effects in the lme are reported so that they are averaged across all tone contrasts and we now use a coding that does this (in the models specified in the plan we erroneously applied a centring function which only works for categorical factors with two levels, or continuous factors). This has been corrected in the current models. Similarly, for the variable “trial type” in the 3 interval oddity task.

5. The error described in (4) above also applied to some of the analyses of pilot data used to get values for informing H1. We have therefore recalculated these values (changes were generally small).

6. We test some additional hypotheses which were not planned but which follow from the patterns we see in the data - these are clearly indicated.

Plan for ongoing data collection

In our original plan we stated that we would increase our sample size as necessary to obtain further evidence for/against our predictions, adding 10 participants at each step up to and including maximum sample of 50 per condition. We also said that we would determine whether or not to continue with the four session version of the study, or move to a two session version. This was dependent of whether there was evidence of a different pattern after four sessions, thus warranting the extra time and effort. In fact, in the data, for training only we do see evidence that a different pattern emerges when all four sessions are considered (if only sessions 1->2 are considered, it appears that there is greater learning in the P condition than the P+T condition; however there is emerging evidence that this benefit does not continue in the later sessions.

Given the above we have decided to continue with the four day procedure.

It will be seen below that not all tests of predictions have substantial evidence for either H0 or H1. Sample size analyses suggest different number of participants could be required to obtain such evidence for different tests and hypotheses (and in some cases, the numbers required are not feasible). In line with our plan, and since we have available resources and a school who are keen to take part, we are going to test another batch of 10 per condition.

We also plan to test a further 30 participants just on the introduction + 2AFC diactrics test (one session, no training).The motivation is to see whether participants are above chance on this test, even without training. This will help us to see the extent to which they actually pay attention to the diacrtics during training.Will test 30 participants (following analyses forAFC Diactrics Prediction 1 below, which suggested need about 27 to see that it is different from chance). We will use the estimate of the intercept for session2 from the model below (afcDIACRITICS.modB) to inform h1 - it will be treated as an estimate (so set sd to be this value) and the BF will compare evidence for the null against evidence for this value. We will also combining this data with the current data for that test to see if there is a main effect of condition. For this we will use the intercept to inform the H1 for the Bayes factor, seeting sd to be equal to the value of the intercept (see AFC picture Prediction see AFC pictures Prediction 2, Bayes Factor, for description of the logic). We acknowledge that it may be harder to get sufficient sample to have a BF <.33 or >3.

Load Packages

rm(list = ls(all = TRUE))
library(languageR)
library(lattice)
library(lme4)

## Loading required package: Matrix

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:stats':
## 
##     sigma

library(plotrix)
library(irr)

## Loading required package: lpSolve

library(plyr)
library(knitr)
library(stringdist)
library(ggplot2)
library(plyr)

Load Helper Functions

SummarySE

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

http://www.cookbook-r.com/Manipulating_data/Summarizing_data/

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

• data: a data frame.
• measurevar: the name of a column that contains the variable to be summarized.
• 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

From that website:

Summarizes data, handling within-subject variables by removing inter-subject variability. It will still work if there are no within-subject variables. Gives count, un-normed mean, normed mean (with same between-group mean), standard deviation, standard error of the mean, and confidence interval. If there are within-subject variables, it calculates adjusted values using the 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 function above. It can be found on the website “Cookbook for R”.

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

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 NAs

normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
                           na.rm=FALSE, .drop=TRUE) {
    #library(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 Jaeger’s 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 first level of the original factor, and the value above 0 will be the second 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 data frame.

• x: data frame
• listfname: a list of the variables to be centered (e.g. list(variable1,variable2))

The output is a copy of the data frame with a column (always a numeric variable) added for each of the centered variables. These columns are labelled with the each column’s previous name, but with “.ct” appended (e.g., “variable1” will become “variable1.ct”).

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 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,])}

Contrasts

lizContrasts

This function can be used used to create two centered dummy variables which stand in place of a three way factor (condition). This allows us to inspect each contrast separately, as well as their interactions with other factors. Other fixed effects in the model can be evaluated as the average effects across all levels of the factor.

The function takes a data frame (d), a factor from that database (condition), which must have three levels, and the name of the level of the factor which is to be used as the baseline for the contrasts (base level).

For example, if d is the data frame with a factor “condition” with three levels “lex_skew” lex_noskew" “mixed” then lizContrasts(d, d$condition, “lex_no_skew”) returns a data frame with two (numeric) columns added labelled “lex_noskew_VERSUS_lex_mixed” and “lex_noskew_VERSUS_lex_skew”. Wherever you would normally use “condition” in a formula in an lme, it can be replaced by (lex_noskew_VERSUS_lex_mixed + “lex_noskew_VERSUS_lex_skew) e.g. ~ (a * condition) becomes ~ (a * (lex_noskew_VERSUS_lex_mixed + lex_noskew_VERSUS_lex_skew)).

lizContrasts= function(d, condition, baselevel) 
{
  
    condition = factor(condition)
  condition = relevel(condition, baselevel)

    a= (contrasts(condition)-apply(contrasts(condition),2,mean))
    d$dummy1[condition== rownames(a)[1]] <- a[1] 
    d$dummy1[condition== rownames(a)[2]] <- a[2] 
    d$dummy1[condition== rownames(a)[3]] <- a[3] 
    
    d$dummy2[condition== rownames(a)[1]] <- a[4] 
    d$dummy2[condition== rownames(a)[2]] <- a[5] 
    d$dummy2[condition== rownames(a)[3]] <- a[6] 

    name1 = paste(baselevel, rownames(a)[2],sep="_VERSUS_")
    name2 = paste(baselevel, rownames(a)[3],sep="_VERSUS_")

    d[name1] = d$dummy1 
    d[name2] = d$dummy2 

    d$dummy1 <-NULL 
    d$dummy2 <-NULL 
    
    return(d)
}

lizContrasts4

This function is like lizContrasts but is used with a factor with four levels

lizContrasts4= function(d, condition, baselevel) 
{

  condition = factor(condition)
  condition = relevel(condition, baselevel)
  
    a= (contrasts(condition)- apply(contrasts(condition),2,mean))
    d$dummy1[condition== rownames(a)[1]] <- a[1] 
    d$dummy1[condition== rownames(a)[2]] <- a[2] 
    d$dummy1[condition== rownames(a)[3]] <- a[3] 
    d$dummy1[condition== rownames(a)[4]] <- a[4] 
    
    d$dummy2[condition== rownames(a)[1]] <- a[5] 
    d$dummy2[condition== rownames(a)[2]] <- a[6] 
    d$dummy2[condition== rownames(a)[3]] <- a[7] 
    d$dummy2[condition== rownames(a)[4]] <- a[8] 
    d$dummy3[condition== rownames(a)[1]] <- a[9] 
    d$dummy3[condition== rownames(a)[2]] <- a[10] 
    d$dummy3[condition== rownames(a)[3]] <- a[11] 
    d$dummy3[condition== rownames(a)[4]] <- a[12] 

    name1 = paste(baselevel, rownames(a)[2],sep="_VERSUS_")
    name2 = paste(baselevel, rownames(a)[3],sep="_VERSUS_")
    name3 = paste(baselevel, rownames(a)[4],sep="_VERSUS_")

    d[name1] = d$dummy1 
    d[name2] = d$dummy2 
    d[name3] = d$dummy3 

    d$dummy1 <-NULL 
    d$dummy2 <-NULL 
    d$dummy3 <-NULL 
    
    return(d)
}

lizContrasts6

This function is like lizContrasts but is used with a factor with six levels

lizContrasts6= function(d, condition, baselevel) 
{
  condition = factor(condition)
  condition = relevel(condition, baselevel)

    a= (contrasts(condition)- apply(contrasts(condition),2,mean))
    d$dummy1[condition== rownames(a)[1]] <- a[1] 
    d$dummy1[condition== rownames(a)[2]] <- a[2] 
    d$dummy1[condition== rownames(a)[3]] <- a[3] 
    d$dummy1[condition== rownames(a)[4]] <- a[4] 
    d$dummy1[condition== rownames(a)[5]] <- a[5] 
    d$dummy1[condition== rownames(a)[6]] <- a[6]    
    
    d$dummy2[condition== rownames(a)[1]] <- a[7] 
    d$dummy2[condition== rownames(a)[2]] <- a[8] 
    d$dummy2[condition== rownames(a)[3]] <- a[9] 
    d$dummy2[condition== rownames(a)[4]] <- a[10] 
    d$dummy2[condition== rownames(a)[5]] <- a[11] 
    d$dummy2[condition== rownames(a)[6]] <- a[12] 
    
    d$dummy3[condition== rownames(a)[1]] <- a[13] 
    d$dummy3[condition== rownames(a)[2]] <- a[14] 
    d$dummy3[condition== rownames(a)[3]] <- a[15] 
    d$dummy3[condition== rownames(a)[4]] <- a[16] 
    d$dummy3[condition== rownames(a)[5]] <- a[17] 
    d$dummy3[condition== rownames(a)[6]] <- a[18]
    
    d$dummy4[condition== rownames(a)[1]] <- a[19] 
    d$dummy4[condition== rownames(a)[2]] <- a[20] 
    d$dummy4[condition== rownames(a)[3]] <- a[21] 
    d$dummy4[condition== rownames(a)[4]] <- a[22] 
    d$dummy4[condition== rownames(a)[5]] <- a[23] 
    d$dummy4[condition== rownames(a)[6]] <- a[24] 
    
    d$dummy5[condition== rownames(a)[1]] <- a[25] 
    d$dummy5[condition== rownames(a)[2]] <- a[26] 
    d$dummy5[condition== rownames(a)[3]] <- a[27] 
    d$dummy5[condition== rownames(a)[4]] <- a[28] 
    d$dummy5[condition== rownames(a)[5]] <- a[29] 
    d$dummy5[condition== rownames(a)[6]] <- a[30] 

    name1 = paste(baselevel, rownames(a)[2],sep="_VERSUS_")
    name2 = paste(baselevel, rownames(a)[3],sep="_VERSUS_")
    name3 = paste(baselevel, rownames(a)[4],sep="_VERSUS_")
  name4 = paste(baselevel, rownames(a)[5],sep="_VERSUS_")
  name5 = paste(baselevel, rownames(a)[6],sep="_VERSUS_")

    d[name1] = d$dummy1 
    d[name2] = d$dummy2 
    d[name3] = d$dummy3 
  d[name4] = d$dummy4   
  d[name5] = d$dummy5   

    d$dummy1 <-NULL 
    d$dummy2 <-NULL 
    d$dummy3 <-NULL 
    d$dummy4 <-NULL 
    d$dummy5 <-NULL 

    return(d)
}

Bayesian Analyses

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=2){
 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_powercalc

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 it also requires the current number of participants (N)). However rather than return 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 some of 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)). However rather than setting sdtheory to a single value, it is set to a vector of values (sdrange) covering the desired range to test over.

Bf_range<-function(sd, obtained, meanoftheory=0, sdtheoryrange, tail=2)
{
 
 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) 
}

logodds

This function takes a percentage p and returns the logodds value

logodds <- function(p){log(p/(100-p))}  

Recalculation of values to be used to inform H1 drawn from pilot data

Load pilot data sets

setwd("pilotdata")
child.train.p = read.csv("kidstrain_clean.csv")
adult.train.p = read.csv ("adultstrain_clean.csv")
child.train.p$agegroup = relevel(child.train.p$agegroup, ref="7years")
adult.train.p$condition = relevel(adult.train.p$condition, ref="i")
adult.discrim.p = read.csv ("adultsdiscrim_clean.csv")

adult.discrim.p$pre_post = relevel(adult.discrim.p$pre_post, ref="pre")
adult.discrim.p$condition = relevel(adult.discrim.p$condition, ref="i")

Training: obtaining values to inform H1 for BF analyses of current training data

Obtaining values for Predictions 1 and 2:

child.train.7.p = subset(child.train.p, agegroup == "7years")
child.train.7.p <- lizCenter(child.train.7.p, list("session"))
child.train.7.p = lizContrasts6(child.train.7.p, child.train.7.p$tonecontrast, "T12")

child.train.7.p.mod = glmer (result ~
        + session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|subject),
        data = child.train.7.p, family = binomial, control=glmerControl(optimizer = "bobyqa"))

round(get_coeffs(child.train.7.p.mod, c("(Intercept)", "session.ct")),3)

##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.157      0.059   2.659    0.008
## session.ct     0.158      0.077   2.060    0.039

train.pred1.H1 = summary(child.train.7.p.mod )$coefficients["(Intercept)", "Estimate"]
train.pred2.H1 = summary(child.train.7.p.mod )$coefficients["session.ct", "Estimate"]

print(train.pred1.H1)

## [1] 0.1573752

print(train.pred2.H1)

## [1] 0.1578488

Note - previously had values pred1 = 0.15723615 pred2 = 0.1579623 (so near identical)

Values for Prediction 3.

adult.train.p <- lizCenter(adult.train.p, list("session", "condition"))
adult.train.p = lizContrasts6(adult.train.p, adult.train.p$tonecontrast, "T12")

adult.train.p.mod = glmer (result ~
        + session.ct * condition.ct
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|subject),
        data = adult.train.p, family = binomial, control=glmerControl(optimizer = "bobyqa"))

round(get_coeffs(adult.train.p.mod, c("condition.ct")),3)

##              Estimate Std. Error z value Pr(>|z|)
## condition.ct    0.431      0.185    2.33     0.02

Note- previously had value pred3= 0.4156183. Value from current model is very close: 0.430513 However as we note above, given that we assume this is a maxiumum values, we are going to be using half this value as sd for H1:

train.pred3.H1 = summary(adult.train.p.mod )$coefficients["condition.ct", "Estimate"]/2
print(train.pred3.H1)

## [1] 0.2152565

3 interval oddity: obtaining values to inform H1 for BF analyses of current training data

adult.discrim.p <- lizCenter(adult.discrim.p, list("pre_post", "condition"))
adult.discrim.p = lizContrasts6(adult.discrim.p, adult.discrim.p$tonecontrast, "T12")
adult.discrim.p = lizContrasts(adult.discrim.p, adult.discrim.p$trialtype, "fff")

adult.discrim.p.mod = glmer (correct ~
        + condition.ct * pre_post.ct
        
        + (fff_VERSUS_fmf + fff_VERSUS_ffm)
       + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (pre_post.ct|participantnumber),
        
         data = adult.discrim.p, family = binomial, control=glmerControl(optimizer = "bobyqa"))

discrim.pred1.H1 = summary(adult.discrim.p.mod )$coefficients["pre_post.ct", "Estimate"]/2

For predition 1. Note previously we had value 0.2868916. The equivalent vlaue from this model is 0.2296138.

We are going to use half this value as sd for H1, as it is clearly a maxium:

discrim.pred1.H1 = summary(adult.discrim.p.mod )$coefficients["pre_post.ct", "Estimate"]/2

Load Main Datasets

data for currrent study

training <- read.csv("training.csv", header=TRUE)
afc <- read.csv("afc.csv", header=TRUE)
discrim <- read.csv("discrim.csv", header=TRUE)
# wordrep <- read.csv("wordrep.csv", header=TRUE) - this data is currently being transcribed by native Mandarin speakers, and is not available for analysis at this point

Training

Inspect the Data

str(training)

## 'data.frame':    19584 obs. of  20 variables:
##  $ Order          : int  97 98 99 100 101 102 103 104 105 106 ...
##  $ pt_N           : Factor w/ 51 levels "Participant 10",..: 10 10 10 10 10 10 10 10 10 10 ...
##  $ program        : Factor w/ 13 levels "S1+2_ver1_P+T_trainingonly.cond\\",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ DEL            : Factor w/ 4 levels "[trainv1$]","[trainv2$]",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ soundfile      : Factor w/ 96 levels "chi_t2_f1a.wav",..: 33 33 93 65 29 81 89 25 25 37 ...
##  $ picture1       : Factor w/ 48 levels "chi_t2_spoon.jpg",..: 19 19 48 36 15 43 48 16 15 19 ...
##  $ picture2       : Factor w/ 48 levels "chi_t2_spoon.jpg",..: 20 20 47 35 16 44 47 15 16 20 ...
##  $ picture_correct: Factor w/ 2 levels "p1","p2": 1 1 1 2 2 1 2 2 1 2 ...
##  $ tone_contrast  : Factor w/ 6 levels "T12","T13","T14",..: 1 1 6 5 2 3 6 2 2 1 ...
##  $ speaker        : Factor w/ 4 levels "tv1","tv2","tv3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ session        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ DEL.1          : logi  NA NA NA NA NA NA ...
##  $ correct        : Factor w/ 2 levels "correct","incorrect": 1 1 2 2 1 2 2 2 1 2 ...
##  $ clicked        : Factor w/ 2 levels "clickedleft",..: 1 1 2 1 2 2 1 1 1 1 ...
##  $ RT             : num  1890 272 2938 3357 2877 ...
##  $ acc            : int  1 1 0 0 1 0 0 0 1 0 ...
##  $ version        : Factor w/ 2 levels "ver1","ver2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ cond           : Factor w/ 2 levels "P","P+T": 2 2 2 2 2 2 2 2 2 2 ...
##  $ item           : Factor w/ 24 levels "chit2","chit3",..: 9 9 24 17 8 21 23 7 7 10 ...
##  $ foil           : Factor w/ 48 levels "chi_t2_spoon.jpg",..: 20 20 47 36 15 44 48 16 16 19 ...

Get Means

t = summarySEwithin(data = aggregate(acc ~ pt_N + cond + session, FUN = mean, data = training), measurevar="acc", betweenvars = c("cond"), withinvars = ("session"))

## Automatically converting the following non-factors to factors: session

kable(t, digits = 4)

cond	session	N	acc	acc_norm	sd	se	ci
P	1	26	0.5337	0.5312	0.0980	0.0192	0.0396
P	2	26	0.6106	0.6081	0.1447	0.0284	0.0584
P	3	26	0.6242	0.6218	0.1488	0.0292	0.0601
P	4	26	0.6306	0.6282	0.1637	0.0321	0.0661
P+T	1	25	0.5596	0.5621	0.0953	0.0191	0.0393
P+T	2	25	0.5904	0.5929	0.1178	0.0236	0.0486
P+T	3	25	0.6121	0.6146	0.1436	0.0287	0.0593
P+T	4	25	0.6171	0.6196	0.1449	0.0290	0.0598

Plot Means

p1 = ggplot(t, aes(x = session, y = acc, color = cond)) + geom_point()
p1 = p1 + geom_errorbar(aes(ymin = acc-ci, ymax = acc+ci), width = .2)
p1 = p1 + xlab("Session") + ylab("Mean accuracy")
p1 = p1 + theme_bw()
p1 = p1 + coord_cartesian(ylim = c(0.5, 0.7))
p1 = p1 + ggtitle('Training')
p1 = p1+ theme(panel.background = element_rect(fill = 'white'),
               panel.grid.major = element_line(colour = "white"))
print(p1)

Main LME

Note on control variables: we include tone contrast (1-2, 1-3, 1-4, 2-3, 2-4, 3-4) and speaker (tv1, tv2, tv4) as fixed effects. Both could contribute to the model but neither is of specific interest. ### All FOUR sessions

training4 <- lizCenter(training, list("cond", "session", "version"))
training4 = lizContrasts6(training4, training4$tone_contrast, "T12")
training4 = lizContrasts4(training4, training4$speaker, "tv1")

training4.mod1 = glmer (acc ~
        + cond.ct * session.ct 
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N),
        
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

# converges - now find slope structure for items

training4.mod2 = glmer (acc ~
        + cond.ct * session.ct 
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct* cond.ct|item),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

training4.mod3 = glmer (acc ~
        + cond.ct * session.ct 
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct* cond.ct||item),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(training4.mod2, training4.mod3)

## Data: training4
## Models:
## training4.mod3: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod3:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod3:     (session.ct * cond.ct || item)
## training4.mod2: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod2:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod2:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod2:     (session.ct * cond.ct | item)
##                Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## training4.mod3 19 25506 25656 -12734    25468                         
## training4.mod2 25 25516 25713 -12733    25466 2.3674      6      0.883

## p>.2 so move to simpler slope structure

training4.mod4 = glmer (acc ~
        + cond.ct * session.ct 
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct+ cond.ct||item),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(training4.mod3, training4.mod4)

## Data: training4
## Models:
## training4.mod4: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod4:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod4:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod4:     (session.ct + cond.ct || item)
## training4.mod3: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod3:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod3:     (session.ct * cond.ct || item)
##                Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## training4.mod4 18 25507 25649 -12735    25471                         
## training4.mod3 19 25506 25656 -12734    25468 2.6461      1     0.1038

## p<.2 so move to stick with more complex slope structure

training4.mod = training4.mod3

round(summary(training4.mod)$coefficients,3)

##                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)           0.434      0.072   6.030    0.000
## cond.ct              -0.035      0.144  -0.247    0.805
## session.ct            0.126      0.026   4.849    0.000
## tv1_VERSUS_tv2       -0.020      0.042  -0.466    0.641
## tv1_VERSUS_tv3        0.013      0.042   0.297    0.766
## tv1_VERSUS_tv4        0.060      0.042   1.403    0.161
## T12_VERSUS_T13        0.230      0.095   2.412    0.016
## T12_VERSUS_T14       -0.037      0.095  -0.390    0.697
## T12_VERSUS_T23        0.198      0.095   2.076    0.038
## T12_VERSUS_T24       -0.010      0.095  -0.106    0.915
## T12_VERSUS_T34        0.162      0.095   1.699    0.089
## cond.ct:session.ct   -0.057      0.053  -1.075    0.282

First TWO sessions only

training2 <- droplevels(subset(training, session != "3"))
training2 <- droplevels(subset(training2, session != "4"))

training2 <- lizCenter(training2, list("cond", "session"))
training2 = lizContrasts6(training2, training2$tone_contrast, "T12")
training2 = lizContrasts4(training2, training2$speaker, "tv1")


training2.mod = glmer (acc ~
         + cond.ct * session.ct 
       + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct+ cond.ct||item),
        data = training2, family = binomial, control = glmerControl(optimizer = "bobyqa"))

round(summary(training2.mod)$coefficients,3)

##                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)           0.317      0.062   5.139    0.000
## cond.ct               0.001      0.124   0.009    0.993
## session.ct            0.249      0.067   3.703    0.000
## tv1_VERSUS_tv2        0.007      0.059   0.118    0.906
## tv1_VERSUS_tv3       -0.021      0.059  -0.354    0.723
## tv1_VERSUS_tv4        0.077      0.059   1.303    0.193
## T12_VERSUS_T13        0.127      0.110   1.148    0.251
## T12_VERSUS_T14       -0.040      0.110  -0.366    0.714
## T12_VERSUS_T23        0.132      0.110   1.199    0.231
## T12_VERSUS_T24       -0.116      0.110  -1.055    0.291
## T12_VERSUS_T34        0.054      0.110   0.490    0.624
## cond.ct:session.ct   -0.206      0.105  -1.970    0.049

Prediction 1: Children will show above chance performance in both condtions

We look across all four sessions, and just for first two sessions.

All FOUR sessions

Frequentist

We refit the models above with separate intercepts for each level of condition.

training4.modA = glmer (acc ~
        + cond
        + cond.ct * session.ct
        - 1
        - cond.ct
        
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct* cond.ct||item),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(training4.mod, training4.modA)

## Data: training4
## Models:
## training4.mod: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod:     (session.ct * cond.ct || item)
## training4.modA: acc ~ +cond + cond.ct * session.ct - 1 - cond.ct + (tv1_VERSUS_tv2 + 
## training4.modA:     tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## training4.modA:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## training4.modA:     pt_N) + (session.ct * cond.ct || item)
##                Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## training4.mod  19 25506 25656 -12734    25468                        
## training4.modA 19 25506 25656 -12734    25468     0      0          1

round(get_coeffs(training4.modA, c("condP", "condP+T")),3)

##         Estimate Std. Error z value Pr(>|z|)
## condP      0.451      0.101   4.479        0
## condP+T    0.416      0.103   4.050        0

Model comparison establishes that we are looking at the same model.

Bayes Factor

In the original plan, we said that we would use estimate from pilot data to inform H1, but that we could also use conditions to inform H1 for each other (and this could be a more accurate estimate). Here we do both for each condition.

1. Pictures-Only Condition

1. Using value from pilot data

meanBF = summary(training4.modA)$coefficients["condP", "Estimate"]
seBF = summary(training4.modA)$coefficients["condP", "Std. Error"]
h1mean = train.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.2316063
## 
## $Likelihoodnull
## [1] 0.0001743901
## 
## $BayesFactor
## [1] 1328.093

2. Using value from the other condition (i.e. Pictures + Diacritics Condition)

meanBF = summary(training4.modA)$coefficients["condP", "Estimate"]
seBF = summary(training4.modA)$coefficients["condP", "Std. Error"]
h1mean = summary(training4.modA)$coefficients["condP+T", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.069674
## 
## $Likelihoodnull
## [1] 0.0001743901
## 
## $BayesFactor
## [1] 6133.797

2. Pictures + Diacritics Condition

1. Using value from pilot data

meanBF = summary(training4.modA)$coefficients["condP+T", "Estimate"]
seBF = summary(training4.modA)$coefficients["condP+T", "Std. Error"]
h1mean = train.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.3674586
## 
## $Likelihoodnull
## [1] 0.001065388
## 
## $BayesFactor
## [1] 344.906

2. Using value from the other condition (i.e. Pictures-Only Condition)

meanBF = summary(training4.modA)$coefficients["condP+T", "Estimate"]
seBF = summary(training4.modA)$coefficients["condP+T", "Std. Error"]
h1mean = summary(training4.modA)$coefficients["condP", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.151943
## 
## $Likelihoodnull
## [1] 0.001065388
## 
## $BayesFactor
## [1] 1081.243

First TWO sessions only

Frequentist

Re fit model with separate interepts for each condition.

training2.modA = glmer (acc ~
        + cond
        + cond.ct * session.ct
        - 1
        - cond.ct
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct*cond.ct||item),
        data = training2, family = binomial, control = glmerControl(optimizer = "bobyqa"))

summary(training2.mod)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 +  
##     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 +  
##     T12_VERSUS_T24 + T12_VERSUS_T34) + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 +  
##     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 +  
##     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) +  
##     (session.ct + cond.ct || item)
##    Data: training2
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##  13112.3  13241.7  -6538.1  13076.3     9774 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3971 -1.0513  0.6245  0.8835  1.3092 
## 
## Random effects:
##  Groups Name        Variance Std.Dev. Corr
##  pt_N   (Intercept) 0.14247  0.3775       
##         session.ct  0.04968  0.2229   0.91
##  item   (Intercept) 0.01378  0.1174       
##  item.1 session.ct  0.04272  0.2067       
##  item.2 cond.ct     0.05583  0.2363       
## Number of obs: 9792, groups:  pt_N, 51; item, 24
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.317305   0.061743   5.139 2.76e-07 ***
## cond.ct             0.001119   0.123588   0.009 0.992777    
## session.ct          0.249179   0.067287   3.703 0.000213 ***
## tv1_VERSUS_tv2      0.006976   0.059029   0.118 0.905933    
## tv1_VERSUS_tv3     -0.020890   0.058981  -0.354 0.723204    
## tv1_VERSUS_tv4      0.077087   0.059175   1.303 0.192679    
## T12_VERSUS_T13      0.126599   0.110243   1.148 0.250819    
## T12_VERSUS_T14     -0.040256   0.110007  -0.366 0.714411    
## T12_VERSUS_T23      0.132188   0.110281   1.199 0.230667    
## T12_VERSUS_T24     -0.115964   0.109876  -1.055 0.291238    
## T12_VERSUS_T34      0.054036   0.110221   0.490 0.623955    
## cond.ct:session.ct -0.205992   0.104572  -1.970 0.048854 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) cnd.ct sssn.c t1_VERSUS_2 t1_VERSUS_3 t1_VERSUS_4
## cond.ct        -0.001                                                  
## session.ct      0.371 -0.003                                           
## t1_VERSUS_2     0.000  0.000  0.000                                    
## t1_VERSUS_3     0.000  0.000  0.000  0.500                             
## t1_VERSUS_4     0.001  0.000  0.001  0.499       0.499                 
## T12_VERSUS_T13  0.001  0.001  0.000  0.000       0.000       0.000     
## T12_VERSUS_T14 -0.001  0.001 -0.001  0.000       0.000       0.000     
## T12_VERSUS_T23  0.001  0.002  0.000  0.000       0.000       0.000     
## T12_VERSUS_T24 -0.002  0.001 -0.002  0.000       0.000       0.000     
## T12_VERSUS_T3   0.001  0.001  0.002  0.000       0.000       0.000     
## cnd.ct:sss.    -0.004  0.476 -0.005  0.000       0.000      -0.001     
##                T12_VERSUS_T13 T12_VERSUS_T14 T12_VERSUS_T23 T12_VERSUS_T24
## cond.ct                                                                   
## session.ct                                                                
## t1_VERSUS_2                                                               
## t1_VERSUS_3                                                               
## t1_VERSUS_4                                                               
## T12_VERSUS_T13                                                            
## T12_VERSUS_T14  0.500                                                     
## T12_VERSUS_T23  0.499          0.500                                      
## T12_VERSUS_T24  0.500          0.501          0.500                       
## T12_VERSUS_T3   0.499          0.500          0.499          0.500        
## cnd.ct:sss.     0.000          0.001          0.001          0.001        
##                T12_VERSUS_T3
## cond.ct                     
## session.ct                  
## t1_VERSUS_2                 
## t1_VERSUS_3                 
## t1_VERSUS_4                 
## T12_VERSUS_T13              
## T12_VERSUS_T14              
## T12_VERSUS_T23              
## T12_VERSUS_T24              
## T12_VERSUS_T3               
## cnd.ct:sss.     0.001

anova(training2.mod, training2.modA)

## Data: training2
## Models:
## training2.mod: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training2.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training2.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training2.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training2.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training2.mod:     (session.ct + cond.ct || item)
## training2.modA: acc ~ +cond + cond.ct * session.ct - 1 - cond.ct + (tv1_VERSUS_tv2 + 
## training2.modA:     tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## training2.modA:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (tv1_VERSUS_tv2 + 
## training2.modA:     tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## training2.modA:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## training2.modA:     pt_N) + (session.ct * cond.ct || item)
##                Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## training2.mod  18 13112 13242 -6538.1    13076                        
## training2.modA 19 13114 13251 -6538.1    13076     0      1          1

round(get_coeffs(training2.modA, c("condP", "condP+T")),3)

##         Estimate Std. Error z value Pr(>|z|)
## condP      0.317      0.087   3.660        0
## condP+T    0.318      0.088   3.606        0

Bayes Factor

In the original plan, we said that we would use estimate from pilot data to inform H1, but we could also use conditions to inform H1 for each other (and this might be a more accurate estimate). Here we do both for each condition.

1. Pictures-Only Condition

1. Using value from pilot data

meanBF = summary(training2.modA)$coefficients["condP", "Estimate"]
seBF = summary(training2.modA)$coefficients["condP", "Std. Error"]
h1mean = train.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.9374232
## 
## $Likelihoodnull
## [1] 0.005693753
## 
## $BayesFactor
## [1] 164.6407

2. Using value from the other condition (i.e. Pictures + Diacritics Condition)

meanBF = summary(training2.modA)$coefficients["condP", "Estimate"]
seBF = summary(training2.modA)$coefficients["condP", "Std. Error"]
h1mean = summary(training2.modA)$coefficients["condP+T", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.525183
## 
## $Likelihoodnull
## [1] 0.005693753
## 
## $BayesFactor
## [1] 267.8696

2. Pictures + Diacritics Condition

1. Using value from pilot data

meanBF = summary(training2.modA)$coefficients["condP+T", "Estimate"]
seBF = summary(training2.modA)$coefficients["condP+T", "Std. Error"]
h1mean =train.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.9355348
## 
## $Likelihoodnull
## [1] 0.006798014
## 
## $BayesFactor
## [1] 137.6188

2. Using value from the other condition (i.e. Pictures-Only Condition)

meanBF = summary(training2.modA)$coefficients["condP+T", "Estimate"]
seBF = summary(training2.modA)$coefficients["condP+T", "Std. Error"]
h1mean = summary(training2.modA)$coefficients["condP", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.520348
## 
## $Likelihoodnull
## [1] 0.006798014
## 
## $BayesFactor
## [1] 223.6458

Prediction 1 Summary

YES. Children show above chance performance in both conditions, both across all four sessions and across just two sessions;
BF > 10 in all cases (strong evidence). ## Prediction 2: Children will show improvement across sessions in both conditions

All FOUR sessions

Frequentist

Refit model with separate slopes for session for each level of condition.

training4.modB = glmer (acc ~
        + session.ct: cond
        + cond.ct * session.ct
        - session.ct
        - session.ct: cond.ct
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct* cond.ct||item),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(training4.mod, training4.modB)

## Data: training4
## Models:
## training4.mod: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training4.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training4.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training4.mod:     (session.ct * cond.ct || item)
## training4.modB: acc ~ +session.ct:cond + cond.ct * session.ct - session.ct - 
## training4.modB:     session.ct:cond.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + 
## training4.modB:     (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + 
## training4.modB:         T12_VERSUS_T34) + (session.ct | pt_N) + (session.ct * 
## training4.modB:     cond.ct || item)
##                Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## training4.mod  19 25506 25656 -12734    25468                        
## training4.modB 19 25506 25656 -12734    25468     0      0          1

round(get_coeffs(training4.modB, c("session.ct:condP+T", "session.ct:condP")),3)

##                    Estimate Std. Error z value Pr(>|z|)
## session.ct:condP+T    0.098      0.037   2.611    0.009
## session.ct:condP      0.154      0.037   4.192    0.000

Bayes Factor

In our original plan, we said that we would use estimate from pilot data to inform H1, but that we could also use conditions to inform H1 for each other (and these could provide a more accurate estimate). Here we do both for each condition.

1. Pictures-Only Condition

1. Using value from pilot data

meanBF = summary(training4.modB)$coefficients["session.ct:condP", "Estimate"]
seBF = summary(training4.modB)$coefficients["session.ct:condP", "Std. Error"]
h1mean = train.pred2.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 3.130413
## 
## $Likelihoodnull
## [1] 0.001654055
## 
## $BayesFactor
## [1] 1892.569

2. Using value from the other condition (i.e. Pictures + Diacritics)

meanBF = summary(training4.modB)$coefficients["session.ct:condP", "Estimate"]
seBF = summary(training4.modB)$coefficients["session.ct:condP", "Std. Error"]
h1mean = summary(training4.modB)$coefficients["session.ct:condP+T", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.562369
## 
## $Likelihoodnull
## [1] 0.001654055
## 
## $BayesFactor
## [1] 1549.144

2. Pictures + Diacritics Condition

1. Using value from pilot data

meanBF = summary(training4.modB)$coefficients["session.ct:condP+T", "Estimate"]
seBF = summary(training4.modB)$coefficients["session.ct:condP+T", "Std. Error"]
h1mean = train.pred2.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 4.081803
## 
## $Likelihoodnull
## [1] 0.3532688
## 
## $BayesFactor
## [1] 11.55438

2. Using value from the other condition (i.e. Pictures-Only Condition)

meanBF = summary(training4.modB)$coefficients["session.ct:condP+T", "Estimate"]
seBF = summary(training4.modB)$coefficients["session.ct:condP+T", "Std. Error"]
h1mean = summary(training4.modB)$coefficients["session.ct:condP", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 4.139983
## 
## $Likelihoodnull
## [1] 0.3532688
## 
## $BayesFactor
## [1] 11.71907

First TWO sessions only

Frequentist

Re fit model with separate slopes for session for each level of condition.

training2.modB = glmer (acc ~
        + session.ct: cond
        + cond.ct * session.ct 
        - session.ct
        - session.ct: cond.ct
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (session.ct* cond.ct||item),
        
        data = training2, family = binomial, control = glmerControl(optimizer = "bobyqa"))


anova(training2.mod, training2.modB)

## Data: training2
## Models:
## training2.mod: acc ~ +cond.ct * session.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training2.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training2.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + 
## training2.mod:     tv1_VERSUS_tv4) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## training2.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## training2.mod:     (session.ct + cond.ct || item)
## training2.modB: acc ~ +session.ct:cond + cond.ct * session.ct - session.ct - 
## training2.modB:     session.ct:cond.ct + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + 
## training2.modB:     (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + 
## training2.modB:         T12_VERSUS_T34) + (session.ct | pt_N) + (session.ct * 
## training2.modB:     cond.ct || item)
##                Df   AIC   BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## training2.mod  18 13112 13242 -6538.1    13076                        
## training2.modB 19 13114 13251 -6538.1    13076     0      1          1

round(get_coeffs(training2.modB, c("session.ct:condP+T", "session.ct:condP")),3)

##                    Estimate Std. Error z value Pr(>|z|)
## session.ct:condP+T    0.144      0.086   1.683    0.092
## session.ct:condP      0.350      0.085   4.130    0.000

Bayes Factor

In the original plan, we said that we would use estimate from pilot data to inform H1, but also that we could use conditions to inform H1 for each other (and that this could be a more accurate estimate). Here we do both for each condition

1. Picture-Only Condition

1. Using value from pilot data

meanBF= summary(training2.modB)$coefficients["session.ct:condP", "Estimate"]
seBF = summary(training2.modB)$coefficients["session.ct:condP", "Std. Error"]
h1mean = train.pred2.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.6596279
## 
## $Likelihoodnull
## [1] 0.000930837
## 
## $BayesFactor
## [1] 708.6395

2. Using value from the other condition (i.e. Pictures + Diacritics Condition)

meanBF = summary(training2.modB)$coefficients["session.ct:condP", "Estimate"]
seBF = summary(training2.modB)$coefficients["session.ct:condP", "Std. Error"]
h1mean = summary(training2.modB)$coefficients["session.ct:condP+T", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.5329234
## 
## $Likelihoodnull
## [1] 0.000930837
## 
## $BayesFactor
## [1] 572.5206

2. Pictures + Diacritics Condition

1. Using value from pilot data

meanBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Estimate"]
seBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Std. Error"]
h1mean = train.pred2.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.99311
## 
## $Likelihoodnull
## [1] 1.129654
## 
## $BayesFactor
## [1] 2.649581

2. Using value from the other condition (i.e. Pictures-Only Condition)

meanBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Estimate"]
seBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Std. Error"]
h1mean = summary(training2.modB)$coefficients["session.ct:condP", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.94129
## 
## $Likelihoodnull
## [1] 1.129654
## 
## $BayesFactor
## [1] 1.718481

For the Pictures + Diacritics condition we don’t yet have substantial evidence evidence for either H1 or H0 when looking at only Sessions 1 and 2 when basing it on value from within experiment (which is a better etimate than the pilot)

How big a sample might we need for this condition (current sample size = 25)?

meanBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Estimate"]
seBF = summary(training2.modB)$coefficients["session.ct:condP+T", "Std. Error"]
h1mean = summary(training2.modB)$coefficients["session.ct:condP", "Estimate"]

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 26, max = 100))

##      x         y
## 1   26  1.791457
## 2   27  1.868266
## 3   28  1.949092
## 4   29  2.034129
## 5   30  2.123585
## 6   31  2.217678
## 7   32  2.316639
## 8   33  2.420713
## 9   34  2.530158
## 10  35  2.645249
## 11  36  2.766272
## 12  37  2.893533
## 13  38  3.027354
## 14  39  3.168072
## 15  40  3.316045
## 16  41  3.471650
## 17  42  3.635284
## 18  43  3.807368
## 19  44  3.988342
## 20  45  4.178672
## 21  46  4.378850
## 22  47  4.589393
## 23  48  4.810847
## 24  49  5.043788
## 25  50  5.288821
## 26  51  5.546587
## 27  52  5.817758
## 28  53  6.103046
## 29  54  6.403200
## 30  55  6.719008
## 31  56  7.051303
## 32  57  7.400962
## 33  58  7.768911
## 34  59  8.156123
## 35  60  8.563628
## 36  61  8.992508
## 37  62  9.443906
## 38  63  9.919026
## 39  64 10.419138
## 40  65 10.945582
## 41  66 11.499769
## 42  67 12.083187
## 43  68 12.697406
## 44  69 13.344080
## 45  70 14.024955
## 46  71 14.741870
## 47  72 15.496766
## 48  73 16.291688
## 49  74 17.128796
## 50  75 18.010363
## 51  76 18.938789
## 52  77 19.916605
## 53  78 20.946480
## 54  79 22.031229
## 55  80 23.173820
## 56  81 24.377383
## 57  82 25.645222
## 58  83 26.980820
## 59  84 28.387849
## 60  85 29.870185
## 61  86 31.431916
## 62  87 33.077353
## 63  88 34.811043
## 64  89 36.637784
## 65  90 38.562637
## 66  91 40.590940
## 67  92 42.728325
## 68  93 44.980733
## 69  94 47.354431
## 70  95 49.856033
## 71  96 52.492512
## 72  97 55.271230
## 73  98 58.199951
## 74  99 61.286868
## 75 100 64.540622

Estimates could have BF>3 when N= 38 (i.e. another 13 in this condition_

Prediction 2 Summary

Pictures-Only Condition: Yes, there is substantial evidence that performance increases over session regardless of whether we consider all four sessions or only data from sessions 1 and 2 (BFs > 10).

Pictures + Diacritics Condition: If we consider data from all four sessions, there is strong evidence for improvement over sessions (BF>10). However, if we consider only sessions 1 and 2 we do not yet have substantial evidence. Sample estimation suggests could get this with another 13 participants in this condition.

Prediction 3: Children will show greater accuracy in the Pictures + Diacritics condition

All FOUR sessions

Frequentist

round(get_coeffs(training4.mod, c("cond.ct")),3)

##         Estimate Std. Error z value Pr(>|z|)
## cond.ct   -0.035      0.144  -0.247    0.805

Bayes Factor

The value we use to inform H1 here comes from prior adult data (since this is a maximum, it has been set to half the value; see above). Note that our mean here is in the opposite value to our prediction.

meanBF = summary(training4.mod)$coefficients["cond.ct", "Estimate"]
seBF = summary(training4.mod)$coefficients["cond.ct", "Std. Error"]
h1mean = train.pred3.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.284046
## 
## $Likelihoodnull
## [1] 2.691868
## 
## $BayesFactor
## [1] 0.4770091

We don’t yet have substantial evidence for either H1 (greater accuracy in the Pictures + Diacritics condition) or H0 (that there is no difference between conditions)

How much power might we need?

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 26, max = 100))

##      x         y
## 1   26 0.4687659
## 2   27 0.4608705
## 3   28 0.4532997
## 4   29 0.4460323
## 5   30 0.4390489
## 6   31 0.4323319
## 7   32 0.4258651
## 8   33 0.4196337
## 9   34 0.4136241
## 10  35 0.4078236
## 11  36 0.4022207
## 12  37 0.3968047
## 13  38 0.3915656
## 14  39 0.3864941
## 15  40 0.3815819
## 16  41 0.3768208
## 17  42 0.3722034
## 18  43 0.3677229
## 19  44 0.3633727
## 20  45 0.3591468
## 21  46 0.3550396
## 22  47 0.3510456
## 23  48 0.3471600
## 24  49 0.3433780
## 25  50 0.3396952
## 26  51 0.3361074
## 27  52 0.3326108
## 28  53 0.3292016
## 29  54 0.3258763
## 30  55 0.3226317
## 31  56 0.3194646
## 32  57 0.3163720
## 33  58 0.3133511
## 34  59 0.3103994
## 35  60 0.3075142
## 36  61 0.3046931
## 37  62 0.3019339
## 38  63 0.2992344
## 39  64 0.2965926
## 40  65 0.2940064
## 41  66 0.2914740
## 42  67 0.2889936
## 43  68 0.2865635
## 44  69 0.2841821
## 45  70 0.2818478
## 46  71 0.2795591
## 47  72 0.2773146
## 48  73 0.2751129
## 49  74 0.2729527
## 50  75 0.2708328
## 51  76 0.2687520
## 52  77 0.2667090
## 53  78 0.2647029
## 54  79 0.2627325
## 55  80 0.2607969
## 56  81 0.2588950
## 57  82 0.2570259
## 58  83 0.2551888
## 59  84 0.2533826
## 60  85 0.2516067
## 61  86 0.2498601
## 62  87 0.2481422
## 63  88 0.2464521
## 64  89 0.2447892
## 65  90 0.2431527
## 66  91 0.2415419
## 67  92 0.2399563
## 68  93 0.2383952
## 69  94 0.2368580
## 70  95 0.2353441
## 71  96 0.2338529
## 72  97 0.2323839
## 73  98 0.2309366
## 74  99 0.2295105
## 75 100 0.2281050

We could get substantial evidence for H0 with ~51 children per condition.

First TWo session Only

Frequentist

round(get_coeffs(training2.mod, c("cond.ct")),3)

##         Estimate Std. Error z value Pr(>|z|)
## cond.ct    0.001      0.124   0.009    0.993

Bayes Factor

meanBF = summary(training2.mod)$coefficients["cond.ct", "Estimate"]
seBF = summary(training2.mod)$coefficients["cond.ct", "Std. Error"]
h1mean = train.pred3.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.623753
## 
## $Likelihoodnull
## [1] 3.227875
## 
## $BayesFactor
## [1] 0.5030409

Don’t yet have substantial evidence for either H1 (that one condition shows greater accuracy than the other) or H0 (that there is no difference between conditions)

How much big a sample might we need?

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 26, max = 100))

##      x         y
## 1   26 0.4957468
## 2   27 0.4887651
## 3   28 0.4820741
## 4   29 0.4756540
## 5   30 0.4694871
## 6   31 0.4635572
## 7   32 0.4578494
## 8   33 0.4523503
## 9   34 0.4470475
## 10  35 0.4419294
## 11  36 0.4369856
## 12  37 0.4322066
## 13  38 0.4275832
## 14  39 0.4231072
## 15  40 0.4187710
## 16  41 0.4145674
## 17  42 0.4104899
## 18  43 0.4065321
## 19  44 0.4026884
## 20  45 0.3989535
## 21  46 0.3953221
## 22  47 0.3917898
## 23  48 0.3883520
## 24  49 0.3850046
## 25  50 0.3817436
## 26  51 0.3785656
## 27  52 0.3754669
## 28  53 0.3724443
## 29  54 0.3694949
## 30  55 0.3666156
## 31  56 0.3638037
## 32  57 0.3610567
## 33  58 0.3583720
## 34  59 0.3557474
## 35  60 0.3531807
## 36  61 0.3506697
## 37  62 0.3482125
## 38  63 0.3458072
## 39  64 0.3434519
## 40  65 0.3411451
## 41  66 0.3388849
## 42  67 0.3366699
## 43  68 0.3344986
## 44  69 0.3323696
## 45  70 0.3302815
## 46  71 0.3282330
## 47  72 0.3262228
## 48  73 0.3242499
## 49  74 0.3223131
## 50  75 0.3204112
## 51  76 0.3185433
## 52  77 0.3167083
## 53  78 0.3149053
## 54  79 0.3131334
## 55  80 0.3113916
## 56  81 0.3096792
## 57  82 0.3079953
## 58  83 0.3063392
## 59  84 0.3047100
## 60  85 0.3031071
## 61  86 0.3015298
## 62  87 0.2999773
## 63  88 0.2984491
## 64  89 0.2969445
## 65  90 0.2954629
## 66  91 0.2940038
## 67  92 0.2925665
## 68  93 0.2911506
## 69  94 0.2897555
## 70  95 0.2883807
## 71  96 0.2870257
## 72  97 0.2856900
## 73  98 0.2843733
## 74  99 0.2830750
## 75 100 0.2817947

Could get substantial evidence for H0 with about 67 per condition

Prediction 3: Summary

There is no evidence to support the theory that accuracy is higher in the Pictures + Diacritics condition, but also no evidence for null. Would need about 51 participants (all FOUR sessions) / 67 participants (first Two sessions only) per condition to get this.

Comment: It may not be feasible/interesting to seek for evidence here, given the interaction with session (see Prediction 4).

Prediction 4: Children will show greater improvement in one of the conditions

All FOUR sessions

Frequentist

round(get_coeffs(training4.mod, c("cond.ct:session.ct")),3)

##                    Estimate Std. Error z value Pr(>|z|)
## cond.ct:session.ct   -0.057      0.053  -1.075    0.282

Bayes Factor

Here we use the value from session to inform H1, as described in the plan. In the most extreme case, there would be improvement in one condition, say conditon A (so positive estimate for session in condition A) and in the other condition B they would not improve (i.e. estimate for session for that condition would be 0). (Note- we assume that they would not get worse). In that case - the maximum case - overall (average) esitmate for session would be equal to approximate half the estimate for session for just condition A, while the value of condition by session would be equal to the value of condition A. Thus in the maximum case, the value of condition by session is approximate twice the value for session. Therefore, for the theory that they improve more in one condition than the other, we set sd= estimate of session (so that twice this is the maximum).

In the plan we considered two theories: (i) They might learn more in the Pictures + Diacritics condition due to it being easier, or (ii) they might learn more in the Pictures-Only condition due to diacritics being hard for children to use. Below we test each of these as different H1s (so one-tailed).

H1= learning more in Pictures+Diactricis condition (note: our coding here means that the negative values is (correctly) in the opposite direction to the prediction)

meanBF = summary(training4.mod)$coefficients["cond.ct:session.ct", "Estimate"]
seBF = summary(training4.mod)$coefficients["cond.ct:session.ct", "Std. Error"]
h1mean = summary(training4.mod)$coefficients["session.ct", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.866629
## 
## $Likelihoodnull
## [1] 4.246048
## 
## $BayesFactor
## [1] 0.2041025

We have substantial evidence that they don’t learn more in the picture + diacritics condition.

H1= learning more in Pictures+only condition (note: our coding here means that value of the coefficient needs to be *-1 since the value is in the same difrection as the prediction)

meanBF = summary(training4.mod)$coefficients["cond.ct:session.ct", "Estimate"]*-1
seBF = summary(training4.mod)$coefficients["cond.ct:session.ct", "Std. Error"]
h1mean = summary(training4.mod)$coefficients["session.ct", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 4.497152
## 
## $Likelihoodnull
## [1] 4.246048
## 
## $BayesFactor
## [1] 1.059138

We don’t have evidence supporting either H1 or H0 that they are better in the Pictures-Only condition.

How big a sample might we need?

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N=25, min = 25, max=100))

##      x        y
## 1   25 1.059138
## 2   26 1.072254
## 3   27 1.085817
## 4   28 1.099822
## 5   29 1.114267
## 6   30 1.129149
## 7   31 1.144466
## 8   32 1.160220
## 9   33 1.176409
## 10  34 1.193036
## 11  35 1.210102
## 12  36 1.227610
## 13  37 1.245562
## 14  38 1.263963
## 15  39 1.282816
## 16  40 1.302126
## 17  41 1.321898
## 18  42 1.342138
## 19  43 1.362851
## 20  44 1.384043
## 21  45 1.405722
## 22  46 1.427893
## 23  47 1.450565
## 24  48 1.473744
## 25  49 1.497439
## 26  50 1.521658
## 27  51 1.546410
## 28  52 1.571703
## 29  53 1.597547
## 30  54 1.623951
## 31  55 1.650925
## 32  56 1.678479
## 33  57 1.706623
## 34  58 1.735369
## 35  59 1.764727
## 36  60 1.794708
## 37  61 1.825325
## 38  62 1.856589
## 39  63 1.888513
## 40  64 1.921109
## 41  65 1.954389
## 42  66 1.988369
## 43  67 2.023060
## 44  68 2.058476
## 45  69 2.094633
## 46  70 2.131545
## 47  71 2.169226
## 48  72 2.207692
## 49  73 2.246958
## 50  74 2.287041
## 51  75 2.327957
## 52  76 2.369723
## 53  77 2.412356
## 54  78 2.455873
## 55  79 2.500293
## 56  80 2.545633
## 57  81 2.591913
## 58  82 2.639152
## 59  83 2.687369
## 60  84 2.736585
## 61  85 2.786819
## 62  86 2.838094
## 63  87 2.890429
## 64  88 2.943849
## 65  89 2.998374
## 66  90 3.054027
## 67  91 3.110832
## 68  92 3.168814
## 69  93 3.227996
## 70  94 3.288403
## 71  95 3.350061
## 72  96 3.412996
## 73  97 3.477234
## 74  98 3.542804
## 75  99 3.609732
## 76 100 3.678047

Estimates we would need 90 per condition to have substnatial evidence for greater effect of session in the pictures-only condition.

First TWO sessions only

Frequentist

round(get_coeffs(training2.mod, c("cond.ct:session.ct")),3)

##                    Estimate Std. Error z value Pr(>|z|)
## cond.ct:session.ct   -0.206      0.105   -1.97    0.049

Bayes factor:

Note: Here we use the value from session to inform H1, as described in the plan. The idea is that one of the conditions might have a much larger effect of session. In contrast to plan I have set the value to be be HALF the value of session (see comments above).

In the plan we considered two theories: (i) They might learn more in the Pictures + Diacritics condition due to it being easier, or (ii) they might learn more in the just Pictures-Only condition due to diacritics being hard for children to use. Below we test each of these as different H1s (so one-tailed).

H1= learning more in Pictures+Diactricis condition (note: our coding here means that the negative values is (correctly) in the opposite direction to the prediction)

meanBF = summary(training2.mod)$coefficients["cond.ct:session.ct", "Estimate"]
seBF = summary(training2.mod)$coefficients["cond.ct:session.ct", "Std. Error"]
h1mean = summary(training2.mod)$coefficients["session.ct", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.07762045
## 
## $Likelihoodnull
## [1] 0.5481327
## 
## $BayesFactor
## [1] 0.1416089

So we have substantial evidence for the null for a theory that says you will learn more in the Pictures + Diacritics condition

H1= learning more in Pictures+only condition (note: our coding here means that the negative values is (correctly) in the opposite direction to the prediction)

meanBF = summary(training2.mod)$coefficients["cond.ct:session.ct", "Estimate"]*-1
seBF = summary(training2.mod)$coefficients["cond.ct:session.ct", "Std. Error"]
h1mean = summary(training2.mod)$coefficients["session.ct", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.1327
## 
## $Likelihoodnull
## [1] 0.5481327
## 
## $BayesFactor
## [1] 3.890846

Here, we also have substantial evidence for the theory that they learn more in the Pictures-only condition.

Prediction 4 Summary:

We tested both the theory that there is more improvement in the Pictures-Only condition VERSUS H0, AND the theory that there is more improvement in the Pictures + Diacritics condition versus H0:

Across the four sessions: - Theory of more improvement in Pictures + Diacritics: substantial evidence for the NULL - Theory of more improvement in Pictures-Only: no evidence either direction (sample analyses say that we could get this with about 90 participants per condition)

Across the first two sessions: - Theory of more improvement in Pictures + Diacritics: substantial evidence for the NULL - Theory of more improvement in Pictures-Only: substantial evidence for the theory

Additional question following up Prediction 4:

Above we found evidence for stronger learning in the Pictures-only compared to the Pictures + Diacritics condition when considering the change in performance from session 1 to session 2, but not when considering across all four sessions. Do we have evidence that the equivalent pattern is absent when considering (i) s2->s3 (ii) s3-> s4

Frequentist

Note that the models here are different (they have different number of dfs) from the original training4.mod above which had session coded as a numeric factor, rather than as contrasts.

training4$sessionb = as.factor(training4$session)
training4$sessionb = revalue(training4$sessionb, c("1"="s1","2"="s2","3"="s3","4"="s4"))

training4 = lizContrasts4(training4, training4$sessionb, "s1")
training4.modC_s1ref = glmer (acc ~
                       + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s3 +s1_VERSUS_s4) 
                        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
                       + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
                         + (session.ct|pt_N)
                         + (session.ct* cond.ct||item),
                       data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

training4 = lizContrasts4(training4, training4$sessionb, "s2")
training4.modC_s2ref = glmer (acc ~
                       + cond.ct * (s2_VERSUS_s1 + s2_VERSUS_s3 +s2_VERSUS_s4) 
                        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
                       + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
                        + (session.ct|pt_N)
                        + (session.ct* cond.ct||item),
                       data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))


training4 = lizContrasts4(training4, training4$sessionb, "s3")
training4.modC_s3ref = glmer (acc ~
                       + cond.ct * (s3_VERSUS_s1 + s3_VERSUS_s2 +s3_VERSUS_s4) 
                        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)
                       + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
                        + (session.ct|pt_N)
                        + (session.ct* cond.ct||item),
                       data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(training4.modC_s1ref , training4.modC_s2ref )

## Data: training4
## Models:
## training4.modC_s1ref: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s3 + s1_VERSUS_s4) + 
## training4.modC_s1ref:     (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + 
## training4.modC_s1ref:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## training4.modC_s1ref:     (session.ct | pt_N) + (session.ct * cond.ct || item)
## training4.modC_s2ref: acc ~ +cond.ct * (s2_VERSUS_s1 + s2_VERSUS_s3 + s2_VERSUS_s4) + 
## training4.modC_s2ref:     (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + 
## training4.modC_s2ref:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## training4.modC_s2ref:     (session.ct | pt_N) + (session.ct * cond.ct || item)
##                      Df   AIC   BIC logLik deviance Chisq Chi Df
## training4.modC_s1ref 23 25500 25682 -12727    25454             
## training4.modC_s2ref 23 25500 25682 -12727    25454     0      0
##                      Pr(>Chisq)
## training4.modC_s1ref           
## training4.modC_s2ref          1

anova(training4.modC_s3ref , training4.modC_s1ref )

## Data: training4
## Models:
## training4.modC_s3ref: acc ~ +cond.ct * (s3_VERSUS_s1 + s3_VERSUS_s2 + s3_VERSUS_s4) + 
## training4.modC_s3ref:     (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + 
## training4.modC_s3ref:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## training4.modC_s3ref:     (session.ct | pt_N) + (session.ct * cond.ct || item)
## training4.modC_s1ref: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s3 + s1_VERSUS_s4) + 
## training4.modC_s1ref:     (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4) + (T12_VERSUS_T13 + 
## training4.modC_s1ref:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## training4.modC_s1ref:     (session.ct | pt_N) + (session.ct * cond.ct || item)
##                      Df   AIC   BIC logLik deviance Chisq Chi Df
## training4.modC_s3ref 23 25500 25682 -12727    25454             
## training4.modC_s1ref 23 25500 25682 -12727    25454     0      0
##                      Pr(>Chisq)    
## training4.modC_s3ref               
## training4.modC_s1ref  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rbind(round(get_coeffs(training4.modC_s1ref, c("cond.ct:s1_VERSUS_s2")),3), round(get_coeffs(training4.modC_s2ref, c("cond.ct:s2_VERSUS_s3")),3), round(get_coeffs(training4.modC_s3ref, c("cond.ct:s3_VERSUS_s4")),3))

##                      Estimate Std. Error z value Pr(>|z|)
## cond.ct:s1_VERSUS_s2   -0.199      0.094  -2.111    0.035
## cond.ct:s2_VERSUS_s3    0.027      0.096   0.277    0.782
## cond.ct:s3_VERSUS_s4   -0.017      0.097  -0.179    0.858

Bayes factor

In each case, use interaction between condition and the contrast between session 1 and session 2 to inform h1 (we are testing the theory that we see an equivalent interaction between the other sessions) - we set sd of theory as equal to this value.

Note- the value for cond.ct:s3_VERSUS_s4 is negative and the formula doesn’t allow this. Therefore times all values by *1 s2->s3

meanBF = summary(training4.modC_s2ref)$coefficients["cond.ct:s2_VERSUS_s3", "Estimate"]*-1
seBF = summary(training4.modC_s2ref)$coefficients["cond.ct:s2_VERSUS_s3", "Std. Error"]
h1mean = summary(training4.modC_s1ref)$coefficients["cond.ct:s1_VERSUS_s2", "Estimate"]*-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.428998
## 
## $Likelihoodnull
## [1] 4.001744
## 
## $BayesFactor
## [1] 0.3570939

Don’t have substantial evidence for either h1 or h0.

Estimate sample size needed:

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 25, max = 100))

##      x         y
## 1   25 0.3570939
## 2   26 0.3499247
## 3   27 0.3431025
## 4   28 0.3366006
## 5   29 0.3303952
## 6   30 0.3244648
## 7   31 0.3187901
## 8   32 0.3133535
## 9   33 0.3081391
## 10  34 0.3031326
## 11  35 0.2983207
## 12  36 0.2936914
## 13  37 0.2892337
## 14  38 0.2849374
## 15  39 0.2807932
## 16  40 0.2767925
## 17  41 0.2729275
## 18  42 0.2691908
## 19  43 0.2655755
## 20  44 0.2620755
## 21  45 0.2586848
## 22  46 0.2553980
## 23  47 0.2522100
## 24  48 0.2491161
## 25  49 0.2461118
## 26  50 0.2431930
## 27  51 0.2403558
## 28  52 0.2375966
## 29  53 0.2349119
## 30  54 0.2322986
## 31  55 0.2297536
## 32  56 0.2272740
## 33  57 0.2248572
## 34  58 0.2225006
## 35  59 0.2202019
## 36  60 0.2179587
## 37  61 0.2157689
## 38  62 0.2136306
## 39  63 0.2115416
## 40  64 0.2095004
## 41  65 0.2075050
## 42  66 0.2055538
## 43  67 0.2036454
## 44  68 0.2017781
## 45  69 0.1999506
## 46  70 0.1981615
## 47  71 0.1964096
## 48  72 0.1946935
## 49  73 0.1930122
## 50  74 0.1913645
## 51  75 0.1897493
## 52  76 0.1881656
## 53  77 0.1866124
## 54  78 0.1850888
## 55  79 0.1835939
## 56  80 0.1821268
## 57  81 0.1806867
## 58  82 0.1792728
## 59  83 0.1778844
## 60  84 0.1765207
## 61  85 0.1751809
## 62  86 0.1738645
## 63  87 0.1725708
## 64  88 0.1712991
## 65  89 0.1700490
## 66  90 0.1688196
## 67  91 0.1676107
## 68  92 0.1664215
## 69  93 0.1652515
## 70  94 0.1641004
## 71  95 0.1629675
## 72  96 0.1618525
## 73  97 0.1607549
## 74  98 0.1596742
## 75  99 0.1586100
## 76 100 0.1575620

Estiamtes we We could have substantial evidence for H0 with 29 participants in each condition (another 7 total)

s3->s4

meanBF = summary(training4.modC_s3ref)$coefficients["cond.ct:s3_VERSUS_s4", "Estimate"]*-1
seBF = summary(training4.modC_s3ref)$coefficients["cond.ct:s3_VERSUS_s4", "Std. Error"]
h1mean = summary(training4.modC_s1ref)$coefficients["cond.ct:s1_VERSUS_s2", "Estimate"]*-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.014423
## 
## $Likelihoodnull
## [1] 4.037084
## 
## $BayesFactor
## [1] 0.4989797

Estimate sample size needed:

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 25, max = 100))

##      x         y
## 1   25 0.4989797
## 2   26 0.4926929
## 3   27 0.4866946
## 4   28 0.4809638
## 5   29 0.4754814
## 6   30 0.4702302
## 7   31 0.4651947
## 8   32 0.4603607
## 9   33 0.4557153
## 10  34 0.4512467
## 11  35 0.4469444
## 12  36 0.4427983
## 13  37 0.4387994
## 14  38 0.4349394
## 15  39 0.4312106
## 16  40 0.4276059
## 17  41 0.4241187
## 18  42 0.4207429
## 19  43 0.4174729
## 20  44 0.4143033
## 21  45 0.4112293
## 22  46 0.4082462
## 23  47 0.4053498
## 24  48 0.4025362
## 25  49 0.3998014
## 26  50 0.3971421
## 27  51 0.3945549
## 28  52 0.3920368
## 29  53 0.3895847
## 30  54 0.3871960
## 31  55 0.3848681
## 32  56 0.3825985
## 33  57 0.3803849
## 34  58 0.3782252
## 35  59 0.3761172
## 36  60 0.3740590
## 37  61 0.3720487
## 38  62 0.3700847
## 39  63 0.3681651
## 40  64 0.3662885
## 41  65 0.3644533
## 42  66 0.3626581
## 43  67 0.3609016
## 44  68 0.3591823
## 45  69 0.3574991
## 46  70 0.3558508
## 47  71 0.3542363
## 48  72 0.3526544
## 49  73 0.3511042
## 50  74 0.3495846
## 51  75 0.3480948
## 52  76 0.3466337
## 53  77 0.3452006
## 54  78 0.3437946
## 55  79 0.3424149
## 56  80 0.3410608
## 57  81 0.3397314
## 58  82 0.3384262
## 59  83 0.3371444
## 60  84 0.3358854
## 61  85 0.3346485
## 62  86 0.3334332
## 63  87 0.3322389
## 64  88 0.3310650
## 65  89 0.3299110
## 66  90 0.3287763
## 67  91 0.3276605
## 68  92 0.3265631
## 69  93 0.3254836
## 70  94 0.3244215
## 71  95 0.3233765
## 72  96 0.3223482
## 73  97 0.3213360
## 74  98 0.3203397
## 75  99 0.3193588
## 76 100 0.3183930

Estimates would need 84 participants per condition (117 more in total).

Additional Question Summary

Additional analyses (unplanned) explored whether the interaction between condition and session which holds from s1->s2 is also relevant from s2->s3 and s3->s4. Currently do not have substantial evidence for either in either direction; for s2->s3 estimate could have substantial evidence with about anoher 7 participants.

2AFC Pictures & 2AFC Diacritics

Inspect Data

str(afc)

## 'data.frame':    4896 obs. of  20 variables:
##  $ sort           : int  49 50 51 52 53 54 55 56 57 58 ...
##  $ half           : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ pt_N           : Factor w/ 51 levels "Participant 10",..: 10 10 10 10 10 10 10 10 10 10 ...
##  $ program        : Factor w/ 8 levels "S2_ver1_P.cond\\",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ task           : Factor w/ 2 levels "[pictest_pic$]",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ soundfile      : Factor w/ 24 levels "chi_t2_f2a.wav",..: 22 8 10 23 17 7 22 7 18 9 ...
##  $ picture1       : Factor w/ 28 levels "chi_t2_spoon.jpg",..: 26 8 10 27 17 7 26 7 17 9 ...
##  $ picture2       : Factor w/ 28 levels "chi_t2_spoon.jpg",..: 25 7 9 28 18 8 25 8 18 10 ...
##  $ picture_correct: Factor w/ 2 levels "p1","p2": 1 1 1 1 1 1 1 1 2 1 ...
##  $ tone_contrast  : Factor w/ 6 levels "T12","T13","T14",..: 3 2 1 6 5 2 3 2 5 1 ...
##  $ speaker        : Factor w/ 1 level "nv1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ session        : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ DEL            : logi  NA NA NA NA NA NA ...
##  $ correct        : Factor w/ 2 levels "correct","incorrect": 2 2 1 2 1 2 1 2 2 1 ...
##  $ side_clicked   : Factor w/ 2 levels "clickedleft",..: 2 2 1 2 1 2 1 2 1 1 ...
##  $ RT             : num  1856 1261 1718 3428 3958 ...
##  $ acc            : int  0 0 1 0 1 0 1 0 0 1 ...
##  $ version        : Factor w/ 2 levels "ver1","ver2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ cond           : Factor w/ 2 levels "P","P+T": 2 2 2 2 2 2 2 2 2 2 ...
##  $ item           : Factor w/ 24 levels "chit2","chit3",..: 22 8 10 23 17 7 22 7 18 9 ...

afc$session = as.factor(afc$session)
afc$half = as.factor(afc$half)

Get Means

a = summarySEwithin(aggregate(acc ~ pt_N + task + cond + session, FUN = mean, data = afc), measurevar="acc", betweenvars = c("task", "cond"), withinvars = c("session"))
kable(a, digits = 2)

task	cond	session	N	acc	acc_norm	sd	se	ci
[pictest_pic$]	P	2	26	0.58	0.57	0.19	0.04	0.08
[pictest_pic$]	P	4	26	0.61	0.60	0.20	0.04	0.08
[pictest_pic$]	P+T	2	25	0.57	0.58	0.25	0.05	0.11
[pictest_pic$]	P+T	4	25	0.58	0.59	0.19	0.04	0.08
[pictest_tone$]	P+T	2	25	0.56	0.56	0.21	0.04	0.09
[pictest_tone$]	P+T	4	25	0.60	0.60	0.21	0.04	0.09
## Get plot

p2 = ggplot(a, aes( x = cond, y = acc, fill = session))
p2 = p2 + facet_wrap(~ task, ncol = 3)
p2 = p2 + geom_bar(stat = "identity", position = "dodge", colour = "blue", size = .2)
p2 = p2 + geom_errorbar(aes(ymin = acc-ci, ymax = acc+ci), width = .2, position = position_dodge(.9))
p2 = p2 + xlab("Condition") + ylab("Mean accuracy")
p2 = p2 + theme_bw()
p2 = p2 + coord_cartesian(ylim = c(0.33, 0.8))
p2 = p2 +  geom_hline(aes(yintercept = .5) )
p2

AFC Pictures

Main LME

Comment on control variables: We also include version (1 vs. 2), and tone contrast (1-2, 1-3, 1-4, 2-3, 2-4, 3-4) as fixed effects. Both are expected to contribute to the model and are thus included but neither is of specific interest.

afcPIC <- droplevels(subset(afc, task == "[pictest_pic$]"))
afcPIC$session = as.factor(afcPIC$session)
afcPIC <- lizCenter(afcPIC, list("cond","version", "session"))
afcPIC = lizContrasts6(afcPIC, afcPIC$tone_contrast, "T12")

afcPIC.mod1 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        ,
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

# converges: now find slope structre for items

afcPIC.mod2 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (session.ct*cond.ct|item),
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00131149 (tol =
## 0.001, component 1)

# doesn't converge, remove corrleations between slopes

afcPIC.mod3 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (session.ct*cond.ct||item),
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

afcPIC.mod4 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (session.ct+cond.ct||item),
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod3 ,afcPIC.mod4)

## Data: afcPIC
## Models:
## afcPIC.mod4: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod4:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod4:     pt_N) + (session.ct + cond.ct || item)
## afcPIC.mod3: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod3:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod3:     pt_N) + (session.ct * cond.ct || item)
##             Df    AIC    BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## afcPIC.mod4 15 4883.9 4977.2  -2427   4853.9                        
## afcPIC.mod3 16 4885.9 4985.4  -2427   4853.9     0      1          1

# p>.2, accept simpler model


afcPIC.mod5 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (session.ct||item),
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))


afcPIC.mod6 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (cond.ct||item),
        
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod4,afcPIC.mod5)

## Data: afcPIC
## Models:
## afcPIC.mod5: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod5:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod5:     pt_N) + (session.ct || item)
## afcPIC.mod4: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod4:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod4:     pt_N) + (session.ct + cond.ct || item)
##             Df    AIC    BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## afcPIC.mod5 14 4881.9 4968.9  -2427   4853.9                        
## afcPIC.mod4 15 4883.9 4977.2  -2427   4853.9     0      1          1

anova(afcPIC.mod4,afcPIC.mod6)

## Data: afcPIC
## Models:
## afcPIC.mod6: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod6:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod6:     pt_N) + (cond.ct || item)
## afcPIC.mod4: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod4:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod4:     pt_N) + (session.ct + cond.ct || item)
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## afcPIC.mod6 14 4883.0 4970.0 -2427.5   4855.0                         
## afcPIC.mod4 15 4883.9 4977.2 -2427.0   4853.9 1.0669      1     0.3017

afcPIC.mod7 = glmer (acc ~
        + cond.ct * session.ct 
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)

        + (session.ct|pt_N)
        + (1|item),
        
data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod4,afcPIC.mod7)

## Data: afcPIC
## Models:
## afcPIC.mod7: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod7:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod7:     pt_N) + (1 | item)
## afcPIC.mod4: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod4:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod4:     pt_N) + (session.ct + cond.ct || item)
##             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## afcPIC.mod7 13 4881.0 4961.8 -2427.5   4855.0                         
## afcPIC.mod4 15 4883.9 4977.2 -2427.0   4853.9 1.0669      2     0.5866

afcPIC.mod = afcPIC.mod7

round(summary(afcPIC.mod)$coefficients,3)

##                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)           0.391      0.084   4.680    0.000
## cond.ct              -0.117      0.162  -0.720    0.472
## session.ct            0.093      0.085   1.101    0.271
## T12_VERSUS_T13       -0.058      0.120  -0.481    0.630
## T12_VERSUS_T14       -0.312      0.119  -2.617    0.009
## T12_VERSUS_T23        0.058      0.121   0.484    0.629
## T12_VERSUS_T24       -0.305      0.119  -2.558    0.011
## T12_VERSUS_T34        0.051      0.121   0.423    0.672
## cond.ct:session.ct   -0.087      0.174  -0.503    0.615

Prediction 1: Performance will be above chance in both conditions

Frequentist

As per plan, we look both at s2 and s4. Here we do it within one model, fitting separate intercepts for each session in each condition.

afcPIC.modA = glmer (acc ~
        + cond:session
        - 1
        - cond.ct
        - session.ct
        
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (1|item),
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod, afcPIC.modA)

## Data: afcPIC
## Models:
## afcPIC.mod: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod:     pt_N) + (1 | item)
## afcPIC.modA: acc ~ +cond:session - 1 - cond.ct - session.ct + (T12_VERSUS_T13 + 
## afcPIC.modA:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## afcPIC.modA:     (session.ct | pt_N) + (1 | item)
##             Df  AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## afcPIC.mod  13 4881 4961.8 -2427.5     4855                            
## afcPIC.modA 13 4881 4961.8 -2427.5     4855     0      0  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

round(get_coeffs(afcPIC.modA, c("condP:session2", "condP+T:session2","condP:session4","condP+T:session4")), 3)

##                  Estimate Std. Error z value Pr(>|z|)
## condP:session2      0.368      0.128   2.875    0.004
## condP+T:session2    0.295      0.144   2.051    0.040
## condP:session4      0.490      0.115   4.245    0.000
## condP+T:session4    0.330      0.132   2.507    0.012

Bayes Factor

For Bayesian analyses we inform our H1 using values from training data obtained in equivalent session. Below we get the values for training for s2 and s4.

Note: it isn’t clear here whether we should treat these values as an estimate, or a maximum. (For P+T condition, it feels like it could be a maximum, since we have removed a cue, however we don’t want to have different standards H1 for the two conditions). Therefore, we try with both.

training4.modD = glmer (acc ~
        -1
        + cond:sessionb
       
        + (tv1_VERSUS_tv2 + tv1_VERSUS_tv3 + tv1_VERSUS_tv4)  
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (cond.ct*session.ct||item)
        + (session.ct|pt_N),
        data = training4, family = binomial, control = glmerControl(optimizer = "bobyqa"))


afc.s4.pred1.condP.H1 = summary(training4.modD)$coefficients["condP:sessionbs4", "Estimate"]
afc.s4.pred1.condPT.H1 = summary(training4.modD)$coefficients["condP+T:sessionbs4", "Estimate"]
afc.s2.pred1.condP.H1 = summary(training4.modD)$coefficients["condP:sessionbs2", "Estimate"]
afc.s2.pred1.condPT.H1 = summary(training4.modD)$coefficients["condP+T:sessionbs2", "Estimate"]

Pictures-Only Condition, Session 4: Using value above as estimate of H1

meanBF = summary(afcPIC.modA)$coefficients["condP:session4", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP:session4", "Std. Error"]
h1mean = afc.s4.pred1.condP.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.9357538
## 
## $Likelihoodnull
## [1] 0.0004227347
## 
## $BayesFactor
## [1] 2213.573

Pictures-Only Condition, Session 4: Using value above as maximum of H1

meanBF = summary(afcPIC.modA)$coefficients["condP:session4", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP:session4", "Std. Error"]
h1mean = afc.s4.pred1.condP.H1/2 
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.8048479
## 
## $Likelihoodnull
## [1] 0.0004227347
## 
## $BayesFactor
## [1] 1903.908

Pictures-Only Condition, Session 2: Using value above as estimate of H1

meanBF = summary(afcPIC.modA)$coefficients["condP:session2", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP:session2", "Std. Error"]
h1mean = afc.s2.pred1.condP.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.217945
## 
## $Likelihoodnull
## [1] 0.0499082
## 
## $BayesFactor
## [1] 24.4037

Pictures-Only Condition, Session 2: Using value above as maximum of H1

meanBF = summary(afcPIC.modA)$coefficients["condP:session2", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP:session2", "Std. Error"]
h1mean = afc.s2.pred1.condP.H1/2 
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.166363
## 
## $Likelihoodnull
## [1] 0.0499082
## 
## $BayesFactor
## [1] 23.37018

Pictures and Diactricis Condition, Session 4: Using value above as estimate of H1

meanBF = summary(afcPIC.modA)$coefficients["condP+T:session4", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP+T:session4", "Std. Error"]
h1mean = afc.s4.pred1.condPT.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.202357
## 
## $Likelihoodnull
## [1] 0.1310912
## 
## $BayesFactor
## [1] 9.171914

Pictures and Diactricis Condition, Session 4: Using value above as maximum of H1

meanBF = summary(afcPIC.modA)$coefficients["condP+T:session4", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP+T:session4", "Std. Error"]
h1mean = afc.s4.pred1.condPT.H1/2 
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.434674
## 
## $Likelihoodnull
## [1] 0.1310912
## 
## $BayesFactor
## [1] 10.94409

Pictures and Diactricis Condition, Session 2: Using value above as estimate of H1

meanBF = summary(afcPIC.modA)$coefficients["condP+T:session2", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP+T:session2", "Std. Error"]
h1mean = afc.s2.pred1.condPT.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.457767
## 
## $Likelihoodnull
## [1] 0.3382534
## 
## $BayesFactor
## [1] 4.30969

Pictures and Diactricis Condition, Session 2: Using value above as maximum of H1

meanBF = summary(afcPIC.modA)$coefficients["condP+T:session2", "Estimate"]
seBF = summary(afcPIC.modA)$coefficients["condP+T:session2", "Std. Error"]
h1mean = afc.s2.pred1.condPT.H1/2 
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.483376
## 
## $Likelihoodnull
## [1] 0.3382534
## 
## $BayesFactor
## [1] 4.385399

Prediction 1 Summary

For Pictures only, at both post-tests separately we have strong evidence for H1 over H0 of chance. For Pictures + Diactricis, at session 2 post test we have substantial evidence that performance is above chance and at post test 4 we have strong evidence of H1 over H0 of chance.

Prediction 2: Performance will be higher in the pictures-only condition

Frequentist

Again look at both s2 and s4, so refit model with effect of session for each:

round(get_coeffs(afcPIC.mod, c("cond.ct")), 3)

##         Estimate Std. Error z value Pr(>|z|)
## cond.ct   -0.117      0.162   -0.72    0.472

afcPIC.modB = glmer (acc ~
        + cond.ct:session
        + session.ct
       
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (1|item)
        + (session.ct|pt_N),
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod, afcPIC.modB)

## Data: afcPIC
## Models:
## afcPIC.mod: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod:     pt_N) + (1 | item)
## afcPIC.modB: acc ~ +cond.ct:session + session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.modB:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (1 | 
## afcPIC.modB:     item) + (session.ct | pt_N)
##             Df  AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## afcPIC.mod  13 4881 4961.8 -2427.5     4855                        
## afcPIC.modB 13 4881 4961.8 -2427.5     4855     0      0          1

round(get_coeffs(afcPIC.modB, c("cond.ct:session2", "cond.ct:session4")), 3)

##                  Estimate Std. Error z value Pr(>|z|)
## cond.ct:session2   -0.073      0.192  -0.379    0.704
## cond.ct:session4   -0.160      0.175  -0.917    0.359

Bayes Factor

As per the plan, using relevant intercept (i.e. from that level of session) to inform H1. The logic is that the maximum value of condition is seen if one condition - say conditionA - shows all the difference from chance and the other condition (condition B) is at chance (assuming that we won’t see below chance behavior). In this situation, the intercept for condition B (in this session) will be 0 (ie log odds 50); thus the the average intercept (in this session) will be approximately half as big the intercept for condition A (in this session), and the effect of condition (in this session) will be equal to the intercept of conditionA. This means that in the maximum case, the effect of condition (in this session) will be equal to two times the intercept (in this session). Thus if we set sd for theory to be equal to the intercept, this will capture the fact that twice this value is the maximum possible value.

In order to compute Bayes factor we need to obtain intercept for each level of session. We do this in the following model:

afcPIC.modC = glmer (acc ~
        - 1
        + session
        + cond.ct
        + cond.ct:session.ct
        + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
        + (session.ct|pt_N)
        + (1|item),
        data = afcPIC, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcPIC.mod, afcPIC.modC)

## Data: afcPIC
## Models:
## afcPIC.mod: acc ~ +cond.ct * session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + 
## afcPIC.mod:     T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | 
## afcPIC.mod:     pt_N) + (1 | item)
## afcPIC.modC: acc ~ -1 + session + cond.ct + cond.ct:session.ct + (T12_VERSUS_T13 + 
## afcPIC.modC:     T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34) + 
## afcPIC.modC:     (session.ct | pt_N) + (1 | item)
##             Df  AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## afcPIC.mod  13 4881 4961.8 -2427.5     4855                            
## afcPIC.modC 13 4881 4961.8 -2427.5     4855     0      0  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(afcPIC.modC)$coefficients

##                       Estimate Std. Error    z value     Pr(>|z|)
## session2            0.34443423 0.09833126  3.5027948 4.604039e-04
## session4            0.43791158 0.08894161  4.9235852 8.497286e-07
## cond.ct            -0.11666293 0.16203486 -0.7199866 4.715332e-01
## T12_VERSUS_T13     -0.05776361 0.12006646 -0.4810970 6.304476e-01
## T12_VERSUS_T14     -0.31209225 0.11924753 -2.6171799 8.865959e-03
## T12_VERSUS_T23      0.05840671 0.12072710  0.4837912 6.285341e-01
## T12_VERSUS_T24     -0.30511179 0.11925772 -2.5584237 1.051479e-02
## T12_VERSUS_T34      0.05106777 0.12067874  0.4231713 6.721703e-01
## cond.ct:session.ct -0.08734335 0.17369005 -0.5028691 6.150563e-01

round(get_coeffs(afcPIC.modC, c("session2","session4")), 3)

##          Estimate Std. Error z value Pr(>|z|)
## session2    0.344      0.098   3.503        0
## session4    0.438      0.089   4.924        0

Session 4:

(Note that performance is numerically higher in the pictures-only condition, so the estimate is in the direction of the theory and we *-1)

meanBF = summary(afcPIC.modB)$coefficients["cond.ct:session4", "Estimate"] *-1
seBF = summary(afcPIC.modB)$coefficients["cond.ct:session4", "Std. Error"]
h1mean =  summary(afcPIC.modC)$coefficients["session4", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.285188
## 
## $Likelihoodnull
## [1] 1.498507
## 
## $BayesFactor
## [1] 0.8576459

The evidence is ambiguous. How big a sample might we need?

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 25, max = 100))

##      x         y
## 1   25 0.8576459
## 2   26 0.8628893
## 3   27 0.8684222
## 4   28 0.8742331
## 5   29 0.8803119
## 6   30 0.8866495
## 7   31 0.8932383
## 8   32 0.9000714
## 9   33 0.9071427
## 10  34 0.9144471
## 11  35 0.9219799
## 12  36 0.9297372
## 13  37 0.9377156
## 14  38 0.9459120
## 15  39 0.9543240
## 16  40 0.9629494
## 17  41 0.9717864
## 18  42 0.9808337
## 19  43 0.9900901
## 20  44 0.9995548
## 21  45 1.0092271
## 22  46 1.0191067
## 23  47 1.0291934
## 24  48 1.0394874
## 25  49 1.0499888
## 26  50 1.0606982
## 27  51 1.0716160
## 28  52 1.0827431
## 29  53 1.0940804
## 30  54 1.1056289
## 31  55 1.1173897
## 32  56 1.1293643
## 33  57 1.1415540
## 34  58 1.1539603
## 35  59 1.1665850
## 36  60 1.1794298
## 37  61 1.1924965
## 38  62 1.2057872
## 39  63 1.2193038
## 40  64 1.2330486
## 41  65 1.2470238
## 42  66 1.2612317
## 43  67 1.2756748
## 44  68 1.2903556
## 45  69 1.3052767
## 46  70 1.3204408
## 47  71 1.3358505
## 48  72 1.3515089
## 49  73 1.3674188
## 50  74 1.3835832
## 51  75 1.4000052
## 52  76 1.4166879
## 53  77 1.4336347
## 54  78 1.4508488
## 55  79 1.4683335
## 56  80 1.4860925
## 57  81 1.5041292
## 58  82 1.5224473
## 59  83 1.5410505
## 60  84 1.5599425
## 61  85 1.5791272
## 62  86 1.5986085
## 63  87 1.6183906
## 64  88 1.6384773
## 65  89 1.6588730
## 66  90 1.6795819
## 67  91 1.7006083
## 68  92 1.7219567
## 69  93 1.7436314
## 70  94 1.7656372
## 71  95 1.7879786
## 72  96 1.8106604
## 73  97 1.8336875
## 74  98 1.8570647
## 75  99 1.8807971
## 76 100 1.9048897

Estimates that even with a sample of 100 participants per condition evidence would still be ambiguous

Session 2

meanBF = summary(afcPIC.modB)$coefficients["cond.ct:session2", "Estimate"] *-1
seBF = summary(afcPIC.modB)$coefficients["cond.ct:session2", "Std. Error"]
h1mean =  summary(afcPIC.modC)$coefficients["session2", "Estimate"]
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.248006
## 
## $Likelihoodnull
## [1] 1.929144
## 
## $BayesFactor
## [1] 0.6469221

The evidence is ambiguous. How big a sample might we need?

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 25, max = 100))

##      x         y
## 1   25 0.6469221
## 2   26 0.6421071
## 3   27 0.6375284
## 4   28 0.6331699
## 5   29 0.6290172
## 6   30 0.6250570
## 7   31 0.6212773
## 8   32 0.6176669
## 9   33 0.6142160
## 10  34 0.6109151
## 11  35 0.6077558
## 12  36 0.6047301
## 13  37 0.6018310
## 14  38 0.5990516
## 15  39 0.5963858
## 16  40 0.5938277
## 17  41 0.5913720
## 18  42 0.5890138
## 19  43 0.5867484
## 20  44 0.5845714
## 21  45 0.5824788
## 22  46 0.5804668
## 23  47 0.5785318
## 24  48 0.5766706
## 25  49 0.5748798
## 26  50 0.5731567
## 27  51 0.5714984
## 28  52 0.5699023
## 29  53 0.5683660
## 30  54 0.5668871
## 31  55 0.5654635
## 32  56 0.5640930
## 33  57 0.5627736
## 34  58 0.5615036
## 35  59 0.5602811
## 36  60 0.5591045
## 37  61 0.5579721
## 38  62 0.5568826
## 39  63 0.5558343
## 40  64 0.5548260
## 41  65 0.5538563
## 42  66 0.5529240
## 43  67 0.5520279
## 44  68 0.5511670
## 45  69 0.5503400
## 46  70 0.5495460
## 47  71 0.5487839
## 48  72 0.5480530
## 49  73 0.5473521
## 50  74 0.5466805
## 51  75 0.5460374
## 52  76 0.5454220
## 53  77 0.5448334
## 54  78 0.5442710
## 55  79 0.5437340
## 56  80 0.5432219
## 57  81 0.5427339
## 58  82 0.5422694
## 59  83 0.5418279
## 60  84 0.5414087
## 61  85 0.5410113
## 62  86 0.5406352
## 63  87 0.5402799
## 64  88 0.5399448
## 65  89 0.5396296
## 66  90 0.5393336
## 67  91 0.5390566
## 68  92 0.5387980
## 69  93 0.5385575
## 70  94 0.5383347
## 71  95 0.5381291
## 72  96 0.5379405
## 73  97 0.5377683
## 74  98 0.5376124
## 75  99 0.5374724
## 76 100 0.5373479

Suggests evidence would be Ambiguous, even with N=100

Prediciton 2 Summary

Evidence is ambiguous both for Session 2 and Session 4 and across sessions. The samples which would be required to look at this are infeasible.

AFC Diacritics Test: Pictures + Diacritics condition only

Main LME

afcDIACRITICS <- droplevels(subset(afc, task == "[pictest_tone$]"))
afcDIACRITICS$session = as.factor(afcDIACRITICS$session) 
afcDIACRITICS <- lizCenter(afcDIACRITICS, list("session", "version"))
afcDIACRITICS <- lizContrasts6(afcDIACRITICS, afcDIACRITICS$tone_contrast, "T12")

afcDIACRITICS.mod1 = glmer (acc ~
                 session.ct 
              + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)        
              + (session.ct|pt_N)
              ,
              data = afcDIACRITICS, family = binomial, control = glmerControl(optimizer = "bobyqa"))

# converges, now find slope structure for item
afcDIACRITICS.mod2 = glmer (acc ~
                 session.ct 
              + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)        
              + (session.ct|pt_N)
              + (session.ct|item),
              data = afcDIACRITICS, family = binomial, control = glmerControl(optimizer = "bobyqa"))

afcDIACRITICS.mod3 = glmer (acc ~
                 session.ct 
              + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)        
              + (session.ct|pt_N)
              + (session.ct||item),
              data = afcDIACRITICS, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcDIACRITICS.mod2,afcDIACRITICS.mod3)

## Data: afcDIACRITICS
## Models:
## afcDIACRITICS.mod3: acc ~ session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.mod3:     (session.ct || item)
## afcDIACRITICS.mod2: acc ~ session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.mod2:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.mod2:     (session.ct | item)
##                    Df    AIC    BIC  logLik deviance Chisq Chi Df
## afcDIACRITICS.mod3 12 1607.7 1668.8 -791.83   1583.7             
## afcDIACRITICS.mod2 13 1609.7 1675.8 -791.83   1583.7     0      1
##                    Pr(>Chisq)
## afcDIACRITICS.mod3           
## afcDIACRITICS.mod2          1

# p>.2 so move to simpler model

afcDIACRITICS.mod4 = glmer (acc ~
                 session.ct 
              + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)        
              + (session.ct|pt_N)
              + (1|item),
              data = afcDIACRITICS, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcDIACRITICS.mod3,afcDIACRITICS.mod2)

## Data: afcDIACRITICS
## Models:
## afcDIACRITICS.mod3: acc ~ session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.mod3:     (session.ct || item)
## afcDIACRITICS.mod2: acc ~ session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.mod2:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.mod2:     (session.ct | item)
##                    Df    AIC    BIC  logLik deviance Chisq Chi Df
## afcDIACRITICS.mod3 12 1607.7 1668.8 -791.83   1583.7             
## afcDIACRITICS.mod2 13 1609.7 1675.8 -791.83   1583.7     0      1
##                    Pr(>Chisq)
## afcDIACRITICS.mod3           
## afcDIACRITICS.mod2          1

# p>.2 so accept simpler model

afcDIACRITICS.mod= afcDIACRITICS.mod4
round(summary(afcDIACRITICS.mod)$coefficients,3)

##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)       0.343      0.124   2.763    0.006
## session.ct        0.187      0.122   1.531    0.126
## T12_VERSUS_T13    0.615      0.211   2.914    0.004
## T12_VERSUS_T14    0.258      0.208   1.244    0.214
## T12_VERSUS_T23    0.456      0.209   2.181    0.029
## T12_VERSUS_T24    0.346      0.208   1.660    0.097
## T12_VERSUS_T34    0.215      0.207   1.036    0.300

Prediction 1: They will be above chance in this condition

Frequentist

Fit separate intercepts for session2 and session4

afcDIACRITICS.modB = glmer (acc ~
              + session 
              - 1
              + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)        
              + (session.ct|pt_N)
              + (1|item),
              data = afcDIACRITICS, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(afcDIACRITICS.modB, afcDIACRITICS.mod)

## Data: afcDIACRITICS
## Models:
## afcDIACRITICS.modB: acc ~ +session - 1 + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.modB:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.modB:     (1 | item)
## afcDIACRITICS.mod: acc ~ session.ct + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## afcDIACRITICS.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + (session.ct | pt_N) + 
## afcDIACRITICS.mod:     (1 | item)
##                    Df    AIC    BIC  logLik deviance Chisq Chi Df
## afcDIACRITICS.modB 11 1605.7 1661.7 -791.83   1583.7             
## afcDIACRITICS.mod  11 1605.7 1661.7 -791.83   1583.7     0      0
##                    Pr(>Chisq)
## afcDIACRITICS.modB           
## afcDIACRITICS.mod           1

round(get_coeffs(afcDIACRITICS.modB, c("session2","session4")),3)

##          Estimate Std. Error z value Pr(>|z|)
## session2    0.250      0.135   1.851    0.064
## session4    0.437      0.142   3.078    0.002

Bayes factor

Bayesian analyses: use the same values as for equivalent prediction for 2AFC pictures. Here it is clearer that value from training should be a maxium, however for consistency we use it both as an estimate and a maximum.

Session 2, using value as an estimate:

meanBF = summary(afcDIACRITICS.modB)$coefficients["session2", "Estimate"]
seBF = summary(afcDIACRITICS.modB)$coefficients["session2", "Std. Error"]
h1mean = afc.s2.pred1.condPT.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.55473
## 
## $Likelihoodnull
## [1] 0.5338273
## 
## $BayesFactor
## [1] 2.912422

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min = 25, max = 100))

##      x          y
## 1   25   2.912422
## 2   26   3.072728
## 3   27   3.242996
## 4   28   3.423846
## 5   29   3.615942
## 6   30   3.819991
## 7   31   4.036747
## 8   32   4.267014
## 9   33   4.511648
## 10  34   4.771561
## 11  35   5.047728
## 12  36   5.341184
## 13  37   5.653035
## 14  38   5.984456
## 15  39   6.336704
## 16  40   6.711115
## 17  41   7.109113
## 18  42   7.532217
## 19  43   7.982044
## 20  44   8.460320
## 21  45   8.968881
## 22  46   9.509685
## 23  47  10.084822
## 24  48  10.696514
## 25  49  11.347134
## 26  50  12.039209
## 27  51  12.775434
## 28  52  13.558681
## 29  53  14.392010
## 30  54  15.278688
## 31  55  16.222192
## 32  56  17.226232
## 33  57  18.294765
## 34  58  19.432006
## 35  59  20.642451
## 36  60  21.930895
## 37  61  23.302450
## 38  62  24.762566
## 39  63  26.317056
## 40  64  27.972118
## 41  65  29.734363
## 42  66  31.610841
## 43  67  33.609068
## 44  68  35.737064
## 45  69  38.003381
## 46  70  40.417141
## 47  71  42.988072
## 48  72  45.726552
## 49  73  48.643650
## 50  74  51.751174
## 51  75  55.061721
## 52  76  58.588727
## 53  77  62.346527
## 54  78  66.350413
## 55  79  70.616698
## 56  80  75.162790
## 57  81  80.007259
## 58  82  85.169916
## 59  83  90.671903
## 60  84  96.535775
## 61  85 102.785600
## 62  86 109.447058
## 63  87 116.547553
## 64  88 124.116327
## 65  89 132.184583
## 66  90 140.785621
## 67  91 149.954976
## 68  92 159.730572
## 69  93 170.152880
## 70  94 181.265094
## 71  95 193.113311
## 72  96 205.746733
## 73  97 219.217869
## 74  98 233.582765
## 75  99 248.901244
## 76 100 265.237155

Estimates could have substnatial evidence with just one more participants

Session 2, using value as a maximum:

meanBF = summary(afcDIACRITICS.modB)$coefficients["session2", "Estimate"]
seBF = summary(afcDIACRITICS.modB)$coefficients["session2", "Std. Error"]
h1mean = afc.s2.pred1.condPT.H1/2
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.806042
## 
## $Likelihoodnull
## [1] 0.5338273
## 
## $BayesFactor
## [1] 3.383195

We already have substantial evidence.

meanBF = summary(afcDIACRITICS.modB)$coefficients["session4", "Estimate"]
seBF = summary(afcDIACRITICS.modB)$coefficients["session4", "Std. Error"]
h1mean = afc.s4.pred1.condPT.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.054526
## 
## $Likelihoodnull
## [1] 0.02467624
## 
## $BayesFactor
## [1] 42.73448

Session 4, using value as a maximum

meanBF = summary(afcDIACRITICS.modB)$coefficients["session4", "Estimate"]
seBF = summary(afcDIACRITICS.modB)$coefficients["session4", "Std. Error"]
h1mean = afc.s4.pred1.condPT.H1/2
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.9279899
## 
## $Likelihoodnull
## [1] 0.02467624
## 
## $BayesFactor
## [1] 37.60662

Prediction 1 Summary

Yes, in session 4 have strong evidence that performance is above chance; in Session 4 almost substantial evidence.

Question: Would they be above chance even without training?

Three-Interval Oddity Discrimination

Inspect Data

str(discrim)

## 'data.frame':    3672 obs. of  20 variables:
##  $ order         : int  25 26 27 28 29 30 31 32 33 34 ...
##  $ pt_N          : Factor w/ 51 levels "Participant 10",..: 10 10 10 10 10 10 10 10 10 10 ...
##  $ program       : Factor w/ 12 levels "S1_ver1_P.cond\\",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ same1         : Factor w/ 24 levels "chi_t2_fn1","chi_t3_fn1",..: 13 16 12 6 15 19 20 3 4 5 ...
##  $ same2         : Factor w/ 24 levels "chi_t2_fn2","chi_t3_fn2",..: 13 16 12 6 15 19 20 3 4 5 ...
##  $ diff          : Factor w/ 24 levels "chi_t2_fn3","chi_t3_fn3",..: 14 15 11 5 16 20 19 4 3 6 ...
##  $ info          : Factor w/ 24 levels "new_fff_t1t3",..: 12 11 24 8 9 18 21 2 7 4 ...
##  $ early_clicks  : Factor w/ 4 levels "also_earlyclick_s1",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ chose         : Factor w/ 3 levels " chosep1"," chosep2",..: 2 3 2 2 3 2 1 1 2 1 ...
##  $ RT            : num  1487 677 1254 5092 439 ...
##  $ response      : Factor w/ 2 levels "chosediff","chosesame": 1 1 1 1 1 2 2 1 2 1 ...
##  $ acc           : int  1 1 1 1 1 0 0 1 0 1 ...
##  $ session       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ oldnew        : Factor w/ 1 level "new": 1 1 1 1 1 1 1 1 1 1 ...
##  $ trial_type    : Factor w/ 3 levels "fff","ffm","fmf": 2 2 3 1 2 3 3 1 1 1 ...
##  $ tone_contrast2: Factor w/ 12 levels "t1t2","t1t3",..: 5 4 12 11 1 2 7 3 10 6 ...
##  $ tone_contrast : Factor w/ 6 levels "T12","T13","T14",..: 4 1 6 5 1 2 2 3 3 5 ...
##  $ version       : Factor w/ 2 levels "ver1","ver2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ cond          : Factor w/ 2 levels "P","P+T": 2 2 2 2 2 2 2 2 2 2 ...
##  $ item          : Factor w/ 24 levels "chit2","chit3",..: 14 15 11 5 16 20 19 4 3 6 ...

discrim$session = as.factor(discrim$session)

Get Means

dmeans = summarySEwithin(aggregate(acc ~ pt_N + cond + session, FUN = mean, data = discrim), measurevar = "acc", betweenvars = c("cond"), withinvars = c("session"))
kable(dmeans, digits = 2)

cond	session	N	acc	acc_norm	sd	se	ci
P	1	26	0.42	0.42	0.13	0.02	0.05
P	2	26	0.39	0.39	0.13	0.03	0.05
P	4	26	0.41	0.41	0.12	0.02	0.05
P+T	1	25	0.43	0.43	0.12	0.02	0.05
P+T	2	25	0.40	0.40	0.14	0.03	0.06
P+T	4	25	0.38	0.38	0.13	0.03	0.05

Get Plot

p2 = ggplot(dmeans, aes( x = cond, y = acc, fill = session))
p2 = p2 + geom_bar(stat = "identity", position = "dodge", colour = "blue", size = .2)
p2 = p2 + geom_errorbar(aes(ymin = acc-ci, ymax = acc+ci), width = .2, position = position_dodge(.9))
p2 = p2 + xlab("Condition") + ylab("Mean accuracy")
p2 = p2 + theme_bw()
p2 = p2 + coord_cartesian(ylim = c(.33, .5))
# = p2 +  geom_hline(aes(yintercept = .3))
p2

Axes from 33% (i.e. chance level)

Eyeballing: they are above chance even at pre-test and they don’t improve. Possible, they get worse, perhaps due to treating the test differently at pre and post test

Main LMEs

Here we look for an effect of condition (Pictures-Only vs. Pictures + Diacritics). We also consider whether performance changes from pretest (session 1) to post-test (either session 2 OR session 4).

We include trial type (fff, ffm, fmf), tone contrast (1-2, 1-3, 1-4, 2-3, 2-4, 3-4) as fixed effects. As before, these factors are expected to contribute to the model and are thus included but none were considered ofa priori interest.

discrim4 <- discrim

discrim4$session2 = as.factor(discrim4$session)
discrim4$session2 = revalue(discrim4$session2, c("1"="s1", "2"="s2","4"="s4"))

discrim4 = lizContrasts(discrim4, discrim4$session2, "s1")
discrim4 = lizContrasts(discrim4, discrim4$trial_type, "fff")
discrim4 = lizContrasts6(discrim4, discrim4$tone_contrast, "T12")
discrim4 = lizCenter(discrim4, c("cond"))


discrim4.mod1 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

# converges, now find slope structure for item

discrim4.mod2 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)*cond.ct|item),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

## Warning in commonArgs(par, fn, control, environment()): maxfun < 10 *
## length(par)^2 is not recommended.

# doesn't converge

discrim4.mod3 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)*cond.ct||item),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(discrim4.mod2,discrim4.mod3)

## Data: discrim4
## Models:
## discrim4.mod3: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod3:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod3:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) * cond.ct || item)
## discrim4.mod2: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod2:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod2:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod2:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) * cond.ct | item)
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## discrim4.mod3 25 4742.7 4897.9 -2346.3   4692.7                         
## discrim4.mod2 40 4765.6 5014.0 -2342.8   4685.6 7.0428     15     0.9564

discrim4.mod4 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)+cond.ct||item),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(discrim4.mod3,discrim4.mod4)

## Data: discrim4
## Models:
## discrim4.mod4: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod4:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod4:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod4:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) + cond.ct || item)
## discrim4.mod3: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod3:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod3:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) * cond.ct || item)
##               Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## discrim4.mod4 23 4738.7 4881.5 -2346.3   4692.7                        
## discrim4.mod3 25 4742.7 4897.9 -2346.3   4692.7     0      2          1

discrim4.mod5 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)||item),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

discrim4.mod6 = glmer (acc ~
          + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
          + (fff_VERSUS_fmf + fff_VERSUS_ffm)
          + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
          + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + (cond.ct||item),
          data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))


anova(discrim4.mod3,discrim4.mod5 )

## Data: discrim4
## Models:
## discrim4.mod5: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod5:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod5:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod5:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) || item)
## discrim4.mod3: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod3:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod3:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) * cond.ct || item)
##               Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## discrim4.mod5 22 4741.7 4878.3 -2348.8   4697.7                        
## discrim4.mod3 25 4742.7 4897.9 -2346.3   4692.7 4.987      3     0.1727

anova(discrim4.mod3, discrim4.mod6)

## Data: discrim4
## Models:
## discrim4.mod6: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod6:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod6:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod6:     pt_N) + (cond.ct || item)
## discrim4.mod3: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod3:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod3:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod3:     pt_N) + ((s1_VERSUS_s2 + s1_VERSUS_s4) * cond.ct || item)
##               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## discrim4.mod6 21 4735.2 4865.6 -2346.6   4693.2                         
## discrim4.mod3 25 4742.7 4897.9 -2346.3   4692.7 0.5494      4     0.9685

# can't remove condition but can remove session 
discrim4.mod = discrim4.mod6
round(summary(discrim4.mod)$coefficients, 3)

##                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)            -0.427      0.114  -3.735    0.000
## cond.ct                 0.010      0.114   0.084    0.933
## s1_VERSUS_s2           -0.157      0.088  -1.790    0.073
## s1_VERSUS_s4           -0.136      0.086  -1.577    0.115
## fff_VERSUS_fmf         -0.132      0.311  -0.425    0.671
## fff_VERSUS_ffm          0.570      0.310   1.837    0.066
## T12_VERSUS_T13         -0.395      0.411  -0.960    0.337
## T12_VERSUS_T14         -0.579      0.410  -1.412    0.158
## T12_VERSUS_T23         -0.005      0.409  -0.013    0.990
## T12_VERSUS_T24          0.043      0.409   0.104    0.917
## T12_VERSUS_T34         -0.239      0.379  -0.631    0.528
## cond.ct:s1_VERSUS_s2    0.019      0.175   0.109    0.913
## cond.ct:s1_VERSUS_s4   -0.162      0.173  -0.935    0.350

Note: intercept z value is incorrect since chance is not compared to 50% (log odds = 0) but to 1/3 (log odds = -0.6931472 Recalculation:

chance = logodds((1/3)*100)
Intercept_z = (summary(discrim4.mod)$coefficients["(Intercept)", "Estimate"] - chance)/summary(discrim4.mod)$coefficients["(Intercept)", "Std. Error"]
Intercept_p = 2*pnorm(-abs(Intercept_z))
round(Intercept_z,3)

## [1] 2.33

round(Intercept_p,6)

## [1] 0.019818

Above chance overall performance.

Prediction 1: There will be an improvement from pre to post test in both conditions

Frequentist

Fit slopes for contrasts between pre and each of the post tests for each condition

discrim4.modA = glmer (acc ~
            + (s1_VERSUS_s2 + s1_VERSUS_s4): cond
            + cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4)
            - (s1_VERSUS_s2 + s1_VERSUS_s4)
            - (s1_VERSUS_s2 + s1_VERSUS_s4): cond.ct
            + (fff_VERSUS_fmf + fff_VERSUS_ffm)
            + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
             + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
          + (cond.ct||item),
            data = discrim4, family = binomial, control = glmerControl(optimizer = "bobyqa"))

anova(discrim4.mod, discrim4.modA)

## Data: discrim4
## Models:
## discrim4.mod: acc ~ +cond.ct * (s1_VERSUS_s2 + s1_VERSUS_s4) + (fff_VERSUS_fmf + 
## discrim4.mod:     fff_VERSUS_ffm) + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + 
## discrim4.mod:     T12_VERSUS_T24 + T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | 
## discrim4.mod:     pt_N) + (cond.ct || item)
## discrim4.modA: acc ~ +(s1_VERSUS_s2 + s1_VERSUS_s4):cond + cond.ct * (s1_VERSUS_s2 + 
## discrim4.modA:     s1_VERSUS_s4) - (s1_VERSUS_s2 + s1_VERSUS_s4) - (s1_VERSUS_s2 + 
## discrim4.modA:     s1_VERSUS_s4):cond.ct + (fff_VERSUS_fmf + fff_VERSUS_ffm) + 
## discrim4.modA:     (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + 
## discrim4.modA:         T12_VERSUS_T34) + ((s1_VERSUS_s2 + s1_VERSUS_s4) | pt_N) + 
## discrim4.modA:     (cond.ct || item)
##               Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## discrim4.mod  21 4735.2 4865.6 -2346.6   4693.2                        
## discrim4.modA 21 4735.2 4865.6 -2346.6   4693.2     0      0          1

get_coeffs(discrim4.modA, c("s1_VERSUS_s2:condP", "s1_VERSUS_s2:condP+T", "s1_VERSUS_s4:condP", "s1_VERSUS_s4:condP+T"))

##                         Estimate Std. Error    z value   Pr(>|z|)
## s1_VERSUS_s2:condP   -0.16643936  0.1232375 -1.3505573 0.17683729
## s1_VERSUS_s2:condP+T -0.14734333  0.1248688 -1.1799854 0.23800601
## s1_VERSUS_s4:condP   -0.05706172  0.1209244 -0.4718792 0.63701298
## s1_VERSUS_s4:condP+T -0.21866628  0.1235134 -1.7703856 0.07666293

Bayes factors:

We use the HALF the value from the adult data (thus treating it as a maximum) a) Pictures-only condition, session 4

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP", "Estimate"]
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP", "Std. Error"]
h1mean = discrim.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.687546
## 
## $Likelihoodnull
## [1] 2.951508
## 
## $BayesFactor
## [1] 0.5717572

Evidence is ambiguous. Estimate of sample needed:

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 26, min=50, max=100))

##      x         y
## 1   50 0.4192691
## 2   51 0.4147885
## 3   52 0.4104108
## 4   53 0.4061325
## 5   54 0.4019501
## 6   55 0.3978602
## 7   56 0.3938598
## 8   57 0.3899458
## 9   58 0.3861153
## 10  59 0.3823655
## 11  60 0.3786940
## 12  61 0.3750981
## 13  62 0.3715755
## 14  63 0.3681238
## 15  64 0.3647409
## 16  65 0.3614247
## 17  66 0.3581731
## 18  67 0.3549841
## 19  68 0.3518560
## 20  69 0.3487870
## 21  70 0.3457753
## 22  71 0.3428192
## 23  72 0.3399172
## 24  73 0.3370679
## 25  74 0.3342696
## 26  75 0.3315210
## 27  76 0.3288207
## 28  77 0.3261675
## 29  78 0.3235600
## 30  79 0.3209971
## 31  80 0.3184776
## 32  81 0.3160003
## 33  82 0.3135641
## 34  83 0.3111681
## 35  84 0.3088111
## 36  85 0.3064923
## 37  86 0.3042106
## 38  87 0.3019651
## 39  88 0.2997550
## 40  89 0.2975794
## 41  90 0.2954375
## 42  91 0.2933284
## 43  92 0.2912514
## 44  93 0.2892058
## 45  94 0.2871908
## 46  95 0.2852057
## 47  96 0.2832497
## 48  97 0.2813224
## 49  98 0.2794230
## 50  99 0.2775509
## 51 100 0.2757054

Esimates would need sample of 75 to have evidence for the null.

2. Pictures-only condition, session 2

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP", "Estimate"]
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP", "Std. Error"]
h1mean = discrim.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.5219185
## 
## $Likelihoodnull
## [1] 1.300437
## 
## $BayesFactor
## [1] 0.4013408

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 26, min=26, max=50))

##     x         y
## 1  26 0.4013408
## 2  27 0.3921107
## 3  28 0.3832900
## 4  29 0.3748524
## 5  30 0.3667738
## 6  31 0.3590322
## 7  32 0.3516070
## 8  33 0.3444797
## 9  34 0.3376328
## 10 35 0.3310503
## 11 36 0.3247174
## 12 37 0.3186203
## 13 38 0.3127463
## 14 39 0.3070833
## 15 40 0.3016205
## 16 41 0.2963475
## 17 42 0.2912547
## 18 43 0.2863330
## 19 44 0.2815741
## 20 45 0.2769701
## 21 46 0.2725137
## 22 47 0.2681979
## 23 48 0.2640162
## 24 49 0.2599626
## 25 50 0.2560313

Estimate could have substantial evidence for null with 35 per condition (another 9 participants)

1. Pictures+ Diacritics only condition, session 4

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP+T", "Estimate"]
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP+T", "Std. Error"]
h1mean = discrim.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.234121
## 
## $Likelihoodnull
## [1] 0.6739056
## 
## $BayesFactor
## [1] 0.3474091

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min=26, max=50))

##     x         y
## 1  26 0.3382246
## 2  27 0.3295020
## 3  28 0.3212083
## 4  29 0.3133131
## 5  30 0.3057889
## 6  31 0.2986109
## 7  32 0.2917560
## 8  33 0.2852033
## 9  34 0.2789337
## 10 35 0.2729295
## 11 36 0.2671745
## 12 37 0.2616538
## 13 38 0.2563536
## 14 39 0.2512611
## 15 40 0.2463646
## 16 41 0.2416531
## 17 42 0.2371165
## 18 43 0.2327454
## 19 44 0.2285311
## 20 45 0.2244652
## 21 46 0.2205404
## 22 47 0.2167494
## 23 48 0.2130856
## 24 49 0.2095427
## 25 50 0.2061151

Esimates could have substantial evidence for H0 with 27 pariticipants (+1)

2. Pictures-only condition, session 2

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP+T", "Estimate"]
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP+T", "Std. Error"]
h1mean = discrim.pred1.H1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.6877976
## 
## $Likelihoodnull
## [1] 1.592604
## 
## $BayesFactor
## [1] 0.4318698

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min=26, max=50))

##     x         y
## 1  26 0.4221433
## 2  27 0.4128423
## 3  28 0.4039396
## 4  29 0.3954105
## 5  30 0.3872323
## 6  31 0.3793838
## 7  32 0.3718458
## 8  33 0.3646003
## 9  34 0.3576308
## 10 35 0.3509219
## 11 36 0.3444594
## 12 37 0.3382301
## 13 38 0.3322217
## 14 39 0.3264227
## 15 40 0.3208225
## 16 41 0.3154110
## 17 42 0.3101790
## 18 43 0.3051177
## 19 44 0.3002189
## 20 45 0.2954751
## 21 46 0.2908790
## 22 47 0.2864238
## 23 48 0.2821033
## 24 49 0.2779114
## 25 50 0.2738425

Esimates could have substantial evidence for H0 with 38 pariticipants (+ 13)

Summary

Pictures-only: evidence is ambiguous both for pre->s2 and pre->s4. For pre-> s2 could potentially get substantial evidence for H0 with 9 more participants. (in feasible numbers needed for pre=>s4)

Pictures + diacritics evidence is ambiguous both for pre->s2 and pre->s4. For pre-> s2 could potentially get substantial evidence for H0 with 13 more participants. for pre->s4 could potentially get substantial evidence with 1 more participant

Prediction 2: There will be greater improvement in one of the conditions (no direction specified)

Frequentist analysis

# Session 1 vs. 4
round(get_coeffs(discrim4.mod, c("cond.ct:s1_VERSUS_s2" , "cond.ct:s1_VERSUS_s4")),3)

##                      Estimate Std. Error z value Pr(>|z|)
## cond.ct:s1_VERSUS_s2    0.019      0.175   0.109    0.913
## cond.ct:s1_VERSUS_s4   -0.162      0.173  -0.935    0.350

Baysean Analysis

Given that we don’t find improvement in either condition, it doesn’t make sense to test this prediction.

Additional question A: Are children actually getting worse in the post-tests (all trial types)?

Begin by looking at this across trials

Frequentist

See previous prediction.

Bayes factors:

We don’t previous data to base H1 on here, but can work out plausible maximum in terms of what we would see if they returned to chance at post test

## get scores  pre-test
pre_test_scores = aggregate(acc ~  + cond + session, FUN = mean, data = aggregate(acc ~ pt_N + cond + session, FUN = mean, data = subset(discrim, session ==1)))
pre_test_scores$logodds = logodds(pre_test_scores$acc*100)
pre_test_scores

##   cond session       acc    logodds
## 1    P       1 0.4214744 -0.3167239
## 2  P+T       1 0.4333333 -0.2682640

## improvement in log odds, if returned to chance - half (since is maxium)  
discrim.newpred.P = (logodds(100/3) - pre_test_scores$logodds[1])/2
discrim.newpred.PT= logodds(100/3) - pre_test_scores$logodds[2]/2
discrim.newpred.P

## [1] -0.1882116

discrim.newpred.PT

## [1] -0.5590152

Note - the theory values here are negative, however this won’t work in BF function, therefore for each h1mean and meanBF *-1

1. Pictures-only condition, pre->session 4

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP", "Estimate"]*-1
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP", "Std. Error"]
h1mean = discrim.newpred.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.264873
## 
## $Likelihoodnull
## [1] 2.951508
## 
## $BayesFactor
## [1] 0.7673613

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 26, min=50, max=100))

##      x         y
## 1   50 0.7276975
## 2   51 0.7271063
## 3   52 0.7265728
## 4   53 0.7260949
## 5   54 0.7256708
## 6   55 0.7252986
## 7   56 0.7249768
## 8   57 0.7247036
## 9   58 0.7244775
## 10  59 0.7242971
## 11  60 0.7241611
## 12  61 0.7240681
## 13  62 0.7240168
## 14  63 0.7240061
## 15  64 0.7240349
## 16  65 0.7241019
## 17  66 0.7242064
## 18  67 0.7243471
## 19  68 0.7245233
## 20  69 0.7247339
## 21  70 0.7249782
## 22  71 0.7252553
## 23  72 0.7255644
## 24  73 0.7259048
## 25  74 0.7262757
## 26  75 0.7266765
## 27  76 0.7271065
## 28  77 0.7275650
## 29  78 0.7280515
## 30  79 0.7285653
## 31  80 0.7291059
## 32  81 0.7296727
## 33  82 0.7302652
## 34  83 0.7308830
## 35  84 0.7315255
## 36  85 0.7321923
## 37  86 0.7328829
## 38  87 0.7335969
## 39  88 0.7343339
## 40  89 0.7350935
## 41  90 0.7358753
## 42  91 0.7366789
## 43  92 0.7375040
## 44  93 0.7383502
## 45  94 0.7392173
## 46  95 0.7401048
## 47  96 0.7410125
## 48  97 0.7419400
## 49  98 0.7428872
## 50  99 0.7438536
## 51 100 0.7448391

Estimates would not have substantial evidence for H1 or H0 even with 100 participants.

2. Pictures-only condition, pre->session 2

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP", "Estimate"]*-1
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP", "Std. Error"]
h1mean = discrim.newpred.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.351407
## 
## $Likelihoodnull
## [1] 1.300437
## 
## $BayesFactor
## [1] 1.808166

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 26, min=27, max=100))

##      x         y
## 1   27  1.854557
## 2   28  1.902385
## 3   29  1.951693
## 4   30  2.002525
## 5   31  2.054928
## 6   32  2.108949
## 7   33  2.164636
## 8   34  2.222042
## 9   35  2.281218
## 10  36  2.342220
## 11  37  2.405103
## 12  38  2.469927
## 13  39  2.536750
## 14  40  2.605636
## 15  41  2.676649
## 16  42  2.749855
## 17  43  2.825323
## 18  44  2.903124
## 19  45  2.983332
## 20  46  3.066022
## 21  47  3.151272
## 22  48  3.239164
## 23  49  3.329781
## 24  50  3.423210
## 25  51  3.519540
## 26  52  3.618862
## 27  53  3.721272
## 28  54  3.826869
## 29  55  3.935754
## 30  56  4.048031
## 31  57  4.163810
## 32  58  4.283202
## 33  59  4.406322
## 34  60  4.533291
## 35  61  4.664231
## 36  62  4.799271
## 37  63  4.938541
## 38  64  5.082177
## 39  65  5.230321
## 40  66  5.383116
## 41  67  5.540713
## 42  68  5.703267
## 43  69  5.870938
## 44  70  6.043890
## 45  71  6.222294
## 46  72  6.406327
## 47  73  6.596170
## 48  74  6.792012
## 49  75  6.994046
## 50  76  7.202474
## 51  77  7.417502
## 52  78  7.639346
## 53  79  7.868225
## 54  80  8.104368
## 55  81  8.348012
## 56  82  8.599400
## 57  83  8.858785
## 58  84  9.126426
## 59  85  9.402592
## 60  86  9.687562
## 61  87  9.981622
## 62  88 10.285068
## 63  89 10.598207
## 64  90 10.921356
## 65  91 11.254841
## 66  92 11.599000
## 67  93 11.954183
## 68  94 12.320749
## 69  95 12.699071
## 70  96 13.089535
## 71  97 13.492538
## 72  98 13.908490
## 73  99 14.337817
## 74 100 14.780957

Estimates would have substantial evidence for H1 with 46 participants.

1. Pictures + diactirics condition, session 4

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP+T", "Estimate"]*-1
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s4:condP+T", "Std. Error"]
h1mean = discrim.newpred.PT *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.242655
## 
## $Likelihoodnull
## [1] 0.6739056
## 
## $BayesFactor
## [1] 1.84396

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min=27, max=100))

##      x          y
## 1   27   2.026674
## 2   28   2.126003
## 3   29   2.231039
## 4   30   2.342099
## 5   31   2.459520
## 6   32   2.583663
## 7   33   2.714907
## 8   34   2.853657
## 9   35   3.000343
## 10  36   3.155419
## 11  37   3.319368
## 12  38   3.492704
## 13  39   3.675967
## 14  40   3.869736
## 15  41   4.074619
## 16  42   4.291265
## 17  43   4.520359
## 18  44   4.762628
## 19  45   5.018845
## 20  46   5.289826
## 21  47   5.576438
## 22  48   5.879600
## 23  49   6.200284
## 24  50   6.539524
## 25  51   6.898414
## 26  52   7.278113
## 27  53   7.679851
## 28  54   8.104933
## 29  55   8.554741
## 30  56   9.030739
## 31  57   9.534482
## 32  58  10.067617
## 33  59  10.631891
## 34  60  11.229156
## 35  61  11.861375
## 36  62  12.530630
## 37  63  13.239130
## 38  64  13.989217
## 39  65  14.783372
## 40  66  15.624229
## 41  67  16.514581
## 42  68  17.457388
## 43  69  18.455792
## 44  70  19.513123
## 45  71  20.632914
## 46  72  21.818913
## 47  73  23.075093
## 48  74  24.405671
## 49  75  25.815118
## 50  76  27.308176
## 51  77  28.889878
## 52  78  30.565560
## 53  79  32.340884
## 54  80  34.221857
## 55  81  36.214851
## 56  82  38.326625
## 57  83  40.564351
## 58  84  42.935638
## 59  85  45.448559
## 60  86  48.111678
## 61  87  50.934083
## 62  88  53.925416
## 63  89  57.095908
## 64  90  60.456415
## 65  91  64.018456
## 66  92  67.794254
## 67  93  71.796780
## 68  94  76.039798
## 69  95  80.537914
## 70  96  85.306627
## 71  97  90.362389
## 72  98  95.722656
## 73  99 101.405953
## 74 100 107.431945

Estimates could have substantial evidence for H1 with 35 participants (another 10)

2. Pictures + diactirics condition, session 2

meanBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP+T", "Estimate"]*-1
seBF = summary(discrim4.modA)$coefficients["s1_VERSUS_s2:condP+T", "Std. Error"]
h1mean = discrim.newpred.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 2.393957
## 
## $Likelihoodnull
## [1] 1.592604
## 
## $BayesFactor
## [1] 1.503171

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min=27, max=100))

##      x        y
## 1   27 1.561405
## 2   28 1.591652
## 3   29 1.622674
## 4   30 1.654487
## 5   31 1.687107
## 6   32 1.720552
## 7   33 1.754840
## 8   34 1.789988
## 9   35 1.826018
## 10  36 1.862947
## 11  37 1.900798
## 12  38 1.939590
## 13  39 1.979346
## 14  40 2.020089
## 15  41 2.061841
## 16  42 2.104626
## 17  43 2.148469
## 18  44 2.193395
## 19  45 2.239429
## 20  46 2.286599
## 21  47 2.334932
## 22  48 2.384456
## 23  49 2.435200
## 24  50 2.487193
## 25  51 2.540467
## 26  52 2.595052
## 27  53 2.650980
## 28  54 2.708285
## 29  55 2.767001
## 30  56 2.827163
## 31  57 2.888805
## 32  58 2.951966
## 33  59 3.016682
## 34  60 3.082993
## 35  61 3.150938
## 36  62 3.220558
## 37  63 3.291895
## 38  64 3.364991
## 39  65 3.439891
## 40  66 3.516640
## 41  67 3.595285
## 42  68 3.675873
## 43  69 3.758452
## 44  70 3.843074
## 45  71 3.929790
## 46  72 4.018652
## 47  73 4.109715
## 48  74 4.203034
## 49  75 4.298666
## 50  76 4.396671
## 51  77 4.497107
## 52  78 4.600038
## 53  79 4.705525
## 54  80 4.813635
## 55  81 4.924433
## 56  82 5.037989
## 57  83 5.154371
## 58  84 5.273653
## 59  85 5.395908
## 60  86 5.521212
## 61  87 5.649642
## 62  88 5.781280
## 63  89 5.916205
## 64  90 6.054503
## 65  91 6.196260
## 66  92 6.341564
## 67  93 6.490507
## 68  94 6.643180
## 69  95 6.799680
## 70  96 6.960106
## 71  97 7.124557
## 72  98 7.293137
## 73  99 7.465952
## 74 100 7.643111

Could have substantial evidnece with 59 participants (34 more).

Additional question A: Summary

Are they getting worse? Note: currently all analyses suggest that the data are ambiguous.

However it seems likely that in fact, they are getting worse specficially on particular trial type where they are ABOVE chance at pre-test

Additional question B: Are children actually getting worse in the post-tests (easy test types)?

Could getting worse depend on TEST TYPE- given that some test items are much easier acoustically

Explore visually:

b = summarySEwithin(aggregate(acc ~ pt_N + cond + session + trial_type, FUN = mean, data = discrim), measurevar = "acc", betweenvars = c("cond"), withinvars = c("session", "trial_type"))
kable(b, digits = 2)

cond	session	trial_type	N	acc	acc_norm	sd	se	ci
P	1	fff	26	0.38	0.39	0.16	0.03	0.06
P	1	ffm	26	0.54	0.54	0.21	0.04	0.08
P	1	fmf	26	0.34	0.34	0.20	0.04	0.08
P	2	fff	26	0.38	0.38	0.14	0.03	0.06
P	2	ffm	26	0.45	0.45	0.20	0.04	0.08
P	2	fmf	26	0.33	0.33	0.18	0.03	0.07
P	4	fff	26	0.36	0.36	0.20	0.04	0.08
P	4	ffm	26	0.48	0.48	0.24	0.05	0.10
P	4	fmf	26	0.39	0.39	0.23	0.04	0.09
P+T	1	fff	25	0.38	0.38	0.16	0.03	0.06
P+T	1	ffm	25	0.57	0.57	0.20	0.04	0.08
P+T	1	fmf	25	0.34	0.34	0.20	0.04	0.08
P+T	2	fff	25	0.34	0.34	0.21	0.04	0.09
P+T	2	ffm	25	0.50	0.50	0.24	0.05	0.10
P+T	2	fmf	25	0.36	0.36	0.22	0.04	0.09
P+T	4	fff	25	0.34	0.34	0.17	0.03	0.07
P+T	4	ffm	25	0.46	0.45	0.25	0.05	0.10
P+T	4	fmf	25	0.36	0.35	0.19	0.04	0.08

p2 = ggplot(b, aes( x = cond, y = acc, fill = session))
p2 = p2 + facet_grid(.~ trial_type)
p2 = p2 + geom_bar(stat = "identity", position = "dodge", colour = "blue", size = .2)
p2 = p2 + geom_errorbar(aes(ymin = acc-ci, ymax = acc+ci), width = .2, position = position_dodge(.9))
p2 = p2 + xlab("Condition") + ylab("Mean accuracy")
p2 = p2 + theme_bw()
p2 = p2 + coord_cartesian(ylim = c(0, 1))
p2 = p2 + geom_hline(aes(yintercept = .3) )
p2

NOTE: the ffm trials and those that the participant can potentially do simply by assuming the male speaker is the odd man. It is possible they get worse specifically on these trials, because they stop doing this.

Frequentist

Look at s1->s2 and s1->s4 just for ffm trials

discrim4 <- discrim
discrim4$session2 = as.factor(discrim4$session)
discrim4$session2 = revalue(discrim4$session2, c("1"="s1", "2"="s2","4"="s4"))
discrim4 = lizContrasts(discrim4, discrim4$session2, "s1")

discrim4 = lizContrasts6(discrim4, discrim4$tone_contrast, "T12")
discrim4 <- lizCenter(discrim4, list("cond", "trial_type", "version"))


discrim4_ffm.mod = glmer (acc ~1
                            +
                (s1_VERSUS_s2 + s1_VERSUS_s4):cond       #+  (s1_VERSUS_s2 + s1_VERSUS_s4)
                + (T12_VERSUS_T13 + T12_VERSUS_T14 + T12_VERSUS_T23 + T12_VERSUS_T24 + T12_VERSUS_T34)
                + ((s1_VERSUS_s2 + s1_VERSUS_s4)|pt_N)
                 + (cond.ct||item),
       , data = subset(discrim4, trial_type=="ffm"), family = binomial, control = glmerControl(optimizer = "bobyqa"))

## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients

round( get_coeffs(discrim4_ffm.mod, c("s1_VERSUS_s2:condP", "s1_VERSUS_s2:condP+T", "s1_VERSUS_s4:condP", "s1_VERSUS_s4:condP+T")) , 3)

##                      Estimate Std. Error z value Pr(>|z|)
## s1_VERSUS_s2:condP     -0.425      0.214  -1.984    0.047
## s1_VERSUS_s2:condP+T   -0.311      0.217  -1.436    0.151
## s1_VERSUS_s4:condP     -0.292      0.225  -1.298    0.194
## s1_VERSUS_s4:condP+T   -0.543      0.229  -2.371    0.018

Bayes factors

use the same method as previous to obtain maxium value of H1

## get scores  pre-test
pre_test_scores = aggregate(acc ~  + cond + session, FUN = mean, data = aggregate(acc ~ pt_N + cond + session, FUN = mean, data = subset(discrim, session ==1 & trial_type == "ffm")))
pre_test_scores$logodds = logodds(pre_test_scores$acc*100)
pre_test_scores

##   cond session       acc   logodds
## 1    P       1 0.5432692 0.1735109
## 2  P+T       1 0.5750000 0.3022809

## improvement in log odds, if returned to change - half (since is maxium)  
discrim.newpredb.P = (logodds(100/3) - pre_test_scores$logodds[1])/2
discrim.newpredb.PT= logodds(100/3) - pre_test_scores$logodds[2]/2
discrim.newpredb.P

## [1] -0.4333291

discrim.newpredb.PT

## [1] -0.8442876

1. Pictures-only condition, session 4

meanBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s4:condP", "Estimate"]*-1
seBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s4:condP", "Std. Error"]
h1mean = discrim.newpredb.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.193849
## 
## $Likelihoodnull
## [1] 0.7627457
## 
## $BayesFactor
## [1] 1.5652

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 26, min=50, max=100))

##      x         y
## 1   50  2.770106
## 2   51  2.840557
## 3   52  2.913043
## 4   53  2.987623
## 5   54  3.064357
## 6   55  3.143309
## 7   56  3.224541
## 8   57  3.308121
## 9   58  3.394117
## 10  59  3.482600
## 11  60  3.573641
## 12  61  3.667317
## 13  62  3.763704
## 14  63  3.862881
## 15  64  3.964932
## 16  65  4.069940
## 17  66  4.177992
## 18  67  4.289179
## 19  68  4.403592
## 20  69  4.521328
## 21  70  4.642484
## 22  71  4.767162
## 23  72  4.895466
## 24  73  5.027504
## 25  74  5.163386
## 26  75  5.303227
## 27  76  5.447145
## 28  77  5.595260
## 29  78  5.747698
## 30  79  5.904587
## 31  80  6.066061
## 32  81  6.232256
## 33  82  6.403312
## 34  83  6.579376
## 35  84  6.760598
## 36  85  6.947130
## 37  86  7.139133
## 38  87  7.336770
## 39  88  7.540209
## 40  89  7.749626
## 41  90  7.965199
## 42  91  8.187114
## 43  92  8.415560
## 44  93  8.650733
## 45  94  8.892837
## 46  95  9.142080
## 47  96  9.398676
## 48  97  9.662847
## 49  98  9.934821
## 50  99 10.214833
## 51 100 10.503126

Estimates would substantial evidence with 53 particiants (27 more)

2. Pictures-only condition, session 2

meanBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s2:condP", "Estimate"]*-1
seBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s2:condP", "Std. Error"]
h1mean = discrim.newpredb.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.079745
## 
## $Likelihoodnull
## [1] 0.2603494
## 
## $BayesFactor
## [1] 4.147294

Have substantial evidence that they get worse.

1. Pictures + diactirics condition, session 4

meanBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s4:condP+T", "Estimate"]*-1
seBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s4:condP+T", "Std. Error"]
h1mean = discrim.newpredb.PT *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 0.744149
## 
## $Likelihoodnull
## [1] 0.1048665
## 
## $BayesFactor
## [1] 7.096156

Estimates could have substantial evidence for H1 with 35 participants (another 10)

2. Pictures + diactirics condition, session 2

meanBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s2:condP+T", "Estimate"]*-1
seBF = summary(discrim4_ffm.mod)$coefficients["s1_VERSUS_s2:condP+T", "Std. Error"]
h1mean = discrim.newpredb.P *-1
Bf(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1)

## $LikelihoodTheory
## [1] 1.20579
## 
## $Likelihoodnull
## [1] 0.6571788
## 
## $BayesFactor
## [1] 1.834798

data.frame(Bf_powercalc(seBF, meanBF, uniform = 0, meanoftheory = 0, sdtheory = h1mean, tail = 1, N = 25, min=25, max=100))

##      x         y
## 1   25  1.834798
## 2   26  1.890122
## 3   27  1.947518
## 4   28  2.007056
## 5   29  2.068812
## 6   30  2.132860
## 7   31  2.199283
## 8   32  2.268163
## 9   33  2.339588
## 10  34  2.413650
## 11  35  2.490442
## 12  36  2.570066
## 13  37  2.652622
## 14  38  2.738219
## 15  39  2.826967
## 16  40  2.918983
## 17  41  3.014386
## 18  42  3.113304
## 19  43  3.215864
## 20  44  3.322204
## 21  45  3.432463
## 22  46  3.546788
## 23  47  3.665331
## 24  48  3.788251
## 25  49  3.915710
## 26  50  4.047879
## 27  51  4.184937
## 28  52  4.327066
## 29  53  4.474459
## 30  54  4.627314
## 31  55  4.785837
## 32  56  4.950243
## 33  57  5.120756
## 34  58  5.297607
## 35  59  5.481036
## 36  60  5.671293
## 37  61  5.868639
## 38  62  6.073344
## 39  63  6.285686
## 40  64  6.505959
## 41  65  6.734465
## 42  66  6.971518
## 43  67  7.217445
## 44  68  7.472585
## 45  69  7.737292
## 46  70  8.011932
## 47  71  8.296886
## 48  72  8.592550
## 49  73  8.899334
## 50  74  9.217666
## 51  75  9.547991
## 52  76  9.890768
## 53  77 10.246479
## 54  78 10.615620
## 55  79 10.998709
## 56  80 11.396285
## 57  81 11.808906
## 58  82 12.237153
## 59  83 12.681630
## 60  84 13.142964
## 61  85 13.621806
## 62  86 14.118835
## 63  87 14.634753
## 64  88 15.170293
## 65  89 15.726215
## 66  90 16.303310
## 67  91 16.902397
## 68  92 17.524332
## 69  93 18.170000
## 70  94 18.840325
## 71  95 19.536263
## 72  96 20.258812
## 73  97 21.009005
## 74  98 21.787920
## 75  99 22.596673
## 76 100 23.436428

Estimates would have substantial evidence for H1 with 41 (16 more)

Additional question B Summary

Just looking at FFM trials (these start off well above chance - presumably children are responding based on the gender of odd man speaker) Pre-test -> session 4 post-test

Pictures-only: evidence is ambiguous (suggests that we need 54 participants, 28 more) Pictures + diacritics: substantial evidence they get worse

Pre-test -> session 2 post-test

Pictures-only: substantial evidence that they get worse Pictures + diacritics: evidence is ambiguous (suggests that we need 44 participants)

Word Repetition

Sound files need to be transcribed and rated by a native speaker before we can analyse this data.