#===============================================================
###########  CONECTAR CON EL DIRECTORIO DE TRABAJO  ############
#===============================================================
	
	setwd("C:/Deforestacion/Choco") 
	dir()

#==========================================================================
###############         CARGAR PAQUETES NECESARIOS          ###############
#==========================================================================
	
	require(maptools)  # Para cargar los datos readpoly
	require(GeoXp)     # Para ejecutar el Histomap
	require(spdep)	 # para matrices de pesos
	require(ade4)		
	require(spacemakeR)# Para hacer el test.W
	require(lmtest)    # Para las pruebas heterocedasticidad
	require(sp)
	require(rgdal)
	require(MASS)
	require(car)	
	require(nortest)	
	

#==========================================================================
###############              LECTURA DE DATOS               ###############
#==========================================================================

	lattice=readShapePoly("choco.shp") 
	xy=coordinates(lattice)
	names(lattice)

#==========================================================================
###############    ESTADÍSTICAS DESCRIPTIVAS            ###############
#==========================================================================

	 # Estadísticas descriptivas básicas #

			# Función para hallar coeficiente de asimetria #
			skewness=function(x) {
				m3=mean((x-mean(x))^3)
				skew=m3/(sd(x)^3)
				skew}
			# Función para hallar coeficiente Kurtosis #
			kurtosis=function(x) {
				m4=mean((x-mean(x))^4) 
				kurt=m4/(sd(x)^4)-3  
				kurt}
			
			summary(lattice$DEF) # Mínimo, máximo, mediana	y cuartiles 
			sd(lattice$DEF) 	# Desviación estándar
			skewness(lattice$DEF) # Asimetria
			kurtosis(lattice$DEF) # Kurtosis 
			
			summary(lattice$MVM) # Mínimo, máximo, mediana	y cuartiles 
			sd(lattice$MVM) 	# Desviación estándar
			skewness(lattice$MVM) # Asimetria
			kurtosis(lattice$MVM) # Kurtosis 
			
			summary(lattice$AAGS) # Mínimo, máximo, mediana	y cuartiles	
			sd(lattice$AAGS) 	# Desviación estándar
			skewness(lattice$AAGS) # Asimetria
			kurtosis(lattice$AAGS) # Kurtosis 
			
			summary(lattice$CC) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$CC) 	# Desviación estándar
			skewness(lattice$CC) # Asimetria
			kurtosis(lattice$CC) # Kurtosis 
			
			summary(lattice$PO) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$PO) 	# Desviación estándar
			skewness(lattice$PO) # Asimetria
			kurtosis(lattice$PO) # Kurtosis 
			
			summary(lattice$PP) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$PP) 	# Desviación estándar
			skewness(lattice$PP) # Asimetria
			kurtosis(lattice$PP) # Kurtosis 
			
			summary(lattice$PPTN) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$PPTN) 	# Desviación estándar
			skewness(lattice$PPTN) # Asimetria
			kurtosis(lattice$PPTN) # Kurtosis 
			
			summary(lattice$CG) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$CG) 	# Desviación estándar
			skewness(lattice$CG) # Asimetria
			kurtosis(lattice$CG) # Kurtosis 
			
			summary(lattice$CPB) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$CPB) 	# Desviación estándar
			skewness(lattice$CPB) # Asimetria
			kurtosis(lattice$CPB) # Kurtosis 
			
			summary(lattice$NBI) # Mínimo, máximo, mediana	y cuartiles
			sd(lattice$NBI) 	# Desviación estándar
			skewness(lattice$NBI) # Asimetria
			kurtosis(lattice$NBI) # Kurtosis 
			
			summary(lattice$PEND) # Frecuencias

	  # Test de correlación (rho) de Spearman #
			
			cor.test(lattice$MVM,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$AAGS,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$CC,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$PO,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$PP,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$PPTN,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$CG,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$CPB,lattice$DEF,alternative="two.sided",method="spearman")
			cor.test(lattice$NBI,lattice$DEF,alternative="two.sided",method="spearman")

	   # Test no parametrico Kruskal-Wallis (comparación de más de dos medianas entre una variable cuantitativa y una cualitativa) #
			
			kruskal.test(DEF~PEND,data=lattice)

	   # Diagramas de dispersión (relación deforestación con v. exogenas-cuatitativas) 
			
		    op=par(mfrow=c(3,3))	
			
			plot(lattice$MVM,lattice$DEF, ylab="Deforestación (DEF)", xlab="Movilización de madera (MVM)") 
			abline(lm(DEF~MVM,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.32 / P-valor = 0.08") 		

			plot(lattice$AAGS,lattice$DEF, ylab="Deforestación (DEF)", xlab="Área agricola sembrada (AAGS)")
			abline(lm(DEF~AAGS,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.51 / P-valor = 0.004") 

			plot(lattice$CC,lattice$DEF, ylab="Deforestación (DEF)", xlab="Cultivos de Coca (CC)")
			abline(lm(DEF~CC,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.48 / P-valor = 0.007") 

			plot(lattice$PO,lattice$DEF, ylab="Deforestación (DEF)", xlab="Producción de oro (PO)") 
			abline(lm(DEF~PO,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.03 / P-valor = 0.88") 

			plot(lattice$PP,lattice$DEF, ylab="Deforestación (DEF)", xlab="Producción de plata (PP)") 
			abline(lm(DEF~PP,lattice), col="red")
			title("","Corr. Spearman - Rho = -0.15 / P-valor = 0.44") 

			plot(lattice$PPTN,lattice$DEF, ylab="Deforestación (DEF)", xlab="Producción de platino (PPTN)") 
			abline(lm(DEF~PPTN,lattice), col="red")
			title("","Corr. Spearman - Rho = -0.19 / P-valor = 0.31") 

			plot(lattice$CG,lattice$DEF, ylab="Deforestación (DEF)", xlab="Inventario Bovino (CG)") 
			abline(lm(DEF~CG,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.16 / P-valor = 0.39") 

			plot(lattice$CPB,lattice$DEF, ylab="Deforestación (DEF)", xlab="Crecimiento poblacional (CPB)") 
			abline(lm(DEF~CPB,lattice), col="red")
			title("","Corr. Spearman - Rho = 0.11 / P-valor = 0.55") 

			plot(lattice$NBI,lattice$DEF, ylab="Deforestación (DEF)", xlab="Necesidades básicas insatisfechas (NBI)") 
			abline(lm(DEF~NBI,lattice), col="red")
    			title("","Corr. Spearman - Rho = 0.35 / P-valor = 0.06") 

		    par(op)

		# Diagrama de caja bidimensional - comparación de las categorias de pendiente con respecto a la deforestación #

			boxplot(DEF~PEND,data=lattice,ylab="Deforestación (DEF)", xlab="Porcentaje de pendiente (PEND)",col=c("brown1","brown2","brown3","brown4"))
			legend(3.5,12000,c("0-3%","3-7%","12-25%","25-50%"), fill=c("brown1","brown2","brown3","brown4"),title="Pendiente")
			title("","Prueba Kruskal-Wallis - Chi-cuadrado = 1.32 / P-valor = 0.72")
			
		# Pruebas de Normalidad para DEF #
			
			hist(lattice$DEF,col="red",xlab="Deforestación (DEF)",ylab="Frecuencia",main="",BREAKS=7)
			legend("topright",c("Shapiro-Wilk","    W = 0.82  /  P-valor = 0.0001","Kolmogorov-Smirnov","     D = 0.24  /  P-valor = 0.0001"),title="Pruebas de normalidad",bty="o")
			shapiro.test(lattice$DEF) # Ho: La distribución es normal
			lillie.test(lattice$DEF) # Ho: La distribución es normal
			
		# Transformación por Box Cox
			
			bc<-boxcox(DEF~CC+PP+NBI,data=lattice, lambda=seq(-2,2,len=10)) #MVM+AAGS+CC+PO+PP+PPTN+CG+CPB+NBI+PEND2+PEND3+PEND4
			bc
			lambda<-bc$x
			lik<-bc$y
			bco=cbind(lambda,lik)
			bco[order(lik),] # Se selecciona el lambda que tenga mayor valor de logaritmo de verosimilitud (log-likelihood)
			
		# DEF transformada #

			hist(lattice$DEF^(0.34),col="red",xlab="Deforestación (T-DEF)",ylab="Frecuencia",main="",breaks=7)
			legend("topright",c("Shapiro-Wilk","    W = 0.96  /  P-valor = 0.25","Kolmogorov-Smirnov","     D = 0.15  /  P-valor = 0.08"),title="Pruebas de normalidad",bty="0")
			shapiro.test(lattice$DEF^(0.34)) # Ho: La distribución es normal
			lillie.test(lattice$DEF^(0.34)) # Ho: La distribución es normal
			
			
#==========================================================================
###############    ANÁLISIS DE DATOS ESPACIALES            ###############
#==========================================================================

	#==================================================
	###########  ANÁLISIS EXPLORATORIO  ###############
	#==================================================

		#### DISTRIBUCIÓN ESPACIAL ####
		#=============================#
			
		   # boxplot #
		
			boxplot(lattice$MVM,plot=T)
			boxplot(lattice$AAGS,plot=T)
			boxplot(lattice$CC,plot=T)
			boxplot(lattice$PO,plot=T)
			boxplot(lattice$PP,plot=T)
			boxplot(lattice$PPTN,plot=T)
			boxplot(lattice$CG,plot=T)
			boxplot(lattice$CPB,plot=T)
			boxplot(lattice$NBI,plot=T)

 		   # Variable endogena (Distribución y box plot) #
		
			#histomap(lattice,"DEF", col="red", nbcol=7, xlab="Deforestación (DEF)", ylab="Frecuencia") # Histograma regional
			#boxplotmap(lattice,"DEF", col="black")# Boxmap  

		#### ASOCIACIÓN ESPACIAL ####
		#============================#
			
		   # MATRICES DE CONTIGÜIDAD #
		   #:::::::::::::::::::::::::#

		    # CRITERIOS DE CONTIGÜIDAD FISICA PARA DATOS DE ÁREA IRREGULARES #  
			
			# Queen #
			queennb <- poly2nb(lattice,queen=TRUE)
			summary(queennb)
				
			# Rook #
			rooknb <- poly2nb(lattice,queen=FALSE)
			summary(queennb)

			op=par(mfrow=c(1,2))
	
				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(queennb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Reina")
				box()
				
				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(queennb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Rook")
				box()

			par(op)
		    
                # CRITERIOS BASADOS EN GRÁFICAS #

			# Triangulación Delaunay #
			trinb=tri2nb(xy)
			summary(trinb)
			
			# Esferas de influencia #
			soinb=graph2nb(soi.graph(trinb,xy))
			summary(soinb)

			op=par(mfrow=c(1,2))

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(trinb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Triangulación Delaunay")
				box()
				
				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(soinb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Esferas de influencia")
				box()

			 par(op)

		    #	CRITERIOS BASADOS EN DISTANCIAS #		

			# K vecinos más cercanos
			knn1_nb <- knn2nb(knearneigh(xy, k = 1)) # Un vecino más cercano
			knn2_nb <- knn2nb(knearneigh(xy, k = 2)) # Dos vecinos más cercanos
			knn3_nb <- knn2nb(knearneigh(xy, k = 3)) # Tres vecinos más cercanos
			knn4_nb <- knn2nb(knearneigh(xy, k = 4)) # Cuatro vecinos más cercanos
			knn5_nb <- knn2nb(knearneigh(xy, k = 5)) # Cinco vecinos más cercanos
			knn6_nb <- knn2nb(knearneigh(xy, k = 6)) # Seis vecinos más cercanos
			knn7_nb <- knn2nb(knearneigh(xy, k = 7)) # Siete vecinos más cercanos
			knn8_nb <- knn2nb(knearneigh(xy, k = 8)) # Ocho vecinos más cercanos
			knn9_nb <- knn2nb(knearneigh(xy, k = 9)) # Nueve vecinos más cercanos
			knn10_nb <- knn2nb(knearneigh(xy, k = 10)) # Diez vecinos más cercanos
				
			summary(knn1_nb)
			summary(knn2_nb)
			summary(knn3_nb)
			summary(knn4_nb)
			summary(knn5_nb)
			summary(knn6_nb)
			summary(knn7_nb)
			summary(knn8_nb)
			summary(knn9_nb)
			summary(knn10_nb)
      		
			op=par(mfrow=c(2,4))

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn1_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Un vecino")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn2_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Dos vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn3_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Tres vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn4_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Cuatro vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn5_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Cinco vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn6_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Seis vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn7_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Siete vecinos")
				box()

				plot(lattice,border="black",axes=F, cex.axis=0.7)
				plot(knn8_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				title(main="Ocho vecinos")
				box()

				#plot(lattice,border="black",axes=F, cex.axis=0.7)
				#plot(knn9_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				#title(main="Nueve vecinos")

				#plot(lattice,border="black",axes=F, cex.axis=0.7)
				#plot(knn10_nb,coordinates(lattice),add=T,col="red",pch="O",cex=0.6)
				#title(main="Diez vecinos")
   	
			par(op)

		   # SELECCIÓN DE MATRICES DE VECINDAD POR PCNM #
		   #::::::::::::::::::::::::::::::::::::::::::::#

				summary(test.W(lattice$DEF,queennb))
				summary(test.W(lattice$DEF,rooknb)) 
				summary(test.W(lattice$DEF,trinb))    
				summary(test.W(lattice$DEF,soinb))	  
				summary(test.W(lattice$DEF,knn1_nb)) 
				summary(test.W(lattice$DEF,knn2_nb)) 
				summary(test.W(lattice$DEF,knn3_nb))
				summary(test.W(lattice$DEF,knn4_nb))
				summary(test.W(lattice$DEF,knn5_nb))
				summary(test.W(lattice$DEF,knn6_nb))
				summary(test.W(lattice$DEF,knn7_nb)) # Matriz con el valos más bajo de AIC
				summary(test.W(lattice$DEF,knn8_nb)) 
				summary(test.W(lattice$DEF,knn9_nb)) 
				summary(test.W(lattice$DEF,knn10_nb)) 
		         
			 # Se selecciona la matriz de contiguidad basada en distancias (Siete vecinos más cercanos)
 			
		   # CONSTRUCCIÓN MATRIZ DE PESOS ESPACIALES (LISTW) #
		   #:::::::::::::::::::::::::::::::::::::::::::::::::#

			#listw k siete vecinos más cercanos #
			knn7_lw_W <- nb2listw(knn7_nb, style = "W") #Tipo estandarizado por filas
			summary(knn7_lw_W)
			plot(lattice,border="black",axes=F, cex.axis=0.7)
			plot(knn7_lw_W,coordinates(lattice),add=T,col="red",pch=15,cex=0.6)
									
			# Creando archivo .gwt para leer la matriz de pesos espaciales en Geoda o ArcGis #
			triGWT <- tempfile("knn7_lw", tmpdir = "C:/Deforestacion", fileext = ".gwt")
			write.sn2gwt(listw2sn(knn7_lw_W), triGWT)


		# INDICADORES DE AUTOCORRELACIÓN ESPACIAL GLOBAL #
		#================================================#

		   # Test globales con matriz W de siete vecinos más cercanos tipo estandarizado por filas#
			
			moran.test(lattice$DEF, knn7_lw_W) #I de Moran global
			geary.test(lattice$DEF, knn7_lw_W) #C de Geary global
			globalG.test(lattice$DEF, knn7_lw_W) #G global Getis-Ord
			moran.plot(as.vector(scale(lattice$DEF)),knn7_lw_W,xlab="DEF",ylab="Rezago espacial DEF",col="red",pch=16)

		
		# INDICADORES DE AUTOCORRELACIÓN ESPACIAL LOCAL #
		#===============================================#
		
		   # MAPA LISA #

			(locm3=localmoran(lattice$DEF^(0.18),knn7_lw_W))

			lattice$sz <- scale(lattice$DEF^(0.18))
			lattice$lag_sz <- lag.listw(knn7_lw_W,lattice$sz)
			lattice$quad_sig <- NA
			lattice@data[(lattice$sz >= 0 & lattice$lag_sz >= 0) & (locm3[, 5] <= 0.05), "quad_sig"] <- 1
			lattice@data[(lattice$sz <= 0 & lattice$lag_sz <= 0) & (locm3[, 5] <= 0.05), "quad_sig"] <- 2
			lattice@data[(lattice$sz >= 0 & lattice$lag_sz <= 0) & (locm3[, 5] <= 0.05), "quad_sig"] <- 3
			lattice@data[(lattice$sz >= 0 & lattice$lag_sz <= 0) & (locm3[, 5] <= 0.05), "quad_sig"] <- 4
			lattice@data[(lattice$sz <= 0 & lattice$lag_sz >= 0) & (locm3[, 5] <= 0.05), "quad_sig"] <- 5

			breaks=seq(1, 5, 1)
			labels=c("Alto-Alto", "Bajo-Bajo", "Alto-Bajo", "Bajo-Alto", "No Signif.")
			np <- findInterval(lattice$quad_sig, breaks)
			colors <- c("red", "blue", "lightpink", "skyblue2", "white")
			
			par(mar=c(3, 2, 1, 0.5))
			plot(lattice, col = colors[np])#colors[np] establece manualmente el color para cada región
			legend("bottom", legend = labels, ncol=5, fill = colors, bty = "n",cex=0.7)
			#box()

	
	#==================================================
	###########  ANÁLISIS CONFIRMATORIO  ##############
	#==================================================


		# AJUSTE MODELOS LINEALES # 
		#=========================#
			
			# MODELO 1 #
			# MODELO LINEAL CON TODAS LAS VARIABLES #
			chlm <- lm(DEF^(0.34)~MVM+AAGS+CC+PO+PP+PPTN+CG+CPB+NBI+PEND1+PEND2+PEND3,lattice) 
			summary(chlm)
			AIC(chlm)	
			# PRUEBA DE HETEROCEDASTICIDAD #
			bptest(chlm) #No rechaza, Ho: Homocedasticidad en los datos
			# PRUEBA DE AUTOCORRELACIÓN ESPACIAL EN LOS RESIDUALES #
			lm.morantest(chlm, knn7_lw_W) #No es pertinente modelar las relaciones espaciales con el modelo estimado
			# PRUEBAS DE LAGRANGE #
			res <- lm.LMtests(chlm, listw =knn7_lw_W, test = "all")
			tres <- t(sapply(res, function(x) c(x$statistic, x$parameter,x$p.value)))
			colnames(tres) <- c("Statistic", "df", "p-value")
			printCoefmat(tres)
			
			# MODELO LINEAL ESPECÍFICADO #
			chlm2<-lm(DEF^(0.34)~CC+PP+NBI,lattice)
			summary(chlm2)
			AIC(chlm2)
			# PRUEBA DE HETEROCEDASTICIDAD #
			bptest(chlm2) #No rechaza, Ho: Homocedasticidad en los datos
			# PRUEBA DE AUTOCORRELACIÓN ESPACIAL EN LOS RESIDUALES #
			lm.morantest(chlm2, knn7_lw_W) #Si es pertinente modelas las relaciones espaciales con el modelo estimado
			# PRUEBAS DE LAGRANGE #
			res <- lm.LMtests(chlm2, listw =knn7_lw_W, test = "all")
			tres <- t(sapply(res, function(x) c(x$statistic, x$parameter,x$p.value)))
			colnames(tres) <- c("Statistic", "df", "p-value")
			printCoefmat(tres)

				#LMerr: Dependencia espacial residual (dependencia espacial en los errores)
				#LMlag: Dependencia espacial sustantiva (omisión de un retardo espacial de la variable endogena) 
				#RLMerr: Dependencia espacial residual - Prueba robusta de multiplicadores de Lagrange en los errores del modelo 
				#RLMlag: Dependencia espacial sustantiva - Prueba robusta de multiplicadores de Lagrange rezagada
				#SARMA: Contrastación simultanea de autocorrelación espacial residual y sustantiva
			
		# ESTIMACIÓN DE MODELOS CON DEPENDENCIA ESPACIAL #
		#================================================#

			# Estimación para la dependencia espacial sustantiva #
			chlag<- lagsarlm(DEF^(0.34)~CC+PP+NBI,lattice, listw=knn7_lw_W, tol.solve=1.0e-13)
			summary(chlag)	
			bptest.sarlm(chlag)	
			
			# Estimación para la dependencia conjunta sustantiva y residual #
			chsarar<-sacsarlm(DEF^(0.34)~CC+PP+NBI,lattice, listw=knn7_lw_W, type="sac")
			summary(chsarar)

			anova(chlag,chsarar)