roe_pres

pacckages

if(!require(googlesheets4)){install.packages("googlesheets4")}
if(!require(AICcmodavg)){install.packages("AICcmodavg")}
if(!require(ggplot2)){install.packages("ggplot2")}
#
library(googlesheets4); gs4_deauth()
library(AICcmodavg)
library(ggplot2)

base de datos link

sheet="factoriesgo"
range="B4:J400"
roe <- read_sheet("https://docs.google.com/spreadsheets/d/16Hq6UTvxoRsa6aSSo9I0x76IaAXWcfdISXnnGeNXfYw/edit?usp=sharing",
                  col_names = TRUE,
                  sheet=sheet,
                  range=range,
                  col_types = NULL,
                  na= "NA")

## Reading from "Factores de riesgo"

## Range "'factoriesgo'!B4:J400"

dim(roe)

## [1] 396   9

head(roe)

## # A tibble: 6 x 9
##   hh.id piso     pared   rest.alim techo sitio pres.roe pres.dog pres.cat
##   <dbl> <chr>    <chr>   <chr>     <chr> <chr>    <dbl>    <dbl>    <dbl>
## 1     1 A.Tyles  1.Block No        zinc  B            0        0        0
## 2     2 Concrete 1.Block Si        zinc  B            1        1        0
## 3     3 Concrete 1.Block Si        zinc  B            1        1        0
## 4     4 Dirt     adobe   Si        zinc  B            0        1        0
## 5     5 Dirt     adobe   No        zinc  B            0        0        1
## 6     6 Dirt     adobe   Si        zinc  B            1        0        0

Explore data

roe$piso <- as.factor(roe$piso)
roe$pared <- as.factor(roe$pared)
roe$rest.alim <- as.factor(roe$rest.alim)
roe$rest.alim.F <- as.factor(roe$rest.alim)
roe$techo <- as.factor(roe$techo)
roe$sitio <- as.factor(roe$sitio)
roe$pres.roe.F <- as.factor(roe$pres.roe)
roe$pres.dog.F <- as.factor(roe$pres.dog)
roe$pres.cat.F <- as.factor(roe$pres.cat)
#
dim(roe)

## [1] 396  13

names(roe)

##  [1] "hh.id"       "piso"        "pared"       "rest.alim"   "techo"      
##  [6] "sitio"       "pres.roe"    "pres.dog"    "pres.cat"    "rest.alim.F"
## [11] "pres.roe.F"  "pres.dog.F"  "pres.cat.F"

summary(roe)

##      hh.id              piso            pared     rest.alim 
##  Min.   :  1.00   A.Tyles : 42   1.Block   :198   No  : 91  
##  1st Qu.: 99.75   Concrete:183   adobe     : 43   Si  :288  
##  Median :198.50   Dirt    :171   Wood-Cane : 50   NA's: 17  
##  Mean   :198.50                  Wood_other: 54             
##  3rd Qu.:297.25                  zinc_other: 51             
##  Max.   :396.00                                             
##             techo     sitio      pres.roe         pres.dog        pres.cat     
##  conc_zinc_teja: 11   A:148   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  zinc          :385   B: 83   1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:0.0000  
##                       C:165   Median :0.0000   Median :1.000   Median :0.0000  
##                               Mean   :0.3005   Mean   :0.697   Mean   :0.4141  
##                               3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.:1.0000  
##                               Max.   :1.0000   Max.   :1.000   Max.   :1.0000  
##  rest.alim.F pres.roe.F pres.dog.F pres.cat.F
##  No  : 91    0:277      0:120      0:232     
##  Si  :288    1:119      1:276      1:164     
##  NA's: 17                                    
##                                              
##                                              
## 

Plot_explore

## PLOT Sitio ####
ggplot(aes(x= sitio, y= as.factor(pres.roe) , 
           fill= as.factor(pres.roe) ), data= roe ) +
  geom_col() 

ggplot(aes(x= sitio, y= as.factor(pres.roe) , fill= as.factor(pres.roe) ), 
       data= roe ) + 
  geom_point(size= 0.5,  position= position_jitter(width= 0.2 , height =0.1 ))

#### chiqs.test
a <- chisq.test(table(roe$sitio, roe$pres.roe));a; a$stdres

## 
##  Pearson's Chi-squared test
## 
## data:  table(roe$sitio, roe$pres.roe)
## X-squared = 12.589, df = 2, p-value = 0.001847

##    
##             0         1
##   A -1.025216  1.025216
##   B -2.708528  2.708528
##   C  3.242188 -3.242188

