## R code to replicate analyses in "Can squirrel monkeys learn an ABnA
## grammar? A re-evaluation of Ravignani et al. (2013)

## load the data.table library and read data:
library(data.table)
d <- fread("Ravignani2013Reanalysis.csv")
## introduce a Grammatical variable, first set No for all Sequences,
## then set to Yes for those that match the grammar:
d[, Grammatical := "No" ]
d[ grepl( "AB+A", Sequence ), Grammatical := "Yes" ]

## the manuscript was prepared in Org Mode for Emacs. this line loads
## the ascii library and sets it to use the appropriate output format:
library(ascii)
options(asciiType="org")

# results of Tests 1 and 2, all data included: 
d.bars <- d[ , mean(Response>0), by=.(Test,Grammatical) ]
setkey( d.bars, Grammatical )
setnames( d.bars, "V1", "Proportion of trials with $>0$ responses" )
ascii( d.bars[ order(Test) ], format=c("d","s","f"), include.rownames=FALSE )

# results of Test 2 after exclusion of AB and BA:
d.bars <- d[ Test==2 & ! Sequence %in% c("AB","BA"), mean(Response>0), by=.(Test,Grammatical) ]
setkey( d.bars, Grammatical )
setnames( d.bars, "V1", "Proportion of trials with $>0$ responses" )
ascii( d.bars, format=c("d","s","f"), include.rownames=FALSE )

## a function to automate the t-tests somewhat. arguments:
##
## test: which test (1 or 2)
##
## dich: whether to dichotomize or not (TRUE or FALSE)
##
## excl: character vector of sequences to exclude
##
## what: what part of the t.test results to return, such as 'estimate'
## for the t value and 'pvalue' for the p value. (note 'pvalue' rather
## than the standard R p.value because this function is called from a
## lisp interface that does not like dots in variable names.)
ravignani.t.test <- function( test, dich, excl, what ) {
    x <- d[ Test==test & ! Sequence %in% excl ]
    if( dich ) {
	x <- x[,
	    mean( Response[Grammatical=="No"]>0 ) -
	    mean( Response[Grammatical=="Yes"]>0 ),
	    by=Subject
	]
    } else {
	x <- x[,
	    mean( Response[Grammatical=="No"] ) -
	    mean( Response[Grammatical=="Yes"] ),
	    by=Subject
	]
    }
    t.out <- t.test( x$V1 )
    t.out$pvalue <- t.out$p.value
    if( what == "statistic" ) {
	d <- 3
    } else {
	d <- 2
    }
    result <- format( t.out[[what]], digits=d )
    if( what == "pvalue" ) {
	if( t.out[[what]] < 0.05 ) {
	    result <- paste0( result, "$^\\star$" )
	}
    }
    result
}

x <- ""
## test 1, dichotomization, no sequences excluded:
ravignani.t.test( 1, TRUE, c(), x )

x <- ""
## test 2, dichotomization, AB and BA excluded:
ravignani.t.test( 2, TRUE, c("AB","BA"), x )

x <- ""
## test 1, no dichotomization, no sequences excluded
ravignani.t.test( 1, FALSE, c(), x )

x <- ""
## test 2, no dichotomization, AB and BA excluded:
ravignani.t.test( 2, FALSE, c("AB","BA"), x )

x <- ""
## test 1, dichotomization, n=1,2,3 sequences excluded:
ravignani.t.test( 1, TRUE, c("ABA","ABBA","ABBBA"), x )

x <- ""
## test 1, dichotomization, n=4,5 sequences excluded:
ravignani.t.test( 1, TRUE, c("ABBBBA","ABBBBBA"), x )

x <- ""
## test 1, no dichotomization, n=1,2,3 sequences excluded:
ravignani.t.test( 1, FALSE, c("ABA","ABBA","ABBBA"), x )

x <- ""
## test 1, no dichotomization, n=4,5 sequences excluded:
ravignani.t.test( 1, FALSE, c("ABBBBA","ABBBBBA"), x )

x <- ""
## test 2, dichotomization, no sequences excluded:
ravignani.t.test( 2, TRUE, c(), x )

x <- ""
## test 2, no dichotomization, no sequences excluded:
ravignani.t.test( 2, FALSE, c(), x )

## ANOVA in the original paper, with a dichotmized response and AB, BA excluded:
library(car)
t <- Anova( lm( I(Response>0) ~ Grammatical * Test, d, subset=!Sequence %in% c("AB","BA")) )
names(t) <- c("Sum of Squares","d.f.","$F$","$p$")
ascii(t, format=c("fg","d","f","fg"), digits=c(2,0,2,2))

## graph of the distribution of responses in both Test 1 and 2:
pdf( "Responses.pdf", height=4, width=10 )
par( mfrow=c(1,2), mar=c(4,4,4,1) )
d.table <- d[ !Sequence %in% c("AB","BA"), .N, by=.(Test,Grammatical,Response) ]
d.table <- d[ , .N, by=.(Test,Grammatical,Response) ]
setkey( d.table, Response )
plot( c(0,4), c(0,25), pch=NA, yaxs="i", xlab="Number of Responses", ylab="Frequency", main="Test 1" )
d.table[ Test==1 & Grammatical=="No", lines( Response, N, type="o", pch=16, xpd=TRUE ) ]
d.table[ Test==1 & Grammatical=="Yes", lines( Response, N, lty=2, type="o", pch=17 ) ]
plot( c(0,4), c(0,25), pch=NA, yaxs="i", xlab="Number of Responses", ylab="", main="Test 2" )
d.table[ Test==2 & Grammatical=="No", lines( Response, N, type="o", pch=16, xpd=TRUE ) ]
d.table[ Test==2 & Grammatical=="Yes", lines( Response, N, lty=2, type="o", pch=17 ) ]
legend( "topright", legend=c("Grammatical","Non-grammatical"), pch=c(17,16), lty=c(2,1), bty="n" )
dev.off()

## ANOVA without dichotomization, but still excluding AB, BA:
t <- Anova( lm( Response ~ Grammatical * Test, d, subset=!Sequence %in% c("AB","BA")) )
names(t) <- c("Sum of Squares","d.f.","$F$","$p$")
ascii(t, format=c("fg","d","f","fg"), digits=c(2,0,2,2))