anova+Tuk2016.R

setwd("C:/Users/AlvaroA/Google Drive/ColaB/BehEco/Claudia")
datosORI <- read.table("DatosDegusTerreno.csv", sep=";", header=TRUE)
datos <- read.table("DatosDegusTerrenoOUT.csv", sep=";", header=TRUE)
datosALT <- read.table("DatosDegusTerrenoALT.csv", sep=";", header=TRUE)

# pairs(~Frec_Peak+Frec_Baja_prin+Frec_alta_prin+Frec_Baja_fin+Frec_alta_fin,
#       data=datosOUT)
# pairs(~Frec_Peak+Frec_Baja_prin+Frec_alta_prin+Frec_Baja_fin+Frec_alta_fin,
#       data=datosALT)


repeats <- 9999

# myFormula <- formula(Dur_silaba ~ stim)
# myFormulR <- formula(sample(Dur_silaba) ~ stim)
# myFormula <- formula(Frec_Peak ~ stim)
# myFormulR <- formula(sample(Frec_Peak) ~ stim)
# myFormula <- formula(ancho_banda ~ stim)
# myFormulR <- formula(sample(ancho_banda) ~ stim)

# myFormula <- formula(acenso ~ stim)
# myFormulR <- formula(sample(acenso) ~ stim)
myFormula <- formula(log(duracion_total) ~ stim)
myFormulR <- formula(sample(log(duracion_total)) ~ stim)

# myFormula <- formula(log(sobretonos) ~ stim)
# myFormulR <- formula(sample(log(sobretonos)) ~ stim)

model1 <- lm(myFormula, data = datos)
anova1 <- anova(model1)

Fa1 <- anova1$F[1]
A <- c(Fa1)
cntA <- 0

Tuk1 <- TukeyHSD(aov(model1))
nTa <- dimnames(Tuk1$stim)[[1]]
dTa <- Tuk1$stim[1:3]
pTa <- Tuk1$stim[10:12]
cTa <- c(1, 1 , 1)
kTa <- c(1, 1 , 1)

start <- Sys.time()
for (i in 1:repeats){
    modelR <- lm(myFormulR, data = datos)
    anovaR <- anova(modelR)
    FaR <- anovaR$F[1]; A[i+1] <- FaR
    if (FaR >= Fa1) {cntA <- cntA+1}

    TukR <- TukeyHSD(aov(modelR))
    dTaR <- TukR$stim[1:3]
    kTa <- kTa + as.integer(abs(dTaR) >= abs(dTa))
    pTaR <- TukR$stim[10:12]
    cTa <- cTa + as.integer(pTaR <= pTa)
    
    }
end <- Sys.time()
dif <- as.double(difftime(end, start, units = "secs"))
h <- dif%/%3600; m <- dif%%3600%/%60; s <- dif%%3600%%60
elapsed <- paste(h,"h ",m,"min ",format(s, digits=5),"sec",sep="")

sink("anova+Tuk2016.txt", append = T)
print(anova1)

## Analysis of Variance Table
## 
## Response: log(duracion_total)
##           Df Sum Sq Mean Sq F value   Pr(>F)   
## stim       2 106.66  53.328  7.8442 0.001532 **
## Residuals 35 237.95   6.798                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat(
    "-----",
    paste("PERMUTACIONES =", repeats),
    paste("p-stim        =", (cntA+1)/(repeats+1)),
    paste("tiempo usado  =", elapsed),
    "--------------------------------------------------------------",
    "",
    sep="\n")

## -----
## PERMUTACIONES = 9999
## p-stim        = 0.0023
## tiempo usado  = 0h 1min 53.481sec
## --------------------------------------------------------------

print(data.frame(contrast = nTa, diff = dTa, p.Tukey = pTa,
  Perm.p = cTa/(repeats+1), Perm.diff = kTa/(repeats+1)), digits=log10(repeats+1))

##       contrast    diff  p.Tukey Perm.p Perm.diff
## 1   fox-flying  4.3968 0.003179 0.0011    0.0016
## 2 human-flying  0.1264 0.990516 0.8946    0.9113
## 3    human-fox -4.2704 0.002052 0.0013    0.0013

cat("",sep="\n", fill=T)

cat(as.character(end))

## 2016-09-08 12:35:13

cat("",sep="\n", fill=T)

cat(rep("=",79),sep="", fill=T)

## ===========================================================================
## ====

sink()


par(mfrow = c(1, 2))

library(yarrr)

## Loading required package: jpeg

## Loading required package: BayesFactor

## Loading required package: coda

## Loading required package: Matrix

## ************
## Welcome to BayesFactor 0.9.12-2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

## ************
## Welcome to YaRrr! The Pirate's Guide to R

## Run yarrr.guide() to see the main package guide

## Email me at yarrr.book@gmail.com with questions, comments, or movie recommendations
## ************

if (toString(myFormula[2]) == "Dur_silaba") {
  ylabel = "Syllable duration"
} else if (toString(myFormula[2]) == "Frec_Peak") {
  ylabel = "Peak frequency"
} else if (toString(myFormula[2]) == "ancho_banda") {
  ylabel = "Bandwidth"
} else if (toString(myFormula[2]) == "acenso") {
  ylabel = "Rise time"
} else if (toString(myFormula[2]) == "log(duracion_total)") {
  ylabel = "log(Total duration)"
} else
  ylabel = "something..."

pirateplot(formula = myFormula,
           data=datos,
           line.o=.7,
           inf.o=.3, 
           ylab=ylabel, 
           xlab="Stimulus")

hist(A, prob = F, 
     main = paste("F = ", as.name(Fa1), ", p = ", as.name((cntA+1)/(repeats+1))), 
     xlab = "Permutation Fs")
abline(v=Fa1)
# lines(density(A))
rug(A)

file.show("anova+Tuk2016.txt", title = "Cecchi's results... by the mighty riverula")