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

Lecture 19 Spatial GLM + Point Reference Spatial Data Colin Rundel 11/09/2017 1

Spatial GLM Models

Scottish Lip Cancer Data 2 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

3 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 4

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 5

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: 6

GLM Fit 7 Observed Cases GLM Predicted Cases 60 ° N 60 ° N 59 ° N 59 ° N 50 58 ° N 58 ° N 30 40 30 20 57 ° N 57 ° N 20 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

8 GLM Fit 50 40 glm_pred 30 20 10 0 0 10 20 30 40 Observed

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 57 ° N 57 ° N 20 −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

16 CAR Predictions 30 car_pred 20 10 0 0 10 20 30 40 Observed

CAR Residuals 17 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 3 58 ° N 58 ° N 2 30 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 () 18

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) }” 19

Model Parameters 20 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 21 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 22 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 40 58 ° N 58 ° N 5 30 0 20 57 ° N 57 ° N −5 −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 () 23

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)) { 24

IAR(2) Parameters 25 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 26 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

Predictions (cont.) 27 40 30 iar_pred 20 10 0 0 10 20 30 40 Observed

Residuals 28 Predicted Cases Residuals 60 ° N 60 ° N 59 ° N 59 ° N 40 58 ° N 58 ° N 5 30 0 20 57 ° N 57 ° N −5 −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 () 29

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 30

Point Referenced Data

Example - PM2.5 from CSN The Chemical Speciation Network are a series of air quality monitors run by 31 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