Lecture 19 Spatial GLM + Point Reference Spatial Data Colin Rundel - PowerPoint PPT Presentation

Lecture 19 Spatial GLM + Point Reference Spatial Data Colin Rundel 04/03/2017 1

Spatial GLM Models 2

Scottish Lip Cancer Data 3 Observed Expected 60 ° N 59 ° N value 80 58 ° N 60 40 57 ° N 20 0 56 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

4 Obs/Exp % Agg Fish Forest 60 ° N 60 ° N 59 ° N 59 ° N 6 20 58 ° N 58 ° N 4 15 10 57 ° N 57 ° N 2 5 0 0 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

Neighborhood / weight matrix 5

Moran’s I ## alternative hypothesis: greater Moran I test under randomisation ## ## data: lip_cancer$Observed/lip_cancer$Expected ## weights: listw ## ## Moran I statistic standard deviate = 8.2916, p-value < 2.2e-16 ## sample estimates: ## ## Moran I statistic Expectation Variance ## 0.589795225 -0.018181818 0.005376506 ## spdep:: moran.test (lip_cancer$Observed / lip_cancer$Expected, listw) spdep:: moran.test (lip_cancer$Observed, listw) ## ## ## Moran I test under randomisation ## ## data: lip_cancer$Observed ## weights: listw ## Moran I statistic standard deviate = 4.5416, p-value = 2.792e-06 0.005284831 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.311975396 -0.018181818 6

GLM ## (Dispersion parameter for poisson family taken to be 1) -7.80 6.21e-15 *** ## pcaff 0.073732 0.005956 12.38 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## ## (Intercept) -0.542268 ## Null deviance: 380.73 on 55 degrees of freedom ## Residual deviance: 238.62 on 54 degrees of freedom ## AIC: 450.6 ## ## Number of Fisher Scoring iterations: 5 0.069525 Estimate Std. Error z value Pr(>|z|) l = glm (Observed ~ offset ( log (Expected)) + pcaff, Min family=”poisson”, data=lip_cancer) summary (l) ## ## Call: ## glm(formula = Observed ~ offset(log(Expected)) + pcaff, family = ”poisson”, ## data = lip_cancer) ## ## Deviance Residuals: ## 1Q ## Median 3Q Max ## -4.7632 -1.2156 0.0967 1.3362 4.7130 ## ## Coefficients: 7

GLM Fit 8 Observed Cases GLM Predicted Cases 60 ° N 60 ° N 59 ° N 59 ° N 50 58 ° N 58 ° N 30 40 30 20 20 57 ° N 57 ° N 10 10 0 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

GLM Residuals 9 GLM Predicted Cases GLM Residuals 60 ° N 60 ° N 59 ° N 59 ° N 50 20 58 ° N 58 ° N 40 10 30 0 20 57 ° N 57 ° N −10 10 −20 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

Model Results ## 0.005323717 -0.018181818 0.333403223 ## Variance Expectation ## Moran I statistic ## sample estimates: ## alternative hypothesis: greater ## Moran I statistic standard deviate = 4.8186, p-value = 7.228e-07 ## weights: listw #RMSE lip_cancer$glm_resid ## data: ## Moran I test under randomisation ## ## spdep:: moran.test (lip_cancer$glm_resid, listw) #Moran's I ## [1] 7.480889 lip_cancer$glm_resid %>% .^2 %>% mean () %>% sqrt () 10

A hierachical model for lip cancer We have observed counts of lip cancer for 56 districts in Scotland. Let 𝑧 𝑗 represent the number of lip cancer for district 𝑗 . log (𝜇 𝑗 ) = log (𝐹 𝑗 ) + 𝑦 𝑗 𝛾 + 𝜕 𝑗 𝝏 ∼ 𝒪(𝟏, 𝜏 2 (𝐄 − 𝜚 𝐗) −1 ) 11 𝑧 𝑗 ∼ Poisson (𝜇 𝑗 ) where 𝐹 𝑗 is the expected counts for each region (and serves as an offet).

