Lecture 19 Fitting CAR and SAR Models Colin Rundel 03/29/2017 1 - PowerPoint PPT Presentation

Lecture 19 Fitting CAR and SAR Models Colin Rundel 03/29/2017 1

Fitting areal models 2

CAR vs SAR • Conditional Autoregressive (CAR) n • Simultaneous Autogressve (SAR) 3 n ∑ y ( s i ) = ϕ W ij y ( s j ) + ϵ j = 1 y ∼ N ( 0 , σ 2 (( I − ϕ W ) − 1 )(( I − ϕ W ) − 1 ) t )   ∑ y ( s i ) | y − s i ∼ N  ϕ W ij y ( s j ) , σ 2  j = 1 y ∼ N ( 0 , σ 2 ( I − ϕ W ) − 1 )

Some specific generalizations Generally speaking we will want to work with a scaled / normalized version where D is as defined above. 1 4 W ij When W is an adjacency matrix we can express this as of the weight matrix, W i · D − 1 W where D = diag ( m i ) and m i = | N ( s i ) | . We can also allow σ 2 to vary between locations, we can define this as D τ = diag ( 1 /σ 2 i ) and most often we use ( ) D τ = diag = D /σ 2 σ 2 / | N ( s i ) |

Revised CAR Model W ij • Joint Model ii • Conditional Model D ii 5 n   ∑ y ( s i ) | y − s i ∼ N  X i · β + ϕ ( y ( s j ) − X j · β ) , σ 2 D − 1  j = 1 y ∼ N ( X β , Σ CAR ) Σ CAR = ( D σ ( I − ϕ D − 1 W )) − 1 = ( 1 /σ 2 D ( I − ϕ D − 1 W )) − 1 = ( 1 /σ 2 ( D − ϕ W )) − 1 = σ 2 ( D − ϕ W ) − 1

