11 So What

Demonstrate that different georeferencing decisions will produce different results in a simple linear regression model in terms of both statistical significance and substantive effects.

rm(list=ls()); gc()
# !diagnostics off
library(MeasuringLandscape)
dir_figures <- glue::glue(getwd(), "/../paper/figures/")
gc()
knitr::opts_knit$set(progress = TRUE, verbose = TRUE)
knitr::opts_chunk$set(fig.width=12, fig.height=8,  warning=FALSE, message=FALSE, cache=TRUE)
options(width = 160)

Create Hexagon Unit of Analysis

covariate_list <-  MeasuringLandscape:::prep_covariates()
region_of_interest_sf <- MeasuringLandscape:::create_roi(bottom_left_x=35.67,
                                        bottom_left_y=-1.43285,
                                        top_right_x=38.19,
                                        top_right_y=0.54543,
                                        crs_out=4326)
hex <- MeasuringLandscape:::create_hexagon_df(cellsize_km=10,
                            region_of_interest=region_of_interest_sf,
                            type="hexagonal",
                            crs=4326 )

[1] 10

Measure covariates at each hexagon

#hex$district <- new_over(st_centroid(hex) , covariate_list[[1]]  , 'name' )
#hex$cadastral <- new_over(st_centroid(hex) , covariate_list[[2]]  , 'name' )
#hex$language <- new_over(st_centroid(hex) , covariate_list[[3]]  , 'LANGUAGE' )
#hex$tribe <- new_over(st_centroid(hex) , covariate_list[[4]]  , 'Tribe' )
hex$rain <- log( MeasuringLandscape:::new_over(st_centroid(hex) , covariate_list[['raster_rain']]  , 'Tribe' ) + 1)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infst_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

hex$population <- log(  MeasuringLandscape:::new_over(st_centroid(hex) , covariate_list[['pop_raster_roi']]  , '' ) + 1)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infst_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

hex$treecover <-log(  MeasuringLandscape:::new_over(st_centroid(hex) , covariate_list[['forest_raster_roi']]  , '' ) + 1)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infst_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

hex$ruggedness <- log(  MeasuringLandscape:::new_over(st_centroid(hex) , covariate_list[['ruggedness_raster_roi']]  , '' ) + 1)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infst_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

hex$roads_distance <-log(  MeasuringLandscape:::new_over(st_centroid(hex) , covariate_list[['roads_distance_to']]  , '' ) + 1)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infst_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

#hex$landuse <- new_over(st_centroid(hex) , covariate_list[[9]]  , 'LANDUSE' )

Tally number of events in each hexagon according to each georeferencing rule and fit a linear model

#Reload from scratch each time in case we subset sometehing weirdly
georef_all_dt <- readRDS(system.file("extdata", "georef_all_dt_recomendations.Rds", package = "MeasuringLandscape")) 
georef_all_dt[!is.finite(distance_km),distance_km:=NA] 
georef_all_dt <- subset(georef_all_dt,  distance_km=!0) #This excludes all self references , we need that here but not necessarily elsewhere
georef_all_dt_events_sf <- georef_all_dt  %>% 
                         filter(!is.na(X1) & !is.na(Y2)) %>% #can drop anything without a coordinate because you'll never be able to aggregate those
                         data.frame() %>% st_as_sf(coords = c("X2","Y2"),  crs = 4326, agr = "constant", remove=F, na.fail =F)
temp <- georef_all_dt_events_sf %>% filter(fuzzy) %>% filter(!duplicated(event_hash))
hex$dv <- log( sapply( st_covers(hex, temp), sum, na.rm=T) + 1)

although coordinates are longitude/latitude, st_covers assumes that they are planar

vars_all <- c("dv","rain","population","treecover","ruggedness","roads_distance") #,"landuse"
#lapply(hex[,vars_all], hist)
lm1 <- lm(dv ~ ., data=as.data.frame(hex)[,vars_all])
temp <- georef_all_dt_events_sf %>% filter(!duplicated(event_hash))
hex$dv <-  sapply( st_covers(hex, temp), length) 

