library(tidyverse)

Objetivo

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.

Datos

"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

Figura 3: Total HDR by pop density

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")

Figura 4: FR HDR by pop density

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")

Figura 5: Total dog pop density by human pop density

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")

Figure 6: FR dog pop density by human pop density

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