Revised SAR Model jj • Formula Model • Joint Model W ij 6 n ) + ϵ i ∑ ( y ( s i ) = X i · β + ϕ y ( s j ) − X j · β D − 1 j = 1 ) + ϵ y = X β + ϕ D − 1 W ( y − X β ) = ϕ D − 1 W ) + ϵ ( ( y − X β y − X β ) ( I − ϕ D − 1 W ) − 1 = ϵ ( y − X β y = X β + ( I − ϕ D − 1 W ) − 1 ϵ ( X β , ( I − ϕ D − 1 W ) − 1 σ 2 D − 1 ( ( I − ϕ D − 1 W ) − 1 ) t ) y ∼ N

2 D Toy CAR Example 0 1 1 0 2 0 1 2 W 1 0 0 0 2 0 7 0 W 0 1 0 1 0 1 0 1 0 D 1 s2 s1 s3

Toy CAR Example 0 1 0 2 0 1 1 0 0 0 2 0 0 1 0 1 1 0 1 0 1 0 7 0 s2 s1 s3     W = D =           − 1 − ϕ Σ = σ 2 ( D − ϕ W ) = σ 2 − ϕ − ϕ     − ϕ

8 [,2] [,3] ## [1,] 1.1666667 0.3333333 0.1666667 ## [2,] 0.3333333 0.6666667 0.3333333 ## [3,] 0.1666667 0.3333333 1.1666667 check_sigma (phi=-0.6) ## [,1] [,3] [,1] ## [1,] 1.28125 -0.46875 0.28125 ## [2,] -0.46875 0.78125 -0.46875 ## [3,] 0.28125 -0.46875 1.28125 [,2] ## check_sigma = function (phi) { 0.0 solve (Sigma_inv) } check_sigma (phi=0) ## [,1] [,2] [,3] ## [1,] 1 0 check_sigma (phi=0.5) ## [2,] 0 0.5 0 ## [3,] 0 0.0 1 When does Σ exist? Sigma_inv = matrix ( c (1,-phi,0,-phi,2,-phi,0,-phi,1), ncol=3, byrow=TRUE)

check_sigma (phi=1) ## 1.363636 -0.6363636 ## [3,] -1.6363636 1.3636364 1.3636364 -1.136364 ## [2,] 1.363636 -1.6363636 ## [1,] -0.6363636 [,3] [,2] [,1] check_sigma (phi=-1.2) ## Error in solve.default(Sigma_inv): Lapack routine dgesv: system is exactly singular: U[3,3] = 0 ## [3,] -1.6363636 -1.363636 -0.6363636 ## [2,] -1.3636364 -1.136364 -1.3636364 ## [1,] -0.6363636 -1.363636 -1.6363636 [,3] [,2] [,1] ## check_sigma (phi=1.2) ## Error in solve.default(Sigma_inv): Lapack routine dgesv: system is exactly singular: U[3,3] = 0 check_sigma (phi=-1) 9

Conclusions Generally speaking just like the AR(1) model for time series we require that autoregressive (IAR) model and they are popular as an improper prior for spatial random effects. An additional sum constraint is necessary for 10 | ϕ | < 1 for the CAR model to be proper. These results for ϕ also apply in the context where σ 2 i is constant across locations (i.e. Σ = ( σ 2 ( I − ϕ D − 1 W )) − 1 ) As a side note, the special case where ϕ = 1 is known as an intrinsic identifiability ( ∑ i = 1 n y ( s i ) = 0).

Example - NC SIDS 11 36.5 ° N BIR79 36 ° N 30000 35.5 ° N 20000 35 ° N 10000 34.5 ° N 34 ° N 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 36.5 ° N SID79 36 ° N 50 35.5 ° N 40 30 35 ° N 20 10 34.5 ° N 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 12 W = st_touches (nc, sparse=FALSE) listW = mat2listw (W)

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

CAR Model ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.06911902 0.67501301 1.5838 0.1132 ## BIR74 0.00175249 0.00010107 17.3401 <2e-16 ## Lambda: 0.13222 LR test value: 8.8654 p-value: 0.0029062 ## Coefficients: ## Numerical Hessian standard error of lambda: 0.030094 ## ## Log likelihood: -275.7655 ## ML residual variance (sigma squared): 13.695, (sigma: 3.7007) ## Number of observations: 100 ## Number of parameters estimated: 4 ## AIC: 559.53 ## ## nc_car = spautolm (formula = SID74 ~ BIR74, data = nc, ## listw = listW, family = ”CAR”) summary (nc_car) ## ## Call: ## spautolm(formula = SID74 ~ BIR74, data = nc, listw = listW, family = ”CAR”) ## ## Residuals: Min 13.54059 1Q Median 3Q Max ## -10.38934 -1.58600 -0.52154 1.14729 14

SAR Model <2e-16 Estimate Std. Error z value Pr(>|z|) ## (Intercept) 1.01971585 0.64910408 1.571 0.1162 ## BIR74 0.00174741 0.00010105 17.292 ## ## Coefficients: ## Lambda: 0.075265 LR test value: 8.4013 p-value: 0.0037495 ## Numerical Hessian standard error of lambda: 0.024085 ## ## Log likelihood: -275.9975 ## ML residual variance (sigma squared): 14.158, (sigma: 3.7627) ## Number of observations: 100 ## Number of parameters estimated: 4 ## AIC: 560 ## ## nc_sar = spautolm (formula = SID74 ~ BIR74, data = nc, ## listw = listW, family = ”SAR”) summary (nc_sar) ## ## Call: ## spautolm(formula = SID74 ~ BIR74, data = nc, listw = listW, family = ”SAR”) ## ## Residuals: Min 14.70511 1Q Median 3Q Max ## -10.94771 -1.72354 -0.56866 1.23273 15

Residuals 16 car_pred sar_pred 36.5 ° N 36.5 ° N 36 ° N 36 ° N 30 30 35.5 ° N 35.5 ° N 35 ° N 20 35 ° N 20 34.5 ° N 34.5 ° N 34 ° N 10 34 ° N 10 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 0 car_resid car_resid 36.5 ° N 36.5 ° N 10 10 36 ° N 36 ° N 35.5 ° N 35.5 ° N 5 5 35 ° N 35 ° N 0 0 34.5 ° N 34.5 ° N 34 ° N 34 ° N −5 −5 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W −10 −10

I agree … 17

Why? 18 Histogram of nc$SID74 CAR Residuals CAR vs SAR Residuals 15 60 10 10 50 Sample Quantiles 40 nc$sar_resid Frequency 5 5 30 0 0 20 −5 −5 10 −10 −10 0 0 10 20 30 40 −2 −1 0 1 2 −10 −5 0 5 10 nc$SID74 Theoretical Quantiles nc$car_resid

Jags CAR Model ## model{ Why don’t we use the conditional definition for the y ’s? x = nc$BIR74 y = nc$SID74 ## } phi ~ dunif(-0.99, 0.99) ## sigma2 ~ dnorm(0, 1/100) T(0,) ## tau <- 1 / sigma2 ## ## beta1 ~ dnorm(0, 1/100) ## beta0 ~ dnorm(0, 1/100) ## ## y_pred ~ dmnorm(beta0 + beta1*x, tau * (D - phi*W)) ## y ~ dmnorm(beta0 + beta1*x, tau * (D - phi*W)) ## 19 W = W * 1L D = diag ( rowSums (W))

Jags CAR Model ## model{ Why don’t we use the conditional definition for the y ’s? x = nc$BIR74 y = nc$SID74 ## } phi ~ dunif(-0.99, 0.99) ## sigma2 ~ dnorm(0, 1/100) T(0,) ## tau <- 1 / sigma2 ## ## beta1 ~ dnorm(0, 1/100) ## beta0 ~ dnorm(0, 1/100) ## ## y_pred ~ dmnorm(beta0 + beta1*x, tau * (D - phi*W)) ## y ~ dmnorm(beta0 + beta1*x, tau * (D - phi*W)) ## 19 W = W * 1L D = diag ( rowSums (W))

Model Results 20 Trace of beta0 Density of beta0 0.4 6 0.2 2 −2 0.0 16000 18000 20000 22000 24000 −2 0 2 4 6 8 10 Iterations N = 1000 Bandwidth = 0.1893 Trace of beta1 Density of beta1 0.0017 2000 0.0013 0 16000 18000 20000 22000 24000 0.0014 0.0016 0.0018 0.0020 Iterations N = 1000 Bandwidth = 2.279e−05

21 Trace of phi Density of phi 1.0 2.0 0.6 1.0 0.2 0.0 16000 18000 20000 22000 24000 0.2 0.4 0.6 0.8 1.0 Iterations N = 1000 Bandwidth = 0.03673 Trace of sigma2 Density of sigma2 0.08 60 0.04 50 40 0.00 16000 18000 20000 22000 24000 40 50 60 70 Iterations N = 1000 Bandwidth = 1.188

Predictions ## [1] 3.72107 nc$jags_pred = y_pred$post_mean ## [1] 3.762664 sqrt ( mean (nc$sar_resid^2)) 22 sqrt ( mean (nc$car_resid^2)) ## [1] 3.987985 sqrt ( mean (nc$jags_resid^2)) nc$jags_resid = nc$SID74 - y_pred$post_mean jags_pred jags_resid 15 36.5 ° N 36.5 ° N 30 36 ° N 36 ° N 10 35.5 ° N 35.5 ° N 5 35 ° N 35 ° N 20 34.5 ° N 34.5 ° N 0 34 ° N 34 ° N 10 −5 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W 84 ° W 82 ° W 80 ° W 78 ° W 76 ° W −10

23 30 jags_pred 20 10 0 0 10 20 30 40 car_pred

1 Brief Aside - SAR Precision Matrix 24 Σ SAR = ( I − ϕ D − 1 W ) − 1 σ 2 D − 1 ( ( I − ϕ D − 1 W ) − 1 ) t ( ( I − ϕ D − 1 W ) − 1 σ 2 D − 1 ( ( I − ϕ D − 1 W ) − 1 ) t ) − 1 Σ − 1 SAR = (( ( I − ϕ D − 1 W ) − 1 ) t ) − 1 1 σ 2 D ( I − ϕ D − 1 W ) = = σ 2 ( I − ϕ D − 1 W ) t D ( I − ϕ D − 1 W )