although coordinates are longitude/latitude, st_covers assumes that they are planar

broom::tidy(lm1)

Fit Linear Models using Each Possible Georeferencing Decision

For each possible georeferencing decision, we fit a seperate linear model predicting the number of events in an area as a function of five covariates, population, rainfall, distance from roads, ruggedness, and treecover. The unit of obsevration is a hexagonal bin of size 10km. The model is a negative binomial regression.

#Hand Rule
lm_results_list0 <- list()
data.table::setkey(georef_all_dt, rule_ensemble)
georef_all_dt_ensemble <- georef_all_dt[,.SD[1], by=list(event_hash) ]
condition <- paste(georef_all_dt_events_sf$event_hash, georef_all_dt_events_sf$place_hash) %in% 
             paste(georef_all_dt_ensemble$event_hash, georef_all_dt_ensemble$place_hash) #table(condition)
temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
hex$dv <- sapply( st_covers(hex, temp), length) 

although coordinates are longitude/latitude, st_covers assumes that they are planar

lm1 <- MASS::glm.nb(dv ~ ., data= na.omit(as.data.frame(hex)[,vars_all]) )
temp_results <- broom::tidy(lm1)
temp_results$Type="Rule"
temp_results$label <- "Hand Rule"
lm_results_list0[[as.character("Hand Rule")]] <- temp_results
#Ensemble Rule
data.table::setkey(georef_all_dt, rule_ensemble)
georef_all_dt_ensemble <- georef_all_dt[,.SD[1], by=list(event_hash) ]
condition <- paste(georef_all_dt_events_sf$event_hash, georef_all_dt_events_sf$place_hash) %in% 
             paste(georef_all_dt_ensemble$event_hash, georef_all_dt_ensemble$place_hash) #table(condition)
temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
hex$dv <- sapply( st_covers(hex, temp), length) 

although coordinates are longitude/latitude, st_covers assumes that they are planar

lm1 <- MASS::glm.nb(dv ~ ., data= na.omit(as.data.frame(hex)[,vars_all]) )
temp_results <- broom::tidy(lm1)
temp_results$Type="Rule"
temp_results$label <- "Ensemble Rule"
lm_results_list0[[as.character("Ensemble Rule")]] <- temp_results
  
#Source Dataset
lm_results_list <- list()
for(q in unique(georef_all_dt_events_sf$source_dataset)){
    print(q)
  try({ #TGN fails
  condition <- georef_all_dt_events_sf$source_dataset==q
  temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
  hex$dv <- sapply( st_covers(hex, temp), length) 
  lm1 <- MASS::glm.nb(dv ~ .,
                data= na.omit(as.data.frame(hex)[,vars_all]) 
                )
  temp_results <- broom::tidy(lm1)
  temp_results$Type="Source Dataset"
  temp_results$label <- q
  lm_results_list[[as.character(q)]] <- temp_results
  })
}

[1] "events"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "historical"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "nga"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "geonames"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "livestock_points"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "bing"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "wikidata"

although coordinates are longitude/latitude, st_covers assumes that they are planar
glm.fit: algorithm did not convergealternation limit reached

[1] "tgn"

although coordinates are longitude/latitude, st_covers assumes that they are planar
step size truncated due to divergenceError in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  : 
  NA/NaN/Inf in 'x'

[1] "google"

although coordinates are longitude/latitude, st_covers assumes that they are planar
glm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not convergealternation limit reached

[1] "openstreetmap"

although coordinates are longitude/latitude, st_covers assumes that they are planar
step size truncated due to divergenceError in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  : 
  NA/NaN/Inf in 'x'

[1] "gadm"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "kenya_cadastral_district"

although coordinates are longitude/latitude, st_covers assumes that they are planar
glm.fit: algorithm did not convergeglm.fit: algorithm did not convergeglm.fit: algorithm did not converge

