"Patrones de herbivoría en Vassobia breviflora (Solanaceae): variación en el daño foliar y selección natural mediada por herbívoros
Valoy, Mariana; Mariano Ordano, Facundo Bernacki, Facundo X. Palacio; Juan Carlos López Acosta, Omar Varela

#Base de datos para el analisis figura en  Apendice.datos.xls
leaf<-read.table ("datos.txt", header=TRUE) #corresponde a base en Apendice.datos.xls
leaf$Herbivoria<-factor(leaf$Herbivoria)
leaf$n.planta<-factor(leaf$n.planta)
leaf$n.hoja<-factor(leaf$n.hoja)

sin.herb<-rownames(leaf[leaf$Herbivoria=="0",])
leaf1<-leaf[-as.numeric(sin.herb),]


#RESUMENES DE DATOS
#X,SD,CV para poblacion
#variable Proporcion de area folir removida (pafr)
tsitio<-with(leaf,tapply(pafr,list(sitio),mean))
tsitiosd<-with(leaf,tapply(pafr,list(sitio),sd))
tsitiomax<-with(leaf,tapply(pafr,list(sitio),max))
cv.cadi<-0.1180/0.08513
cv.lahorqueta<-0.1656/0.1445
cv.Monti<-0.1408/0.1991
cv.san.ignacio<-0.1125/0.1192

t1<-with(leaf,tapply(pafr,list(Tipo.herb,sitio),mean))
t2<-with(leaf,tapply(pafr,list(Tipo.herb,sitio,edad.hoja),mean))
t3<-with(leaf,tapply(pafr,list(sitio,edad.hoja),mean))
#desvios
t1sd<-with(leaf,tapply(pafr,list(Tipo.herb,sitio),sd))
t2sd<-with(leaf,tapply(pafr,list(Tipo.herb,sitio,edad.hoja),sd))
t3sd<-with(leaf,tapply(pafr,list(sitio,edad.hoja),sd))

#variable area foliar (af)
tsitio<-with(leaf,tapply(afl,list(sitio),mean))
tsitiosd<-with(leaf,tapply(af,list(sitio),sd))
tsitiomax<-with(leaf,tapply(afr,list(sitio),max))

tsitio<-with(leaf,tapply(af,list(sitio),mean))
tsitiosd<-with(leaf,tapply(af,list(sitio),sd))

tarea<-with(leaf,tapply(af,list(Tipo.herb,sitio),mean))
tarea2<-with(leaf,tapply(af,list(Tipo.herb,sitio,edad.hoja),mean))
#desvio
tareasd<-with(leaf,tapply(af,list(Tipo.herb,sitio),sd))
tarea2sd<-with(leaf,tapply(af,list(Tipo.herb,sitio,edad.hoja),sd))

#variable laf (relacion longitud/ancho)
tsitio<-with(leaf,tapply(laf,list(sitio),mean))
tsitiosd<-with(leaf,tapply(laf,list(sitio),sd))
tlonga<-with(leaf,tapply(laf,list(Tipo.herb,sitio),mean))
tlonga2<-with(leaf,tapply(laf,list(Tipo.herb,sitio,edad.hoja),mean))
#desvio
tlongasd<-with(leaf,tapply(laf,list(Tipo.herb,sitio),sd))
tlonga2sd<-with(leaf,tapply(laf,list(Tipo.herb,sitio,edad.hoja),sd))



#cambio nombre de sitio para graficar
levels(leaf2$sitio)[levels(leaf2$sitio)=="Cadillal"]<-"A.Chiquitas"
levels(leaf2$sitio)[levels(leaf2$sitio)=="lahorqueta"]<-"La Horqueta"
levels(leaf2$sitio)[levels(leaf2$sitio)=="Monti"]<-"V.P.Monti"
levels(leaf2$sitio)[levels(leaf2$sitio)=="san.ignacio"]<-"S.Ignacio"
library(lattice)