## Según la tabla de valores esperados de la Chi2 arriba, los sitios B y C tienen más y menos (respectivamente) proporción de lo "esperado" bajo un modelo nulo.

## PLOT restos alimentos ####
ggplot(aes(x= rest.alim.F, y= as.factor(pres.roe) , 
           fill= as.factor(pres.roe) ), data= roe ) +
  geom_col() 

ggplot(aes(x= rest.alim.F, y= as.factor(pres.roe) , fill= as.factor(pres.roe) ), 
       data= roe ) + 
  geom_point(size= 0.5,  position= position_jitter(width= 0.2 , height =0.1 ))

##  PLOT piso    ####

ggplot(aes(x= piso, y= as.factor(pres.roe) , 
           fill= as.factor(pres.roe) ), data= roe ) +
  geom_col() 

ggplot(aes(x= piso, y= as.factor(pres.roe) , fill= as.factor(pres.roe) ), 
       data= roe ) + 
  geom_point(size= 0.5,  position= position_jitter(width= 0.2 , height =0.1 ))

## PLOT pared ####

ggplot(aes(x= pared, y= as.factor(pres.roe) , 
           fill= as.factor(pres.roe) ), data= roe ) +
  geom_col() 

ggplot(aes(x= pared, y= as.factor(pres.roe) , fill= as.factor(pres.roe) ), 
       data= roe ) + 
  geom_point(size= 0.5,  position= position_jitter(width= 0.2 , height =0.1 ))

# chiqs.test
a <- chisq.test(table(roe$pared, roe$pres.roe));a; a$stdres

## 
##  Pearson's Chi-squared test
## 
## data:  table(roe$pared, roe$pres.roe)
## X-squared = 9.639, df = 4, p-value = 0.04697

##             
##                       0          1
##   1.Block    -0.3288179  0.3288179
##   adobe      -1.7890592  1.7890592
##   Wood-Cane   0.9983159 -0.9983159
##   Wood_other  2.3083089 -2.3083089
##   zinc_other -1.2022716  1.2022716

## PLOT techo ####

ggplot(aes(x= techo, y= as.factor(pres.roe) , 
           fill= as.factor(pres.roe) ), data= roe ) +
  geom_col() 

ggplot(aes(x= techo, y= as.factor(pres.roe) , fill= as.factor(pres.roe) ), 
       data= roe ) + 
  geom_point(size= 0.5,  position= position_jitter(width= 0.2 , height =0.1 ))

#### Asociaci[on entre la presencia de gatos y el sitio 
a <- chisq.test(table(roe$sitio, roe$pres.cat)); a

## 
##  Pearson's Chi-squared test
## 
## data:  table(roe$sitio, roe$pres.cat)
## X-squared = 43.636, df = 2, p-value = 3.347e-10

a$stdres

##    
##             0         1
##   A  3.646605 -3.646605
##   B  3.602760 -3.602760
##   C -6.552865  6.552865

GLM

cor (roe[ , c(7:9)])    ### Chequeando por colinearidad de variables cuantitativas

## Sitio ####
m.sitio <- glm( pres.roe ~ sitio, family= binomial, data=roe)
summary(m.sitio)

#### pres.ROe ~ BUILDING  ####
# piso
m.piso <- glm( pres.roe ~ piso, family= binomial, data=roe)  ### piso
summary(m.piso)
# pared
m.pared <- glm( pres.roe ~ pared, family= binomial, data=roe) ### pared
summary(m.pared)
# techo
m.techo <- glm( pres.roe ~ techo, family= binomial, data=roe) ### techo
summary(m.techo)
# restos de alimento
m.rest.alim <- glm( pres.roe ~ rest.alim, family= binomial, data=roe)   ### rest.alim
summary(m.rest.alim)

m.house <- glm( pres.roe ~ piso + pared + rest.alim + roe$techo, family= binomial, data=roe)
summary(m.house)

### pres.roe ~ DOMESTIC_ANIMAL ####

m.cat <- glm( pres.roe ~ pres.cat, family= binomial, data=roe)    ### gatos
summary(m.cat)

m.dog <- glm( pres.roe ~ pres.dog, family= binomial, data=roe)    ### perros
summary(m.dog)

m.dom.anim <- glm( pres.roe ~ pres.dog * pres.cat, family= binomial, data=roe)
summary(m.dom.anim)

### combinaciones ####
m.pared.cat <- glm( pres.roe ~ pared + pres.cat        , family= binomial, data=roe)
summary(m.pared.cat)
#
m.pared.sitio <- glm( pres.roe ~ pared +            sitio, family= binomial, data=roe)
summary(m.pared.sitio)
#
m.cat.sitio <- glm( pres.roe ~         + pres.cat + sitio, family= binomial, data=roe)
summary(m.cat.sitio)

