#############################################################################################
### Analisis de Abundancia de aves en parches de Vegetacion de Victoria, Australia        ###
### Loyn 1987   Base de datos usada por Quinn & Keough (2002) y  Zuur cols (2009)         ###
#############################################################################################

R

library(bbmle)
library(car)


#Utilizando la funcion read.table

abav <- read.table("abun_aves.csv", header=T, sep=",")
str(abav)

#Forrajeo es una variable predictora categorica
abav$forrajeo <- factor(abav$forrajeo)

#Inspeccion visual
par(mfrow=c(2,3))
dotchart(abav$area, xlab="area")
dotchart(abav$aislamiento, xlab="aislamiento")
dotchart(abav$distpar, xlab="Distancia al parche mas cercano")
dotchart(abav$distparmay, xlab= "distancia a un parche mayor")
dotchart(abav$altitud, xlab="altitud")

boxplot(abundancia~forrajeo,data=abav)


#Si continuaramos sin transformar las variables solo con fines demonstrativos

m1 <- lm(abundancia~area + aislamiento + distpar + distparmay + forrajeo + altitud, data=abav)
par(mfrow=c(2,2))
plot(m1)

summary(m1)


#Transformaciones necesarias en las variables predictoras

abav$larea <- log(abav$area)
abav$ldistpar <- log(abav$distpar)
abav$ldistparmay <- log(abav$distparmay)

par(mfrow=c(1,3))
dotchart(abav$larea, xlab="area")
dotchart(abav$ldistpar, xlab="Distancia al parche mas cercano")
dotchart(abav$ldistparmay, xlab= "distancia a un parche mayor")


m2<-lm(abundancia~ larea + aislamiento + ldistpar + ldistparmay + forrajeo + altitud, data=abav)
par(mfrow=c(2,2))
plot(m2) #los residuos no presentan dist normal

summary(m2)
anova(m1);anova(m2) #los valores de probabilidad cambiaron

#Evaluando el grado de autocorrelacion entre las variables predictoras
scatterplotMatrix(matriz)


matriz <- cbind(abav$larea,abav$aislamiento,abav$ldistpar,abav$ldistparmay,abav$altitud)
cor(matriz)


m2b<-lm(abundancia~ larea + aislamiento + ldistpar + ldistparmay + altitud, data=abav) #Sin forrajeo porque es un factor
#Evaluando los factores de inflacion de la varianza, 
vif(m2b) #criterio para eliminar VIF>3, Zuur cols.


#Seleccion de modelos, "stepwise selection" usando el valor de AIC como criterio
step(m2)

m3 <- lm(abundancia ~larea+forrajeo, data=abav)
vif(m3)
summary(m3)
anova(m3)

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

#Distancia de Cook, punto de quiebra propuesto por John Fox

#Distancias de Cook, cutoff mas estricto: 4/(n-p), segun Fox (2002) chp 6.1
coff<-4/(dim(abav)[1]-length(coef(m3)))
coff


#Graficando las distancias de Cook
plot(cookd(m3))
abline(h=coff)

dcook <-cookd(m3)
which(dcook>coff)

m3b<- lm(abundancia ~larea+forrajeo, data=abav, subset=c(-9,-18,-31,-56))
anova(m3);anova(m3b)
#Cambio en los valores de F

##Seleccion de modelos
m4 <- update(m2,.~.-altitud)
m5 <-update(m4,.~.-forrajeo)
m6 <- update(m5,.~.-ldistparmay)
m7 <- update(m6,.~.-aislamiento)
m8 <- update(m7,.~.-larea)

#Seleccion de modelos ordenados, y grado de soporte relativo
AICtab(m2,m3,m4,m5,m6,m7,m8, base=TRUE, weights=TRUE)

#Prueba de razon de verosimilitud (o de F), entre el modelo global y el m4 que no cuenta con el efecto del forrajeo
anova(m3,m5)

#Ejercicio, analisis de la base de datos divusa que analiza la tasa de divorcios en los estados unidos ubicada en el paquetefaraway, 
#detectar la presencia de posibles variables correlacinadas, realizar la seleccion y validacion del modelo.