Lecture 18 Fitting CAR and SAR Models Colin Rundel 11/07/2018 1 - PowerPoint PPT Presentation

Lecture 18 Fitting CAR and SAR Models Colin Rundel 11/07/2018 1

Fitting areal models

Revised SAR Model • Formula Model 𝑧(𝑡 𝑗 ) = 𝑌 𝑗⋅ 𝛾 + 𝜚 𝑜 ∑ 𝑘=1 𝐸 −1 𝝑 ∼ 𝒪(𝟏, 𝜏 2 𝐄 −1 ) • Joint Model 𝐳 ∼ 𝒪 (𝐘𝜸, (𝐉 − 𝜚𝐄 −1 𝐁) −1 𝜏 2 𝐄 −1 ((𝐉 − 𝜚𝐄 −1 𝐁) −1 ) 𝑢 ) 2 𝑘𝑘 𝐵 𝑗𝑘 (𝑧(𝑡 𝑘 ) − 𝑌 𝑘⋅ 𝛾) + 𝜗 𝑗

Revised CAR Model • Conditional Model 𝑜 ∑ 𝑘=1 𝐵 𝑗𝑘 𝐸 𝑗𝑗 (𝑧(𝑡 𝑘 ) − 𝑌 𝑘⋅ 𝛾), 𝜏 2 𝐸 −1 • Joint Model 𝐳 ∼ 𝒪(𝐘𝜸, 𝜏 2 (𝐄 − 𝜚𝐁) −1 ) 3 𝑧(𝑡 𝑗 )|𝐳 −𝑡 𝑗 ∼ 𝒪 (𝑌 𝑗⋅ 𝛾 + 𝜚 𝑗𝑗 )

Example - NC SIDS 4 36.5 ° N BIR74 36 ° N 20000 35.5 ° N 15000 35 ° N 10000 34.5 ° N 5000 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 36.5 ° N SID74 36 ° N 40 35.5 ° N 30 35 ° N 20 34.5 ° N 10 0 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W

ggplot () + geom_sf (data=nc, aes (fill=SID74/BIR74*1000)) 5 36.5 ° N SID74/BIR74 * 1000 36 ° N 35.5 ° N 7.5 5.0 35 ° N 2.5 34.5 ° N 0.0 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W

Using spautolm from spdep ## Weights style: M ## M 100 10000 490 980 10696 S2 S1 S0 nn n ## ## Weights constants summary: ## library (spdep) ## Average number of links: 4.9 ## Percentage nonzero weights: 4.9 ## Number of nonzero links: 490 ## Number of regions: 100 ## Neighbour list object: ## Characteristics of weights list object: listW 6 A = st_touches (nc, sparse=FALSE) listW = mat2listw (A)

nc_coords = nc %>% st_centroid () %>% st_coordinates () plot ( st_geometry (nc)) plot (listW, nc_coords, add=TRUE, col=”blue”, pch=16) 7

Moran’s I ## alternative hypothesis: greater Moran I test under randomisation ## ## data: 1000 * nc$SID74/nc$BIR74 ## weights: listW ## ## Moran I statistic standard deviate = 3.6355, p-value = 0.0001387 ## sample estimates: ## ## Moran I statistic Expectation Variance ## 0.210046454 -0.010101010 0.003666802 ## spdep:: moran.test (1000*nc$SID74/nc$BIR74, listW) spdep:: moran.test (nc$SID74, listW) ## ## ## Moran I test under randomisation ## ## data: nc$SID74 ## weights: listW ## Moran I statistic standard deviate = 2.1707, p-value = 0.01498 0.003542176 ## alternative hypothesis: greater ## sample estimates: ## Moran I statistic Expectation Variance ## 0.119089049 -0.010101010 8

Geary’s C ## alternative hypothesis: Expectation greater than statistic Geary C test under randomisation ## ## data: nc$SID74/nc$BIR74 ## weights: listW ## ## Geary C statistic standard deviate = 3.0989, p-value = 0.0009711 ## sample estimates: ## ## Geary C statistic Expectation Variance ## 0.67796679 1.00000000 0.01079878 ## spdep:: geary.test (nc$SID74/nc$BIR74, listW) spdep:: geary.test (nc$SID74, listW) ## ## ## Geary C test under randomisation ## ## data: nc$SID74 ## weights: listW ## Geary C statistic standard deviate = 0.91949, p-value = 0.1789 0.01434105 ## alternative hypothesis: Expectation greater than statistic ## sample estimates: ## Geary C statistic Expectation Variance ## 0.88988684 1.00000000 9

