library(tidyverse)
Este es el script de las gráficas de los anáilis comparando las diferentes variables de demografía canina y humana en los diferentes sitios del estudio de población canina en Guatemala. Es una modificación basada en el scritp original creado por Laura Grajeda el 16-04-2020.
"DT.csv".
dt<-read.csv("DT.csv")
Revisar variables de densidad poblacional:
dt %>% count(Human.Population.Density, SITE)
dt$Human.Population.Density <- as.numeric(gsub(",", "", dt$Human.Population.Density)) # quitar las "," y darle fomato numerico
Revisar variables a usar en figura 3:
#dt %>% count(REVISED.TOTAL.HDR, SITE)
#dt %>% count(Revised.LO.1, SITE)
#dt %>% count(Revised.HI.1, SITE)
dt %>% select(Human.Population.Density, REVISED.TOTAL.HDR, Revised.LO.1, Revised.HI.1)
Calcular Linea de tendencia para REVISED.TOTAL.HDR:
m <-lm(log(REVISED.TOTAL.HDR) ~ log(Human.Population.Density), data=dt)
summary(m)
##
## Call:
## lm(formula = log(REVISED.TOTAL.HDR) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61196 -0.20391 0.03631 0.13873 0.65088
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.39206 0.63177 -2.203 0.044811 *
## log(Human.Population.Density) 0.37976 0.09015 4.213 0.000869 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3834 on 14 degrees of freedom
## Multiple R-squared: 0.559, Adjusted R-squared: 0.5275
## F-statistic: 17.75 on 1 and 14 DF, p-value: 0.0008689
exp(m$coefficients[1]) # Esto es a
## (Intercept)
## 0.2485625
m$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.379759
summary(m)$adj.r.squared # Esto es R2
## [1] 0.5275083
# Crear funcion que luego sera utilizada en el grafico
fun <- function(x) exp(m$coefficients[1])* x^m$coefficients[2]
Calcular Linea de tendencia para Revised.LO.1:
mLO <-lm(log(Revised.LO.1) ~ log(Human.Population.Density), data=dt)
summary(mLO)
##
## Call:
## lm(formula = log(Revised.LO.1) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69382 -0.23849 0.05639 0.19143 0.75108
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8178 0.7103 -1.151 0.26885
## log(Human.Population.Density) 0.3293 0.1014 3.249 0.00582 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4311 on 14 degrees of freedom
## Multiple R-squared: 0.4299, Adjusted R-squared: 0.3892
## F-statistic: 10.56 on 1 and 14 DF, p-value: 0.005823
exp(mLO$coefficients[1]) # Esto es a
## (Intercept)
## 0.4413855
mLO$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.3293033
# Crear funcion que luego sera utilizada en el grafico
funLO <- function(x) exp(mLO$coefficients[1])* x^mLO$coefficients[2]
Calcular linea de tendencia para Revised.HI.1:
mHI <-lm(log(Revised.HI.1) ~ log(Human.Population.Density), data=dt)
summary(mHI)
##
## Call:
## lm(formula = log(Revised.HI.1) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54877 -0.27509 0.03142 0.14588 0.59771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.81010 0.60526 -2.991 0.00973 **
## log(Human.Population.Density) 0.41502 0.08636 4.805 0.00028 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3674 on 14 degrees of freedom
## Multiple R-squared: 0.6226, Adjusted R-squared: 0.5956
## F-statistic: 23.09 on 1 and 14 DF, p-value: 0.0002797
exp(mHI$coefficients[1]) # Esto es a
## (Intercept)
## 0.1636372
mHI$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.4150229
# Crear funcion que luego sera utilizada en el grafico
funHI <- function(x) exp(mHI$coefficients[1])* x^mHI$coefficients[2]
Hacer figura:
Figura3<-ggplot(dt, aes(x=Human.Population.Density))+
geom_pointrange(aes(ymin = Revised.LO.1, ymax = Revised.HI.1, y=REVISED.TOTAL.HDR))+
stat_function(fun = fun)+
stat_function(fun = funLO, colour="gray")+
stat_function(fun = funHI, colour="gray") +
geom_text(x = 1400, y = 12,
label = as.expression(bquote(~y==.(as.character(round(exp(m$coefficients[1]),4)))
~"x"^{.(as.character(round(m$coefficients[2],4)))}~", "~
R^2==.(round(summary(m)$adj.r.squared,4))))) +
labs(x="Human population density", y="Free Roaming Human dog ratio") +
scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14))+
scale_x_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000))+
theme_light()
Figura3
ggsave(plot=Figura3, "Figure 3_TOTAL HDR by POP DENSITY.png")
Revisar variables que se usaran en la figura 4:
#dt %>% count(Free.Roaming.HDR, SITE)
#dt %>% count(FR.HDR.95..Low, SITE)
dt$Free.Roaming.HDR_LO <- dt$FR.HDR.95..High # High y low estan intercambiadas
#dt %>% count(FR.HDR.95..High, SITE)
dt$Free.Roaming.HDR_HI <- dt$FR.HDR.95..Low # High y low estan intercambiadas
dt %>% select(Human.Population.Density, Free.Roaming.HDR_LO, Free.Roaming.HDR, Free.Roaming.HDR_HI)
Calcular linea de tendencia para Free.Roaming.HDR:
m <-lm(log(Free.Roaming.HDR) ~ log(Human.Population.Density), data=dt)
summary(m)
##
## Call:
## lm(formula = log(Free.Roaming.HDR) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61911 -0.21081 -0.00767 0.20189 0.64153
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.25559 0.60042 -2.091 0.055227 .
## log(Human.Population.Density) 0.39759 0.08567 4.641 0.000382 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3644 on 14 degrees of freedom
## Multiple R-squared: 0.606, Adjusted R-squared: 0.5779
## F-statistic: 21.54 on 1 and 14 DF, p-value: 0.0003818
exp(m$coefficients[1]) # Esto es a
## (Intercept)
## 0.284909
m$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.3975875
summary(m)$adj.r.squared # Esto es R2
## [1] 0.5778932
# Crear funcion que luego sera utilizada en el grafico
fun <- function(x) exp(m$coefficients[1])* x^m$coefficients[2]
Calcular linea de tendencia para Free.Roaming.HDR_LO:
mLO <-lm(log(Free.Roaming.HDR_LO) ~ log(Human.Population.Density), data=dt)
summary(mLO)
##
## Call:
## lm(formula = log(Free.Roaming.HDR_LO) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59031 -0.17313 -0.02902 0.12983 0.54787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.4851 0.5551 -2.676 0.018103 *
## log(Human.Population.Density) 0.4115 0.0792 5.196 0.000136 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3369 on 14 degrees of freedom
## Multiple R-squared: 0.6585, Adjusted R-squared: 0.6341
## F-statistic: 26.99 on 1 and 14 DF, p-value: 0.0001357
exp(mLO$coefficients[1]) # Esto es a
## (Intercept)
## 0.2264873
mLO$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.4114849
# Crear funcion que luego sera utilizada en el grafico
funLO <- function(x) exp(mLO$coefficients[1])* x^mLO$coefficients[2]
Calcular linea de tendencia para Free.Roaming.HDR_HI:
mHI <-lm(log(Free.Roaming.HDR_HI) ~ log(Human.Population.Density), data=dt)
summary(mHI)
##
## Call:
## lm(formula = log(Free.Roaming.HDR_HI) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68985 -0.31158 -0.04257 0.32772 0.73867
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.83436 0.66683 -1.251 0.23135
## log(Human.Population.Density) 0.36689 0.09515 3.856 0.00175 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4047 on 14 degrees of freedom
## Multiple R-squared: 0.515, Adjusted R-squared: 0.4804
## F-statistic: 14.87 on 1 and 14 DF, p-value: 0.001747
exp(mHI$coefficients[1]) # Esto es a
## (Intercept)
## 0.4341514
mHI$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.3668875
# Crear funcion que luego sera utilizada en el grafico
funHI <- function(x) exp(mHI$coefficients[1])* x^mHI$coefficients[2]
Hacer figura:
Figura4<-ggplot(dt, aes(x=Human.Population.Density))+
geom_pointrange(aes(ymin = Free.Roaming.HDR_LO,
ymax = Free.Roaming.HDR_HI,
y=Free.Roaming.HDR))+
stat_function(fun = fun)+
stat_function(fun = funLO, colour="gray")+
stat_function(fun = funHI, colour="gray") +
geom_text(x = 1400, y = 16,
label = as.expression(bquote(~y==.(as.character(round(exp(m$coefficients[1]),4)))~"x"^{.(as.character(round(m$coefficients[2],4)))}~", "~R^2==.(round(summary(m)$adj.r.squared,4))))) +
labs(x="Human population density", y="Human dog ratio") +
scale_y_continuous(breaks = c(0,2,4,6,8,10,12,14,16,18))+
scale_x_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000))+
theme_light()
Figura4
ggsave(plot=Figura4, "Figure 4_FR HDR by POP DENSITY.png")
Revisar variables que seran usados en la figura 5:
#dt %>% count(Dogs.Per.KM2, SITE)
dt$Dogs.Per.KM2 <- as.numeric(gsub(",", "", dt$Dogs.Per.KM2)) # quitar las "," y darle fomato numerico
#dt %>% count(Dogs.per.KM2.Lo, SITE)
#dt %>% count(Dogs.per.KM2.Hi, SITE)
dt$Dogs.per.KM2.Hi <- as.numeric(gsub(",", "", dt$Dogs.per.KM2.Hi)) # quitar las "," y darle fomato numerico
dt %>% select(Human.Population.Density, Dogs.per.KM2.Lo, Dogs.Per.KM2, Dogs.per.KM2.Hi)
Calcular linea de tendencia para Dogs.Per.KM2:
m <-lm(log(Dogs.Per.KM2) ~ log(Human.Population.Density), data=dt)
summary(m)
##
## Call:
## lm(formula = log(Dogs.Per.KM2) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64959 -0.13543 -0.03765 0.20089 0.60453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.37488 0.62606 2.196 0.0454 *
## log(Human.Population.Density) 0.62239 0.08933 6.967 6.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.38 on 14 degrees of freedom
## Multiple R-squared: 0.7761, Adjusted R-squared: 0.7602
## F-statistic: 48.54 on 1 and 14 DF, p-value: 6.582e-06
exp(m$coefficients[1]) # Esto es a
## (Intercept)
## 3.954595
m$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.6223888
summary(m)$adj.r.squared # Esto es R2
## [1] 0.7601551
# Crear funcion que luego sera utilizada en el grafico
fun <- function(x) exp(m$coefficients[1])* x^m$coefficients[2]
Calcular linea de tendencia para Dogs.Per.KM2.Lo:
mLO <-lm(log(Dogs.per.KM2.Lo) ~ log(Human.Population.Density), data=dt)
summary(mLO)
##
## Call:
## lm(formula = log(Dogs.per.KM2.Lo) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74953 -0.18248 -0.06224 0.23588 0.69653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8067 0.7113 1.134 0.276
## log(Human.Population.Density) 0.6719 0.1015 6.620 1.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4317 on 14 degrees of freedom
## Multiple R-squared: 0.7579, Adjusted R-squared: 0.7406
## F-statistic: 43.83 on 1 and 14 DF, p-value: 1.15e-05
exp(mLO$coefficients[1]) # Esto es a
## (Intercept)
## 2.240487
mLO$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.6718976
# Crear funcion que luego sera utilizada en el grafico
funLO <- function(x) exp(mLO$coefficients[1])* x^mLO$coefficients[2]
Calcular linea de tendencia para Dogs.per.KM2.Lo.Hi:
mHI <-lm(log(Dogs.per.KM2.Hi) ~ log(Human.Population.Density), data=dt)
summary(mHI)
##
## Call:
## lm(formula = log(Dogs.per.KM2.Hi) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59940 -0.15214 -0.02143 0.26558 0.55576
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7837 0.6062 2.943 0.0107 *
## log(Human.Population.Density) 0.5886 0.0865 6.805 8.53e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3679 on 14 degrees of freedom
## Multiple R-squared: 0.7679, Adjusted R-squared: 0.7513
## F-statistic: 46.31 on 1 and 14 DF, p-value: 8.53e-06
exp(mHI$coefficients[1]) # Esto es a
## (Intercept)
## 5.952119
mHI$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.5886015
# Crear funcion que luego sera utilizada en el grafico
funHI <- function(x) exp(mHI$coefficients[1])* x^mHI$coefficients[2]
Hacer figura:
Figura5<-ggplot(dt, aes(x=Human.Population.Density))+
geom_pointrange(aes(ymin = Dogs.per.KM2.Lo,
ymax = Dogs.per.KM2.Hi,
y=Dogs.Per.KM2))+
stat_function(fun = fun)+
stat_function(fun = funLO, colour="gray")+
stat_function(fun = funHI, colour="gray") +
geom_text(x = 1400, y = 1150,
label = as.expression(bquote(~y==.(as.character(round(exp(m$coefficients[1]),4)))~"x"^{.(as.character(round(m$coefficients[2],4)))}~", "~R^2==.(round(summary(m)$adj.r.squared,4))))) +
labs(x="Human population density", y=expression(paste("Dogs per Km"^2))) +
scale_y_continuous(breaks = c(0,200,400,600,800,1000,1200))+
scale_x_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000))+
theme_light()
Figura5
ggsave(plot=Figura5, "Figure 5_Total dog pop density by human pop density.png")
Revisar variables que seran usadas en la figura 6:
#dt %>% count(FR.Dogs.per.KM2, SITE)
#dt %>% count(FR.Dogs.perKM2.Lo, SITE)
#dt %>% count(FR.Dogs.per.KM2.Hi, SITE)
dt %>% select(Human.Population.Density, FR.Dogs.perKM2.Lo, FR.Dogs.per.KM2, FR.Dogs.per.KM2.Hi)
Calcular linea de tendencia para FR.Dogs.per.KM2:
m <-lm(log(FR.Dogs.per.KM2) ~ log(Human.Population.Density), data=dt)
summary(m)
##
## Call:
## lm(formula = log(FR.Dogs.per.KM2) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64275 -0.21511 0.00905 0.21427 0.60416
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.24469 0.59926 2.077 0.0567 .
## log(Human.Population.Density) 0.60393 0.08551 7.063 5.66e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3637 on 14 degrees of freedom
## Multiple R-squared: 0.7809, Adjusted R-squared: 0.7652
## F-statistic: 49.88 on 1 and 14 DF, p-value: 5.658e-06
exp(m$coefficients[1]) # Esto es a
## (Intercept)
## 3.471852
m$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.6039273
summary(m)$adj.r.squared # Esto es R2
## [1] 0.7651978
# Crear funcion que luego sera utilizada en el grafico
fun <- function(x) exp(m$coefficients[1])* x^m$coefficients[2]
Calcular linea de tendencia para FR.Dogs.perKM2.Lo:
mLO <-lm(log(FR.Dogs.perKM2.Lo) ~ log(Human.Population.Density), data=dt)
summary(mLO)
##
## Call:
## lm(formula = log(FR.Dogs.perKM2.Lo) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73774 -0.32615 0.04776 0.30250 0.67930
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.84196 0.66007 1.276 0.223
## log(Human.Population.Density) 0.63210 0.09418 6.711 9.92e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4006 on 14 degrees of freedom
## Multiple R-squared: 0.7629, Adjusted R-squared: 0.7459
## F-statistic: 45.04 on 1 and 14 DF, p-value: 9.923e-06
exp(mLO$coefficients[1]) # Esto es a
## (Intercept)
## 2.32091
mLO$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.6320976
# Crear funcion que luego sera utilizada en el grafico
funLO <- function(x) exp(mLO$coefficients[1])* x^mLO$coefficients[2]
Calcular linea de tendencia para Dogs.per.KM2.Lo.Hi:
mHI <-lm(log(FR.Dogs.per.KM2.Hi) ~ log(Human.Population.Density), data=dt)
summary(mHI)
##
## Call:
## lm(formula = log(FR.Dogs.per.KM2.Hi) ~ log(Human.Population.Density),
## data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54613 -0.13423 0.04103 0.16812 0.58286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.50466 0.55212 2.725 0.0164 *
## log(Human.Population.Density) 0.58565 0.07878 7.434 3.18e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3351 on 14 degrees of freedom
## Multiple R-squared: 0.7979, Adjusted R-squared: 0.7834
## F-statistic: 55.26 on 1 and 14 DF, p-value: 3.184e-06
exp(mHI$coefficients[1]) # Esto es a
## (Intercept)
## 4.502608
mHI$coefficients[2] # Esto es b
## log(Human.Population.Density)
## 0.5856474
# Crear funcion que luego sera utilizada en el grafico
funHI <- function(x) exp(mHI$coefficients[1])* x^mHI$coefficients[2]
Crear figura:
Figura6<-ggplot(dt, aes(x=Human.Population.Density))+
geom_pointrange(aes(ymin = FR.Dogs.perKM2.Lo,
ymax = FR.Dogs.per.KM2.Hi,
y=FR.Dogs.per.KM2))+
stat_function(fun = fun)+
stat_function(fun = funLO, colour="gray")+
stat_function(fun = funHI, colour="gray") +
geom_text(x = 1400, y = 900,
label = as.expression(bquote(~y==.(as.character(round(exp(m$coefficients[1]),4)))~"x"^{.(as.character(round(m$coefficients[2],4)))}~", "~R^2==.(round(summary(m)$adj.r.squared,4))))) +
labs(x="Human population density", y=expression(paste("Free roaming dogs per Km"^2))) +
scale_y_continuous(breaks = c(0,200,400,600,800,1000,1200))+
scale_x_continuous(breaks = c(0, 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000))+
theme_light()
Figura6
ggsave(plot=Figura6, "Figure 6_FR dog pop density by human pop density.png")
FIN