R
library(ISwR)
library(faraway)
library(car)
library(nlme)


#Crecimiento de plantulas en condiciones de longnosidad contrastantes, una semana, medicion de la longitud
#fijar la generacion de valores pseudoaleatorios para fines de contar con los mismos estadisticos y valores de p. en clase
set.seed(123)

#observaciones por categoria, ahora con desviaciones estandar diferentes
n<-40
l<-rnorm(n,mean=8.5,sd=1.6)
s<-rnorm(n,mean=8.0,sd=0.7)

long<-c(l,s)
cat <-c(rep("l",n),rep("s",n))
plantula <- data.frame(long,cat)

tapply(plantula$long,plantula$cat,mean)
tapply(plantula$long,plantula$cat,sd)

#Inspeccion visual
boxplot(long~cat, data=plantula,ylab="Longitud de la plantula (mm)")

#Considerando la varianza es homogenea para cada grupo
tig <-t.test(long~cat, data=plantula, var.equal=TRUE)
tig

#Varianzas diferentes
tdesig <-t.test(long~cat, data=plantula)
tdesig

#Analisis de varianza, niveles de folato en globulos rojos, diferentes tipos de vent
data(red.cell.folate)
help(red.cell.folate)
str(red.cell.folate)

#Inspeccion visual
rdf <-red.cell.folate

plot(folate~ventilation,rdf)

#Ajuste del modelo
m1 <- lm(folate~ventilation, data=rdf)


#Validacion
par(mfrow=c(2,2))
plot(m1)


boxplot(resid(m1)~ rdf$ventilation)

anova(m1)

#Usando la funcion gls

m1gls <- gls(folate~ventilation, data=rdf)
summary(m1);summary(m1gls)


#Especificando varianzas diferentes para cada grupo
vfi <- varIdent(form= ~1 |ventilation)
m2gls <- gls(folate~ventilation,weights=vfi,data=rdf)

#Comparando los dos modelos

anova(m1gls,m2gls)

#El modelo con mejor soporte es aquel que considera las varianzas homogeneas


#Concentracion de ozono en diferentes jardines, tomado de Crawley 2007 R Book
jrd <- read.table("jardines.txt", header=T, sep= "\t")
str(jrd)
dim(jrd)

ozono <- c(jrd$gardenA, jrd$gardenB, jrd$gardenC)
jarcat <-rep(c("A","B","C"), each=10)
jarcat

gard <- data.frame(ozono,jarcat)
head(gard)

#Inspeccion visual
boxplot(ozono~jarcat,data=gard)

#cuando tenemos pocos puntos, alternativa a la grafica de cajas
stripchart(ozono~jarcat,data=gard,method="jitter",vertical=TRUE, pch=19)

#Pruebas para determinar si la varianza es clasica, pero sensibles a la normalidad
#Zuur cols (2009) desaconsejan su uso, Prueba de Barlett

bartlett.test(ozono~jarcat, data=gard)

#Ajustando el modelo de regresion lineal y su correspondiente gls

m3 <- lm(ozono~jarcat, data=gard)

ncvTest(m3) #Prueba de Cook y Weisberg para probar la homogeneidad de varianza

par(mfrow=c(2,2))
plot(m3)

m3gls <- gls(ozono~jarcat, data=gard)
anova(m3gls)

#Ajustando un modelo con varianzas diferents para cada grupo
vfi2 <- varIdent(form = ~1|jarcat)
m3glsvi <- update(m3gls, weights=vfi2)

anova(m3gls,m3glsvi)
AICtab(m3gls,m3glsvi, base=TRUE,weights=TRUE)

summary(m3gls);summary(m3glsvi)
anova(m3gls)