# Libraries

library(stats)

library(rgdal)

library(Rcpp)

library(sf)

library(pwr)

library(gvlma)

library(lmtest)

library(dplyr)

library(ggplot2)

library(stargazer)

library(pastecs)

library(psych)

library(sjPlot)

library(sjmisc)

library(sjlabelled)

library(bbmle)

library(ggpubr)

library(fitdistrplus)

library(rcompanion)

library(mgcv)

library(stats)

library(lmtest)

library(jtools)

library(car)

library(drc)

library(nlme)

library(remotes)

# Datos y mapa

cursomodelos.shape <- readOGR(dsn = "...", layer = "X")

plot(cursomodelos.shape)

names(cursomodelos.shape)

summary(cursomodelos.shape)

# Descriptivos

## Variables

### Sacar área

cursomodelos.shape <- st_read("...")

cursomodelos.shape$area <- st_area(cursomodelos.shape)

names(cursomodelos.shape)

View(cursomodelos.shape)

area = cursomodelos.shape$area

sum(area)

cursomodelos.shape$area2 <- cursomodelos.shape$area/1000000

sum(area2)

attach(cursomodelos.shape)

### Definir variables de análisis

cursomodelos.shape$pobtotal <- cursomodelos.shape$POBTOT_sum

cursomodelos.shape$tasahoms <- cursomodelos.shape$tothoms/cursomodelos.shape$pobtotal *1000

cursomodelos.shape$densidad <- cursomodelos.shape$pobtotal/cursomodelos.shape$area2

cursomodelos.shape$hacina <- cursomodelos.shape$OCUPVIVPAR/cursomodelos.shape$VIVPAR_HAB

cursomodelos.shape$migrantes <- cursomodelos.shape$PRESOE15_s/cursomodelos.shape$pobtotal

tasahoms = tasahoms[is.na(tasahoms)] <- 0

hacina = hacina[is.na(hacina)] <- 0

pobtotal = as.numeric(cursomodelos.shape$pobtotal)

tasahoms = as.numeric(cursomodelos.shape$tasahoms)

densidad = as.numeric(cursomodelos.shape$densidad)

hacina = as.numeric(cursomodelos.shape$hacina)

migrantes = as.numeric(cursomodelos.shape$migrantes)

summary(tothoms)

summary(pobtotal)

summary(tasahoms)

summary(densidad)

summary(hacina)

summary(migrantes)

summary(tasahoms)

describe(tasahoms)

summary(hacina)

describe(hacina)

names(cursomodelos.shape)

# Mapear homicidios

plot(cursomodelos.shape["tasahoms"])

ggplot(cursomodelos.shape) +

geom_sf(aes(fill=tasahoms))

histograma = hist(tasahoms,xlab = "tasahoms",col = "red",border = "black", breaks = 10)

# Modelo 1: OLS

reg.eq1=tasahoms ~ scale(densidad) + scale(hacina) + scale(migrantes)

reg1=lm(reg.eq1,data=cursomodelos.shape)

summary(reg1)

tab_model(reg1, p.style = "stars", show.aic = T)

stargazer(reg1, title="Regression Results", align=TRUE)

(est <- cbind(Estimate = coef(reg1), confint(reg1)))

## Diagnósticos

ols.res1 = resid(reg1)

hist(ols.res1)

qqnorm(ols.res1);qqline(ols.res1)

reg1 <- gvlma(reg1)

summary(reg1)

plot.gvlma(reg1)

mean(fitted(reg1))

hist(fitted(reg1), main="Reg1: Valores Predichos de tasa de Homicidio")

plot(lm(tasahoms ~ scale(densidad) + scale(hacina) + scale(migrantes),

data=cursomodelos.shape))

## Heterocedasticidad residuales

bptest(reg1,studentize=TRUE)

## Normalidad residuales

normality.test <- shapiro.test(ols.res1)

print(normality.test)

## Colinealidad

car::vif(reg1)

## Prueba de Ramsey

Ramsey2 <- resettest(reg1, power = 2, type = "fitted")

Ramsey2

Ramsey3 <- resettest(reg1, power = 3, type = "fitted")

Ramsey3

## Exploracion de formas funcionales

descdist(tasahoms)

plotdist(tasahoms, histo=TRUE, demp=TRUE)

# Modelo 2: Log-log (Gamma)

describe(tasahoms)

describe(densidad)

describe(hacina)

describe(migrantes)

tasahomslog = log(1+tasahoms)

densidadlog = log(1+densidad)

hacinalog = log(1+hacina)

migranteslog = log(1+migrantes)

descdist(tasahomslog)

plotdist(tasahoms, histo=TRUE, demp=TRUE)

reg.eq2=tasahomslog ~ densidadlog + hacinalog + migranteslog

reg2=lm(reg.eq2,data=cursomodelos.shape)

summary(reg2)

tab_model(reg1, reg2, p.style = "stars", show.aic = T)

stargazer(reg1, reg2, title="Regression Results", align=TRUE, type="text")

stargazer(reg1, reg2, title="Regression Results", align=TRUE, type="latex")

AIC(reg1, reg2)

(est <- cbind(Estimate = coef(reg2), confint(reg2)))

## Diagnósticos

ols.res2 = resid(reg2)

hist(ols.res2)

qqnorm(ols.res2);qqline(ols.res4)

reg2 <- gvlma(reg2)

summary(reg2)

plot.gvlma(reg2)

mean(fitted(reg2))

hist(fitted(reg2), main="Log-log: Valores Predichos de tasa de Homicidio")

## Heterocedasticidad

bptest(reg2,studentize=TRUE)

## Normalidad residuales

normality.test <- shapiro.test(ols.res2)

print(normality.test)

### A mayor detalle

ggqqplot(tasahomslog, ylab = "tasahomslog")

plot(densidadlog, tasahomslog)

abline(fit, col = "red")

AIC(reg1, reg2)

cor.test(densidadlog,tasahomslog, method = c("pearson"))

cor.test(hacinalog,tasahomslog, method = c("pearson"))

cor.test(migranteslog,tasahomslog, method = c("pearson"))

# Modelo 3: Exponencial

tasahomslog = log(1+tasahoms)

descdist(tasahomslog)

plotdist(tasahoms, histo=TRUE, demp=TRUE)

reg.eq3=tasahomslog ~ densidad + hacina + migrantes

reg3=lm(reg.eq3,data=cursomodelos.shape)

summary(reg3)

tab_model(reg1, reg2, reg3, p.style = "stars", show.aic = T)

stargazer(reg1, reg2, title="Regression Results", align=TRUE, type="text")

stargazer(reg1, reg2, title="Regression Results", align=TRUE)

## Diagnósticos

reg3 <- gvlma(reg3)

summary(reg3)

plot.gvlma(reg3)

plot(lm(tasahomslog ~ densidad + hacina + migrantes,

data=cursomodelos.shape))

# Modelo 4: Log-normal con interacciones

reg.eq4= tasahomslog ~ log(pobtotal) + log(densidad) + (log(pobtotal)*log(densidad))

reg4=lm(reg.eq4,data=cursomodelos.shape)

summary(reg4)

tab_model(reg4, p.style = "stars")

AIC(reg1, reg2, reg3, reg4)

## Diagnósticos

reg4 <- gvlma(reg4)

summary(reg4)

plot.gvlma(reg4)