Lecture 20 Point referenced data (pt. 2) Colin Rundel 11/14/2018 - PowerPoint PPT Presentation

Lecture 20 Point referenced data (pt. 2) Colin Rundel 11/14/2018 1

Loa Loa Example

Loa Loa 2

Data 0.8 8 ## 0.78 0.502 503 11 163 4.88 11.4 16 7 7 ## 0.436 217 909 3 66 5.36 8.93 116 6 6 ## 0.85 0.415 109 3 8 8.07 5.10 0.8 stdev9901 <dbl> ## # ## # ... with 187 more rows, and 2 more variables: min9901 <dbl>, 0.84 0.487 268 4 57 6.00 9.31 104 10 ## 10 0.481 5.90 751 4 30 5.59 9.02 112 9 9 ## 0.69 0.373 103 0 83 167 8.18 loaloa = tbl_df (PrevMap::loaloa) %>% setNames (., tolower ( names (.))) %>% ## 2 2 ## 0.69 0.439 108 0 162 5.74 8.04 214 1 1 <dbl> 8.00 <dbl> <int> <int> <int> <dbl> <dbl> <int> <int> ## elev mean9901 max9901 row villcode longitude latitude no_exam no_inf ## ## # A tibble: 197 x 11 loaloa rename (elev=elevation) 215 5.68 212 ## 5 5 ## 0.67 0.432 104 5 62 5.92 8.10 219 4 4 0.79 167 0.491 783 5 88 5.35 8.91 118 3 3 ## 0.74 0.426 99 1 3

Spatial Distribution 4 12 ° N no_inf/no_exam 0.5 10 ° N 0.4 0.3 0.2 8 ° N 0.1 latitude 0.0 6 ° N no_exam 100 200 4 ° N 300 400 2 ° N 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude

Normalized Difference Vegetation Index (NDVI) 5 12˚N 11˚N 10˚N 9˚N 8˚N Latitude 7˚N 6˚N 5˚N 4˚N 3˚N 2˚N 8˚E 10˚E 12˚E 14˚E 16˚E Longitude -0.2 0 0.2 0.4 0.6 0.8 1 USGS LandDAAC MODIS version_005 WAF NDVI

Paper / Data summary Loa loa risk: decision making under uncertainty . Annals of Tropical Medicine and Parasitology, 101, 499-509. • no_exam and no_inf - Collected between 1991 and 2001 by NGOs (original paper mentions 168 villages and 21,938 observations) • elev - USGS gtopo30 (1km resolution) • mean9901 to stdev9901 - aggregated data from 1999 to 2001 from the Flemish Institute for Technological Research (1 km resolution) 6 Original paper - Diggle, et. al. (2007). Spatial modelling and prediction of

Diggle’s Model log ( 𝑞(𝑡) 1 − 𝑞(𝑡)) = 𝛽 + 𝑔 1 ( elev (𝑡)) + 𝑔 2 ( MAX.NDVI (𝑡)) + 𝑔 3 ( SD.NDVI (𝑡)) + 𝑥(𝑡) where 𝑥(𝑡) ∼ 𝒪(0, Σ) 7 {Σ} 𝑗𝑘 = 𝜏 2 exp (−𝑒 𝜚)

EDA 8 0 logit_prop −1 −2 −3 −4 −5 0 500 1000 1500 elev 0 logit_prop −1 −2 −3 −4 −5 0.7 0.8 0.9 max9901 0 logit_prop −1 −2 −3 −4 −5 0.12 0.15 0.18 0.21 stdev9901

Diggle’s EDA 9 t of

Data Augmentation 0.78 (0,0.8] 5 109 (0,1000] 0.85 (0.8,1] ## 6 909 (0,1000] 0.8 (0,0.8] ## 7 503 (0,1000] ## 0.67 (0,0.8] 8 103 (0,1000] 0.69 (0,0.8] ## 9 751 (0,1000] 0.8 (0,0.8] ## 10 268 (0,1000] 0.84 (0.8,1] ## # ... with 187 more rows ## 104 (0,1000] loaloa = loaloa %>% <dbl> <fct> mutate ( elev_f = cut (elev, breaks= c (0,1000,1300,2000), dig.lab=5), max_f = cut (max9901, breaks= c (0,0.8,1)) ) ## # A tibble: 197 x 4 ## elev elev_f max9901 max_f ## <int> <fct> ## 4 1 108 (0,1000] 0.69 (0,0.8] ## 2 99 (0,1000] 0.74 (0,0.8] ## 3 783 (0,1000] 0.79 (0,0.8] ## 10 loaloa %>% select (elev, elev_f, max9901, max_f)