[1] "kenya_district1962"

although coordinates are longitude/latitude, st_covers assumes that they are planar
glm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergeglm.fit: algorithm did not convergeiteration limit reachedNaNs producediteration limit reachedNaNs producedglm.fit: algorithm did not convergealternation limit reached

[1] "livestock_boundaries"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "kenya_cadastral"

although coordinates are longitude/latitude, st_covers assumes that they are planar
Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  : 
  NA/NaN/Inf in 'x'

#Geometry Type
#Half the time, linestring doesn't converge. Exclude from the plots downstream anyway
lm_results_list2 <- list()
for(q in unique(georef_all_dt_events_sf$geometry_type)){
  print(q)
  try({
    condition <- georef_all_dt_events_sf$geometry_type==q
    temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
    hex$dv <- sapply( st_covers(hex, temp), length) 
    lm1 <- MASS::glm.nb(dv ~ ., data= na.omit(as.data.frame(hex)[,vars_all]) ,
                  control=glm.control(maxit=100) #Up the number of iterations so it converges
                  ) 
    temp_results <- broom::tidy(lm1)
    temp_results$Type="Geometry Type"
    temp_results$label <- q
    lm_results_list2[[as.character(q)]] <- temp_results
  })
}

[1] "POINT"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "MULTIPOLYGON"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "POLYGON"

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] "LINESTRING"

although coordinates are longitude/latitude, st_covers assumes that they are planar
step size truncated due to divergenceError in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  : 
  NA/NaN/Inf in 'x'

#Fuzzy
lm_results_list3 <- list()
for(q in unique(georef_all_dt_events_sf$fuzzy)){
  print(q)
  condition <- georef_all_dt_events_sf$fuzzy==q
  temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
  hex$dv <- sapply( st_covers(hex, temp), length) 
  lm1 <- MASS::glm.nb(dv ~ ., data= na.omit(as.data.frame(hex)[,vars_all]) )
  temp_results <- broom::tidy(lm1)
  temp_results$Type="Match Type"
  temp_results$label <- q
  lm_results_list3[[as.character(q)]] <- temp_results
}

[1] FALSE

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] TRUE

although coordinates are longitude/latitude, st_covers assumes that they are planar

#Self Reference
lm_results_list4 <- list()
for(q in unique(georef_all_dt_events_sf$SelfReference)){
  print(q)
  condition <- georef_all_dt_events_sf$SelfReference==q
  temp <- georef_all_dt_events_sf[condition,] %>% filter(!duplicated(event_hash))
  hex$dv <- sapply( st_covers(hex, temp), length) 
  lm1 <- MASS::glm.nb(dv ~ ., data= na.omit(as.data.frame(hex)[,vars_all]) )
  temp_results <- broom::tidy(lm1)
  temp_results$Type="SelfReference"
  temp_results$label <- q
  lm_results_list4[[as.character(q)]] <- temp_results
}

[1] TRUE

although coordinates are longitude/latitude, st_covers assumes that they are planar

[1] FALSE

although coordinates are longitude/latitude, st_covers assumes that they are planar
glm.fit: algorithm did not convergealternation limit reached

Plot the coefficients for each model

(Figure 5)

lm_results_dt <- data.table::rbindlist( list(  data.table::rbindlist(lm_results_list0),
                                   data.table::rbindlist(lm_results_list),
                                   data.table::rbindlist(lm_results_list2),
                                   data.table::rbindlist(lm_results_list3),
                                   data.table::rbindlist(lm_results_list4)
)
                           )
lm_results_dt <- lm_results_dt %>% mutate(rate= round(exp(estimate),2) )
sentence_case <- function(x) stringr::str_to_sentence(tolower(gsub("_"," ",x)))
#install.packages("extrafont");
#library(extrafont)
#library(extrafont)
#extrafont::font_import(prompt=F )
#capabilities()
#windowsFonts()
#sort(as.vector(unlist(windowsFonts())))
fonts <- c('Times New Roman',
           'Calibri',
           'Courier New',
           "Georgia",
           "Tunga")
           #'serif','Helvetica','Bookman','Palatino')