Data prep & CAR model for(i in 1:2) { }” phi ~ dunif(0,0.99) tau ~ dgamma(2, 2) sigma2 = 1/tau omega ~ dmnorm(rep(0,length(y)), tau * (D - phi*W)) } beta[i] ~ dnorm(0,1) } D = diag ( rowSums (W)) log(lambda[i]) = log_offset[i] + X[i,] %*% beta + omega[i] y_pred[i] ~ dpois(lambda[i]) y[i] ~ dpois(lambda[i]) for(i in 1:length(y)) { car_model = ”model{ y = lip_cancer$Observed 12 X = model.matrix (~ scale (lip_cancer$pcaff)) log_offset = log (lip_cancer$Expected)

CAR Results 13 0.5 beta[1] 0.0 −0.5 estimate 0.6 0.4 beta[2] 0.2 0.0 0 250 500 750 1000 .iteration

14 1.0 0.8 phi 0.6 estimate 0.4 2.0 1.5 sigma2 1.0 0.5 0 250 500 750 1000 .iteration

CAR Predictions 15 Observed Cases Predicted Cases 60 ° N 60 ° N 59 ° N 59 ° N 58 ° N 58 ° N 30 30 20 20 57 ° N 57 ° N 10 10 0 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

CAR Residuals 16 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 3 2 58 ° N 30 58 ° N 1 20 0 57 ° N 57 ° N −1 10 −2 −3 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

CAR Results ## 0.005658742 -0.018181818 0.061804482 ## Variance Expectation ## Moran I statistic ## sample estimates: ## alternative hypothesis: greater ## Moran I statistic standard deviate = 1.0633, p-value = 0.1438 ## weights: listw #RMSE lip_cancer$car_resid ## data: ## Moran I test under randomisation ## ## spdep:: moran.test (lip_cancer$car_resid, listw) #Moran's I ## [1] 1.586241 lip_cancer$car_resid %>% .^2 %>% mean () %>% sqrt () 17

Intrinsic Autoregressive Model iar_model = ”model{ for(i in 1:length(y)) { y[i] ~ dpois(lambda[i]) y_pred[i] ~ dpois(lambda[i]) log(lambda[i]) = log_offset[i] + X[i,] %*% beta + omega[i] } for(i in 1:2) { beta[i] ~ dnorm(0,1) } omega_free ~ dmnorm(rep(0,length(y)), tau * (D - W)) omega = omega_free - mean(omega_free) sigma2 = 1/tau tau ~ dgamma(2, 2) }” 18

Model Parameters 19 0.2 beta[1] 0.1 0.0 0.5 estimate 0.4 beta[2] 0.3 0.2 0.06055 0.06050 sigma2 0.06045 0.06040 0.06035 0 250 500 750 1000 .iteration

Predictions 20 Observed Cases Predicted Cases 60 ° N 60 ° N 59 ° N 59 ° N 40 58 ° N 58 ° N 30 30 20 20 57 ° N 57 ° N 10 10 0 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

Residuals 21 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 40 5 58 ° N 58 ° N 30 0 20 −5 57 ° N 57 ° N −10 10 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

IAR Results ## 0.005375965 -0.018181818 0.175216401 ## Variance Expectation ## Moran I statistic ## sample estimates: ## alternative hypothesis: greater ## Moran I statistic standard deviate = 2.6377, p-value = 0.004174 ## weights: listw #RMSE lip_cancer$iar_resid ## data: ## Moran I test under randomisation ## ## spdep:: moran.test (lip_cancer$iar_resid, listw) #Moran's I ## [1] 4.473069 lip_cancer$iar_resid %>% .^2 %>% mean () %>% sqrt () 22

Intrinsic Autoregressive Model - Reparameterized } }” tau ~ dgamma(2, 2) sigma = sqrt(sigma2) sigma2 = 1/tau omega = omega_free - mean(omega_free) omega_free ~ dmnorm(rep(0,length(y)), (D - W)) beta[i] ~ dnorm(0,1) iar_model2 = ”model{ for(i in 1:2) { } log(lambda[i]) = log_offset[i] + X[i,] %*% beta + sigma * omega[i] y_pred[i] ~ dpois(lambda[i]) y[i] ~ dpois(lambda[i]) for(i in 1:length(y)) { 23

IAR(2) Parameters 24 0.2 beta[1] 0.1 0.0 0.5 0.4 estimate 0.3 beta[2] 0.2 0.1 0.0 −0.1 1.6 1.2 sigma2 0.8 0.4 0 250 500 750 1000 .iteration

Predictions 25 Observed Cases Predicted Cases 60 ° N 60 ° N 59 ° N 59 ° N 40 58 ° N 58 ° N 30 30 20 20 57 ° N 57 ° N 10 10 0 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

Residuals 26 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 40 5 58 ° N 58 ° N 30 0 20 −5 57 ° N 57 ° N −10 10 56 ° N 56 ° N 55 ° N 55 ° N 8 ° W 6 ° W 4 ° W 2 ° W 8 ° W 6 ° W 4 ° W 2 ° W

IAR(2) Results ## 0.005617437 -0.018181818 0.011188088 ## Variance Expectation ## Moran I statistic ## sample estimates: ## alternative hypothesis: greater ## Moran I statistic standard deviate = 0.39186, p-value = 0.3476 ## weights: listw #RMSE lip_cancer$iar2_resid ## data: ## Moran I test under randomisation ## ## spdep:: moran.test (lip_cancer$iar2_resid, listw) #Moran's I ## [1] 1.654235 lip_cancer$iar2_resid %>% .^2 %>% mean () %>% sqrt () 27

Overall Results 0.0618 0.0112 1.6542 iar2 0.1752 4.4731 iar 1.5862 model car 0.3334 7.4809 glm moran rmse 28

Point Referenced Data 29

Example - PM2.5 from CSN The Chemical Speciation Network are a series of air quality monitors run by 30 2007 (n=191) for just PM2.5. EPA (221 locations in 2007). We’ll look at a subset of the data from Nov 11th, 50 45 pm25 40 40 lat 30 20 35 10 30 25 −120 −100 −80 long

csn ## 32.5 2007-11-14 00:00:00 22.6 ## 6 40139997 -112. 33.5 2007-11-14 00:00:00 12.3 ## 7 40191028 -111. 32.3 2007-11-14 00:00:00 7.20 8 51190007 5 11130001 -92.3 34.8 2007-11-14 00:00:00 12.7 ## 9 60070002 -122. 39.8 2007-11-14 00:00:00 10.0 ## 10 60190008 -120. 36.8 2007-11-14 00:00:00 32.3 ## # ... with 181 more rows -85.0 ## ## # A tibble: 191 x 5 -86.8 ## site longitude latitude date pm25 ## <int> <dbl> <dbl> <dttm> <dbl> ## 1 10730023 33.6 2007-11-14 00:00:00 19.4 32.4 2007-11-14 00:00:00 19.7 ## 2 10732003 -86.9 33.5 2007-11-14 00:00:00 26.4 ## 3 10890014 -86.6 34.7 2007-11-14 00:00:00 13.4 ## 4 11011002 -86.3 31

Aside - Splines 32