CAR Model ## Lambda: 0.13062 LR test value: 8.6251 p-value: 0.0033157 ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 0.00200203 0.00024272 8.2484 2.22e-16 ## ## Numerical Hessian standard error of lambda: 0.030472 0.00768640 ## ## Log likelihood: 508.3767 ## ML residual variance (sigma squared): 2.1205e-06, (sigma: 0.0014562) ## Number of observations: 100 ## Number of parameters estimated: 3 ## AIC: -1010.8 ## 0.00055014 nc_car = spautolm (formula = SID74/BIR74 ~ 1, data = nc, ## listw = listW, family = ”CAR”) summary (nc_car) ## ## Call: spautolm(formula = SID74/BIR74 ~ 1, data = nc, listw = listW, ## family = ”CAR”) ## Residuals: ## -0.00213872 -0.00083535 -0.00022355 ## Min 1Q Median 3Q Max 10

SAR Model ## Lambda: 0.079934 LR test value: 8.8911 p-value: 0.0028657 ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 0.00201084 0.00023622 8.5127 < 2.2e-16 ## ## Numerical Hessian standard error of lambda: 0.0246 ## ## ## Log likelihood: 508.5097 ## ML residual variance (sigma squared): 2.1622e-06, (sigma: 0.0014704) ## Number of observations: 100 ## Number of parameters estimated: 3 ## AIC: -1011 ## Coefficients: 0.00762830 nc_sar = spautolm (formula = SID74/BIR74 ~ 1, data = nc, ## listw = listW, family = ”SAR”) summary (nc_sar) ## ## Call: spautolm(formula = SID74/BIR74 ~ 1, data = nc, listw = listW, ## family = ”SAR”) ## Residuals: 0.00051156 ## Min 1Q Median 3Q Max ## -0.00209307 -0.00087039 -0.00020274 11

Predictions 12 36.5 ° N car_pred 36 ° N 0.0035 35.5 ° N 0.0030 0.0025 35 ° N 0.0020 0.0015 34.5 ° N 0.0010 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 36.5 ° N sar_pred 36 ° N 0.0030 35.5 ° N 0.0025 35 ° N 0.0020 34.5 ° N 0.0015 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W

Residuals 13 36.5 ° N car_resid 36 ° N 0.0075 35.5 ° N 0.0050 35 ° N 0.0025 34.5 ° N 0.0000 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 36.5 ° N sar_resid 36 ° N 0.0075 35.5 ° N 0.0050 35 ° N 0.0025 34.5 ° N 0.0000 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W

Predicted vs Observed 14 0.003 car_pred 0.002 0.001 0.0000 0.0025 0.0050 0.0075 0.0100 SID74/BIR74 0.0030 0.0025 sar_pred 0.0020 0.0015 0.0000 0.0025 0.0050 0.0075 0.0100 SID74/BIR74

What’s wrong? 15 Histogram of nc$SID74 Histogram of nc$SID74/nc$BIR74 CAR Residuals 35 60 30 0.006 50 25 Sample Quantiles 40 Frequency Frequency 20 30 0.002 15 20 10 10 5 −0.002 0 0 0 10 20 30 40 0.000 0.004 0.008 −2 −1 0 1 2 nc$SID74 nc$SID74/nc$BIR74 Theoretical Quantiles

Comparing CAR vs SAR. 16 CAR vs SAR Residuals 0.006 nc$sar_resid 0.002 −0.002 −0.002 0.000 0.002 0.004 0.006 0.008 nc$car_resid