#qplot(1:10)+theme(text=element_text(family="Gill Sans Ultra Bold"))
#qplot(1:10)+theme(text=element_text(family=fonts[1]))
#qplot(1:10)+theme(text=element_text(family=fonts[2]))
#qplot(1:10)+theme(text=element_text(family=fonts[3]))
#qplot(1:10)+theme(text=element_text(family=fonts[4]))
#qplot(1:10)+theme(text=element_text(family=fonts[5]))
lm_results_dt <- data.table::as.data.table(lm_results_dt)
lm_results_dt[Type=="Match Type" & label=="FALSE", label:= "Exact",]
lm_results_dt[Type=="Match Type" & label=="TRUE", label:= "Fuzzy",]
lm_results_dt[Type=="SelfReference" & label=="FALSE", label:= "Self Ref.",]
lm_results_dt[Type=="SelfReference" & label=="TRUE", label:= "No Self Ref.",]
lm_results_dt[term=="treecover" , term:="tree cover",]
fontfaces <- factor(c("plain","bold","italic","bold.italic","plain"))
#p_load(ggrepel, tools)
p_lm <- lm_results_dt %>% 
        filter(!(label %in% c('kenya_district1962','kenya_cadastral','kenya_cadastral_district',"LINESTRING"))) %>%
        filter(term != "(Intercept)") %>%
        #mutate(label[Type=="Match Type"]=gsub("FALSE", "Exact", label[Type=="Match Type"])) %>%
        #mutate(label[Type=="Match Type"]=gsub("True", "Fuzzy", label[Type=="Match Type"])) %>%
        
        mutate(Type=sentence_case(Type),
               term=sentence_case(term),
               label=sentence_case(label)) %>%
  ggplot(aes(x=rate,
             y=p.value,
             color=Type,
             label=label,
             family = fonts[as.numeric(as.factor(Type))],
             fontface= fontfaces[as.numeric(as.factor(Type))]
             ))  + 
  facet_wrap(~term, scales="free") + 
  theme_bw() +
  geom_vline(xintercept = 1, , color="grey") + 
  geom_hline(yintercept = 0.01, linetype=3, color="grey") + 
  ggrepel::geom_text_repel(size=2) + 
  xlab(sentence_case("Coefficient Estimates (Change in relative rate of Violent Events)")) +
  ylab(sentence_case("Coefficient P-Value Estimates (0.01 Threshold Indicated)")) +
  theme(
    legend.position = c(0.9, 0.1), # c(0,0) bottom left, c(1,1) top-right.
  )
p_lm

 #+ coord_cartesian(y="log")
ggsave(
  filename = glue::glue(dir_figures, "p_lm.pdf"),
  plot = p_lm, width = 9, height = 6,
  device = cairo_pdf #have to use cairo to correctly embed the fonts
)

Interpretation

We evaluated five coefficients in a typical linear model predicting the number of violent events to occur in geographic area. The only thing we varied from model to model was the method by which we georeferenced the underlying events. We varied them across four main decisions, fuzzy/exact matching, geometry type, whether to allow self reference, and the source gazeteer. We further included two proposed cominations of decisions, one based on fixed hand tailored rules of when to prefer on decision to another, and another based on a supervsied learned model whose rules can vary based on specific context.

We found variation in our coefficient estimates across the different geographic decisions. The coefficeints varried in sign, magnitude, and statistical significance. We found variation in the stability of coefficients by each covariate. Both population and rainfall had stable estimates, with consistent sign, significance, and substantive effects with a few outliars. Distance, ruggedness, and treecover had unstable coefficients, with variation in significance and occasionaly the direction of the effect. In a few outliars, we observed a statistically significant effect in the opposite direction.