Model Matrix {.} ## 109 0 0 ## 6 909 0 0 ## 7 503 0 0 8 ## 103 0 0 ## 9 751 0 0 ## 10 268 0 0 ## # ... with 187 more rows 5 0 model.matrix ( 0 ~ elev:elev_f - 1, data = loaloa as_data_frame () ## # A tibble: 197 x 3 ## `elev:elev_f(0,1000]` `elev:elev_f(1000,1300]` `elev:elev_f(1300,2000]` ## <dbl> <dbl> <dbl> ## 1 108 0 0 ## 2 99 0 0 ## 3 783 0 0 ## 4 104 11 ) %>%

OOS Validation loaloa_test = loaloa %>% sample_frac (0.20) 12 loaloa = anti_join (loaloa, loaloa_test, quiet=TRUE) Training 7 ° N no_inf/no_exam 6 ° N 0.5 latitude 0.4 5 ° N 0.3 4 ° N 0.2 3 ° N 0.1 2 ° N 0.0 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude Testing 7 ° N no_inf/no_exam 6 ° N 0.4 latitude 5 ° N 0.3 4 ° N 0.2 3 ° N 0.1 2 ° N 0.0 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude

Model -6.750 1.48e-11 *** 5.2632126 ## max9901:max_f(0.8,1] 7.617 2.60e-14 *** 0.6918702 5.2697375 ## max9901:max_f(0,0.8] 0.0002513 8.273 ## elev:elev_f(1300,2000] -0.0016964 3.507 0.000453 *** 0.0000953 0.0003343 ## elev:elev_f(1000,1300] < 2e-16 *** 0.6362108 < 2e-16 *** 0.0001018 degrees of freedom ## Number of Fisher Scoring iterations: 5 ## ## AIC: 2181 degrees of freedom on 151 ## Residual deviance: 1557.4 on 157 ## --- Null deviance: 3210.2 ## ## ## (Dispersion parameter for binomial family taken to be 1) ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Signif. codes: 15.660 0.0015951 g = glm (no_inf/no_exam ~ elev:elev_f + max9901:max_f + stdev9901, ## 3Q Median 1Q Min ## ## Deviance Residuals: stdev9901, family = binomial, data = loaloa, weights = loaloa$no_exam) ## -6.9522 ## ## glm(formula = no_inf/no_exam ~ elev:elev_f + max9901:max_f + ## Call: ## summary (g) data=loaloa, family=binomial, weights=loaloa$no_exam) Max -2.5662 ## elev:elev_f(0,1000] 0.5333413 -16.075 < 2e-16 *** 9.116 1.3070028 11.9141737 ## stdev9901 < 2e-16 *** -8.5735389 -0.4621 ## (Intercept) Estimate Std. Error z value Pr(>|z|) ## ## Coefficients: ## 10.1809 1.6720 13

Predictions - Training 14 Data 7 ° N no_inf/no_exam 6 ° N 0.5 latitude 0.4 5 ° N 0.3 4 ° N 0.2 3 ° N 0.1 2 ° N 0.0 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude GLM Prediction 7 ° N glm_pred 6 ° N 0.4 latitude 5 ° N 0.3 4 ° N 0.2 3 ° N 0.1 2 ° N 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude

Predictions - Testing 15 Data 7 ° N no_inf/no_exam 6 ° N 0.4 latitude 5 ° N 0.3 4 ° N 0.2 3 ° N 0.1 2 ° N 0.0 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude GLM Prediction 7 ° N glm_pred 6 ° N latitude 0.3 5 ° N 0.2 4 ° N 3 ° N 0.1 2 ° N 8 ° E 10 ° E 12 ° E 14 ° E 16 ° E longitude

Fit - Training 16 0.4 0.3 glm_pred 0.2 0.1 0.0 0.0 0.2 0.4 no_inf/no_exam