Stan CAR Model transformed parameters { ” } y ~ normal(beta+w_s, sigma2_w); sigma2_w ~ cauchy(0,5); sigma2 ~ cauchy(0,5); beta ~ normal(0,10); w_s ~ multi_normal_prec(rep_vector(0,N), Sigma_inv); matrix[N,N] Sigma_inv = (D - phi * A) / sigma2; model { } vector[N] y_pred = beta + w_s; } car_model = ” real<lower=0,upper=1> phi; real<lower=0> sigma2_w; real<lower=0> sigma2; real beta; vector[N] w_s; parameters { } matrix[N,N] D; matrix[N,N] A; vector[N] y; int<lower=0> N; data { 17

data = list ( N = nrow (nc), y = nc$SID74 / nc$BIR74, A = A * 1, D = diag ( rowSums (A)) ) car_fit = rstan:: stan ( model_code = car_model, data = data, iter = 10000, chains = 1, thin=20 ) Why don’t we use the conditional definition for the 𝑧 ’s? 18

data = list ( N = nrow (nc), y = nc$SID74 / nc$BIR74, A = A * 1, D = diag ( rowSums (A)) ) car_fit = rstan:: stan ( model_code = car_model, data = data, iter = 10000, chains = 1, thin=20 ) Why don’t we use the conditional definition for the 𝑧 ’s? 18

Model Results 19 0.0030 0.0025 beta 0.0020 0.0015 1.0 0.8 phi 0.6 estimate 0.4 sigma2 1e−05 5e−06 1e−03 sigma2_w 5e−04 0e+00 0 50 100 150 200 250 .iteration

Predictions 20 36.5 ° N bayes_car_pred 36 ° N 0.008 35.5 ° N 0.006 35 ° N 0.004 34.5 ° N 0.002 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 36.5 ° N bayes_car_resid 36 ° N 1e−03 35.5 ° N 5e−04 35 ° N 0e+00 34.5 ° N 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W

21 0.003 0.002 bayes_car_pred 0.001 0.000 0.000 0.001 0.002 0.003 car_pred

Brief Aside - SAR Precision Matrix 𝑢 Σ −1 𝑢 ) −1 𝑢 ) 22 Σ 𝑇𝐵𝑆 = (𝐉 − 𝜚𝐄 −1 𝐁) −1 𝜏 2 𝐄 −1 ((𝐉 − 𝜚𝐄 −1 𝐁) −1 ) 𝑇𝐵𝑆 = ((𝐉 − 𝜚𝐄 −1 𝐁) −1 𝜏 2 𝐄 −1 ((𝐉 − 𝜚𝐄 −1 𝐁) −1 ) −1 1 = (((𝐉 − 𝜚𝐄 −1 𝐁) −1 ) 𝜏 2 𝐄 (𝐉 − 𝜚𝐄 −1 𝐁) = 1 𝜏 2 (𝐉 − 𝜚𝐄 −1 𝐁) 𝑢 𝐄 (𝐉 − 𝜚𝐄 −1 𝐁)

Jags SAR Model N = nrow (nc), w_s ~ multi_normal_prec(rep_vector(0,N), Sigma_inv); beta ~ normal(0,10); sigma2 ~ cauchy(0,5); sigma2_w ~ cauchy(0,5); y ~ normal(beta + w_s, sigma2_w); } ” D = diag ( rowSums (A)) y = nc$SID74 / nc$BIR74, sar_model = ” x = rep (1, nrow (nc)), D_inv = D_inv, W_tilde = D_inv %*% A ) sar_fit = rstan:: stan ( model_code = sar_model, data = data, iter = 10000, chains = 1, thin=20 ) matrix[N,N] Sigma_inv = C' * D * C / sigma2; matrix[N,N] C = I - phi * W_tilde; model { parameters { data { int<lower=0> N; vector[N] y; matrix[N,N] W_tilde; matrix[N,N] D; } transformed data { matrix[N,N] I = diag_matrix(rep_vector(1, N)); } 23 } vector[N] w_s; real beta; real<lower=0> sigma2; real<lower=0> sigma2_w; real<lower=0,upper=1> phi; } transformed parameters { vector[N] y_pred = beta + w_s; D_inv = diag (1/ diag (D)) data = list (

Model Results 24 0.0032 0.0028 0.0024 beta 0.0020 0.0016 0.0012 0.75 phi 0.50 estimate 0.25 1e−05 sigma2 5e−06 0.0015 sigma2_w 0.0010 0.0005 0.0000 0 50 100 150 200 250 .iteration