#Reduccion de los niveles del factor Tipo.herb
levels(Tipo.herb)
summary(Tipo.herb)

leaf2<-leaf1
levels(leaf2$Tipo.herb)
levels(leaf2$Tipo.herb)[levels(leaf2$Tipo.herb)=="cortadorypunt"]<-"mix"
levels(leaf2$Tipo.herb)[levels(leaf2$Tipo.herb)=="minadorypunt"]<-"mix"
levels(leaf2$Tipo.herb)[levels(leaf2$Tipo.herb)=="daño?"]<-"mix"
levels(leaf2$Tipo.herb)

#Elimino san ignacio por que no tiene ninguna de las categorias de daño que quedan
levels(leaf2$Tipo.herb)[levels(leaf2$Tipo.herb)=="minador"]<-"mix"
levels(leaf2$Tipo.herb)
leaf3<-subset(leaf2,!sitio=="san.ignacio")
m<-with(leaf3,tapply(prop.afr,list(Tipo.herb,sitio,edad.hoja),length))
leaf4<-subset(leaf3,!Tipo.herb=="mix")

#Modelos
#Regresion Beta
library(betareg)
betamod<-betareg(pafr~edad.hoja + sitio*Tipo.herb, data=leaf2)
summary(betamod)

#Modelos Lineales Generales Mixtos
library(nlme)
lme1<-lme(af~afr*edad.hoja*Tipo.herb, data=leaf4, random= ~1|sitio/n.planta)
summary(lme1)

lme2<-lme(laf~pafr*Tipo.herb*edad.hoja, data=leaf4, random= ~1|sitio/n.planta)
summary(lme2)

#Graficos

boxplot(prop.afr~ factor(sitio),
        varwidth = TRUE, xlab = "Sitio",
        main = "Boxplot del Area foliar removida sobre\
        Sitio", ylab = "Area foliar removida", data = leaf)

boxplot(af~ factor(sitio),
        varwidth = TRUE, xlab = "Sitio",
        main = "Boxplot del Area foliar removida sobre\
        Sitio", ylab = "Area foliar removida", data = leaf)

boxplot(afr~ factor(sitio),
        varwidth = TRUE, xlab = "Sitio",
        main = "Boxplot del Area foliar removida sobre\
        Sitio", ylab = "Area foliar removida", data = leaf)

boxplot(laf~ factor(sitio),
        varwidth = TRUE, xlab = "Sitio",
        main = "Boxplot del Area foliar removida sobre\
        Sitio", ylab = "Area foliar removida", data = leaf)



library(lattice)
par(cex.lab=3,cex.axis=6,cex.sub=3)
main<-(par(cex.lab=3,cex.axis=6,cex.sub=3))

lattice.options(default.theme = standard.theme(color = FALSE))
sb <- trellis.par.get("strip.background") 
sb[["col"]][1] <- "lightgray"
trellis.par.set("strip.background", sb)

bwplot(areatotal~Tipo.herb|sitio,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Pafr",col="black", par.strip.text = list(cex=2),layout=c(1,4),                                                                                       
        auto.key=T)
bwplot(af~Tipo.herb|edad.hoja,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Pafr",col="black",par.strip.text = list(cex=2),layout=c(1,2), auto.key=T)

bwplot(pafr~Tipo.herb|sitio,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Pafr",col="black", par.strip.text = list(cex=2),layout=c(1,4),auto.key=T)

bwplot(pafr~Tipo.herb|edad.hoja,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Pafr",col="black",par.strip.text = list(cex=2),layout=c(1,2), auto.key=T)

bwplot(laf~Tipo.herb|sitio,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Longitud/ancho foliar",col="black",par.strip.text = list(cex=2),layout=c(1,4), auto.key=T)

bwplot(laf~Tipo.herb|edad.hoja,data=leaf2, xlab="Tipo de herbivoria",
       ylab="Longitud/ancho foliar",col="black",par.strip.text = list(cex=2),layout=c(1,2), auto.key=T)