Fit - Testing 17 0.4 0.3 glm_pred 0.2 0.1 0.0 0.0 0.1 0.2 0.3 0.4 no_inf/no_exam

Spatial Structure? ## variog: computing omnidirectional variogram geoR:: variog (coords = cbind (loaloa$longitude, loaloa$latitude), 18 uvec = seq (0, 4, length.out = 50)) %>% plot () data = loaloa$prop - loaloa$glm_pred, 0.015 semivariance 0.010 0.005 0.000 0 1 2 3 4 distance

spBayes GLM Model spg = spBayes:: spGLM ( no_inf/no_exam ~ elev:elev_f + max9901:max_f + stdev9901, data=loaloa, family=”binomial”, weights=loaloa$no_exam, coords= cbind (loaloa$longitude, loaloa$latitude), cov.model=”exponential”, n.samples=20000, starting= list (beta= rep (0,7), phi=3, sigma.sq=1, w=0), priors= list (phi.unif= c (0.1, 10), sigma.sq.ig= c (2, 2)), amcmc= list (n.batch=1000, batch.length=20, accept.rate=0.43)) save (spg, loaloa, file=”loaloa.Rdata”) 19

spg$p.beta.theta.samples %>% -0.25813 -0.01064 -0.04942 -0.00030 max9901:max_f(0,0.8] 0.08517 -0.78200 -6.96111 9.06059 max9901:max_f(0.8,1] 0.69926 -5.79400 elev:elev_f(1300,2000] 9.08833 sigma.sq 0.45277 0.39071 0.14322 1.17856 phi 2.12385 1.44856 0.12026 8.46872 -0.01448 0.00176 post_summary () %>% stdev9901 knitr:: kable (digits=5) param post_mean post_med post_lower post_upper (Intercept) -7.62467 -7.10607 -15.33201 -1.56786 1.77896 -0.00471 -0.26705 -19.15846 24.59887 elev:elev_f(0,1000] 0.00010 0.00065 -0.00780 0.00316 elev:elev_f(1000,1300] -0.00059 -0.00035 20

Prediction 21 0.002 pred_spg_mean 0.001 0.000 0.0 0.2 0.4 prop

spBayes GLM Model - Fixed? spg_fix = spBayes:: spGLM ( data=loaloa, family=”binomial”, weights=loaloa$no_exam, coords= cbind (loaloa$longitude, loaloa$latitude), cov.model=”exponential”, n.samples=20000, starting= list (beta= rep (0,7), phi=3, sigma.sq=1, w=0), priors= list (phi.unif= c (0.1, 10), sigma.sq.ig= c (2, 2)), amcmc= list (n.batch=1000, batch.length=20, accept.rate=0.43) ) save (spg_fix, loaloa, file=”loaloa_fix.Rdata”) 22 no_inf ~ elev:elev_f + max9901:max_f + stdev9901,

param -0.18829 -0.00310 -0.00131 max9901:max_f(0,0.8] 0.74129 0.55728 -0.98971 2.78417 max9901:max_f(0.8,1] 1.15469 0.92740 2.89406 -0.00209 sigma.sq 1.26052 1.21204 0.32891 2.36502 phi 2.51439 2.38441 1.08064 4.86766 -0.00206 elev:elev_f(1300,2000] post_mean 1.02957 post_med post_lower post_upper (Intercept) -3.14223 -3.43877 -4.38140 -1.01108 stdev9901 1.88811 -5.28818 0.00039 9.04674 elev:elev_f(0,1000] 0.00036 0.00048 -0.00069 0.00114 elev:elev_f(1000,1300] -0.00036 -0.00031 -0.00127 23

Fit - Training 24 0.4 pred_spg_mean 0.2 0.0 0.0 0.2 0.4 no_inf/no_exam

25 Fit - Testing 0.3 pred_spg_mean 0.2 0.1 0.0 0.0 0.1 0.2 0.3 0.4 no_inf/no_exam

Diggle’s Predictive Surface 26 FIG. 2. Point estimates of the prevalence of Loa loa microfilaraemia, over-laid with the prevalences observed in field studies.

Exceedance Probability - Posterior Summary 27 Village 8 Village 127 60 40 20 village 0 Village 8 density Village 127 Village 104 Village 112 Village 104 60 Village 112 40 20 0 0.0 0.1 0.2 0.0 0.1 0.2 p