## 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", "Fraction 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", "Mean" )
ascii( d.bars, format=c("d","s","f"), include.rownames=FALSE )

## results of t-test for Test 1:
d.mean <- d[Test==1,
            mean( Response[ Grammatical=="No"  ] > 0 ) -
            mean( Response[ Grammatical=="Yes" ] > 0),
            by=Subject ]
t.test( d.mean$V1 )

## results of t-test for Test 2 after excluding AB and BA:
d.mean <- d[Test==2 & ! Sequence %in% c("AB","BA"),
            mean( Response[ Grammatical=="No"  ] > 0 ) -
            mean( Response[ Grammatical=="Yes" ] > 0),
            by=Subject ]
t.test( d.mean$V1 )

## 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()

## t-test for Test 1, without dichotomization:
d.mean2 <- d[ Test==1,
             mean(Response[Grammatical=="Yes"]) -
             mean(Response[Grammatical=="No"]),
             by=Subject ]
t.test( d.mean2$V1 )

## t-test for Test 1, without dichotomization, but excluding AB and BA:
d.mean2 <- d[ Test==2 & !Sequence %in% c("AB","BA"),
             mean(Response[Grammatical=="Yes"]) -
             mean(Response[Grammatical=="No"]),
             by=Subject ]
t.test( d.mean2$V1 )

## 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))

## t-test for Test 2, including AB, BA, with dichotomization:
d.mean2 <- d[ Test==2,
             mean(Response[Grammatical=="Yes"]>0) -
             mean(Response[Grammatical=="No"]>0),
             by=Subject ]
t.test( d.mean2$V1 )

## t-test for Test 2, including AB, BA, and without dichotomization:
d.mean2 <- d[ Test==2,
             mean(Response[Grammatical=="Yes"]) -
             mean(Response[Grammatical=="No"]),
             by=Subject ]
t.test( d.mean2$V1 )

## table of fraction of responses and mean of responses to grammatical
## and non-grammatical patterns in Test 1 and 2:
i <- d[, ! (Sequence %in% c("AB","BA") & Test==2) ]
d.sum <- d[ i, .(mean(Response>0), mean(Response)), by=.(Test,Grammatical) ]
setkey( d.sum, Grammatical )
setnames( d.sum, c("V1","V2"), c("Fraction of Trials", "Mean Responses") )
ascii( d.sum[order(Test)], include.rownames=FALSE, format=c("d","s","f","f") )

## logistic regression with dichotomized responses, excluding AB, BA:
mod <- glm( I(Response>0) ~ Grammatical * Test,
                 d, subset=!Sequence %in% c("AB","BA"), family="binomial" )
ascii( Anova( mod ), format=c("f","d","fg"), digits=c(2,0,2) )

## Poisson regression (non-dichotomized), excluding AB, BA:
mod <- glm( Response ~ Grammatical * Test,
                 d, subset=!Sequence %in% c("AB","BA"), family="poisson" )
ascii( Anova( mod ), format=c("f","d","fg"), digits=c(2,0,2) )