### MODELO GLOBAL ####
global <- glm( pres.roe ~ pared + pres.cat + sitio, family= binomial, data=roe)
summary(global)

AIC, AICc, Wi

## models

n.= length(roe$pres.roe)
## extracting AICc
extractAIC(global)[2]
k=AICc(global, return.K = TRUE)
AICc.global <- extractAIC(global)[2] + (((2*k)*(k+1))/(n.-k-1))
k1=AICc(global, return.K = TRUE)

extractAIC(m.cat)[2]
k=AICc(m.cat, return.K = TRUE)
AICc.m.cat <- extractAIC(m.cat)[2] + (((2*k)*(k+1))/(n.-k-1))
k2=AICc(m.cat, return.K = TRUE)

extractAIC(m.pared)[2]
k=AICc(m.pared, return.K = TRUE)
AICc.m.pared <- extractAIC(m.pared)[2] + (((2*k)*(k+1))/(n.-k-1))
k3=AICc(m.pared, return.K = TRUE)

extractAIC(m.sitio)[2]
k=AICc(m.sitio, return.K = TRUE)
AICc.m.sitio <- extractAIC(m.sitio)[2] + (((2*k)*(k+1))/(n.-k-1))
k4=AICc(m.sitio, return.K = TRUE)

extractAIC(m.pared.cat)[2]
k=AICc(m.pared.cat, return.K = TRUE)
AICc.m.pared.cat <- extractAIC(m.pared.cat)[2] + (((2*k)*(k+1))/(n.-k-1))
k5=AICc(m.pared.cat, return.K = TRUE)

extractAIC(m.pared.sitio)[2]
k=AICc(m.pared.sitio, return.K = TRUE)
AICc.m.pared.sitio <- extractAIC(m.pared.sitio)[2] + (((2*k)*(k+1))/(n.-k-1))
k6=AICc(m.pared.sitio, return.K = TRUE)

extractAIC(m.cat.sitio)[2]
k=AICc(m.cat.sitio, return.K = TRUE)
AICc.m.cat.sitio <- extractAIC(m.cat.sitio)[2] + (((2*k)*(k+1))/(n.-k-1))
k7=AICc(m.cat.sitio, return.K = TRUE)

########## AIC table #####
aics<- data.frame(paste("m",c("Global", "Gato", "Pared", "Sitio", "pared.gato", "pared.sitio", "gato.sitio"),sep=""),
                  c(AICc.global,AICc.m.cat,AICc.m.pared, AICc.m.sitio, AICc.m.pared.cat, AICc.m.pared.sitio, AICc.m.cat.sitio),
                  c(k1,k2,k3,k4,k5,k6,k7),
                  row.names=NULL) 
colnames(aics) <- c("model","AICc", "k") 
# sort according to AIC
aics<-aics[order(aics$AICc),] 

wi (weights)

for(i in 1:dim(aics)[1]){
  aics$deltaAIC[i]<-aics$AICc[1]-aics$AICc[i]}
wi<-exp(aics$deltaAIC/2)
wi <- wi/sum(wi)
aics$wi <- round(wi/sum(wi), 3)
aics

Delta AICs. Cuánto “pierde” en ajuste (AICc) el modelo completo Gato+Pared+Sitio sacando (ó incluyendo) cada factor?

Deviance Gain Method [ loss of GOF]

###
global$formula; AICc(global)      ###### GATO + PARED 
##########################  Contribution
m.pared.sitio$formula; AICc(m.pared.sitio); AICc(global) - AICc(m.pared.sitio)   ### 'GATO' contribution
m.cat.sitio$formula; AICc(m.cat.sitio); AICc(global)- AICc(m.cat.sitio)          ### 'Pared' contribution
m.pared.cat$formula; AICc(m.pared.cat); AICc(global)- AICc(m.pared.cat)          ### 'Sitio' contribution

#
DeltaAics <- as.data.frame(cbind("Model"= c("global","Gato", "Pared", "Sitio"),
      "Delta.AIC"= c(0,AICc(global) - AICc(m.pared.sitio), AICc(global)- AICc(m.cat.sitio),AICc(global)- AICc(m.pared.cat))))

DeltaAics$Delta.AIC <- as.numeric(DeltaAics$Delta.AIC)
DeltaAics$Delta.AIC <- round(DeltaAics$Delta.AIC,2)
DeltaAics

El incluir pared y Sitio de hecho hace que el modelo empeore (aumenta el AICc)