[PPT] - Bayesian Spatial Analysis Dr. Jarad Niemi STAT 615 - Iowa State PowerPoint Presentation

SLIDE 1

Bayesian Spatial Analysis

Dr. Jarad Niemi

STAT 615 - Iowa State University

November 9, 2017

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 1 / 51

SLIDE 2

Spatial modeling

Three main types of spatial data: Point/geo-referenced Areal-referenced Point process/pattern

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 2 / 51

SLIDE 3

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 4

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 5

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Some quantity, Y (s), is measured at each location s ∈ D Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 6

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Some quantity, Y (s), is measured at each location s ∈ D Examples Air quality monitoring

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 7

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Some quantity, Y (s), is measured at each location s ∈ D Examples Air quality monitoring Coastal tide level monitoring

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 8

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Some quantity, Y (s), is measured at each location s ∈ D Examples Air quality monitoring Coastal tide level monitoring Earthquake monitoring

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 9

Spatial modeling

Point-referenced spatial data

Features Some spatial domain D is under study Measured spatial locations s ∈ D are pre-determined Some quantity, Y (s), is measured at each location s ∈ D Examples Air quality monitoring Coastal tide level monitoring Earthquake monitoring Bird point counts

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 3 / 51

SLIDE 10

Spatial modeling

Areal-referenced spatial data

Features Some set of spatial regions 1, . . . , S are pre-determined Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 4 / 51

SLIDE 11

Spatial modeling

Areal-referenced spatial data

Features Some set of spatial regions 1, . . . , S are pre-determined Some quantity, Ys, is measured as an aggregate over that region Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 4 / 51

SLIDE 12

Spatial modeling

Areal-referenced spatial data

Features Some set of spatial regions 1, . . . , S are pre-determined Some quantity, Ys, is measured as an aggregate over that region Examples Disease occurrence per county

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 4 / 51

SLIDE 13

Spatial modeling

Areal-referenced spatial data

Features Some set of spatial regions 1, . . . , S are pre-determined Some quantity, Ys, is measured as an aggregate over that region Examples Disease occurrence per county Unemployment rate per state

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 4 / 51

SLIDE 14

Spatial modeling

Areal-referenced spatial data

Features Some set of spatial regions 1, . . . , S are pre-determined Some quantity, Ys, is measured as an aggregate over that region Examples Disease occurrence per county Unemployment rate per state Inflation per country

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 4 / 51

SLIDE 15

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 16

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Spatial locations s ∈ S ⊂ D are random Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 17

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Spatial locations s ∈ S ⊂ D are random Y (s) = 1 indicates an occurence of the event Examples

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 18

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Spatial locations s ∈ S ⊂ D are random Y (s) = 1 indicates an occurence of the event Examples Locations of Mayan ruins

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 19

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Spatial locations s ∈ S ⊂ D are random Y (s) = 1 indicates an occurence of the event Examples Locations of Mayan ruins Locations of invasive species

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 20

Spatial modeling

Point-process spatial data

Features Some spatial domain D is under study Spatial locations s ∈ S ⊂ D are random Y (s) = 1 indicates an occurence of the event Examples Locations of Mayan ruins Locations of invasive species Locations of caught Lingcod

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 5 / 51

SLIDE 21

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 22

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 23

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 24

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Intrinsic stationarity

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 25

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Intrinsic stationarity Weak stationarity

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 26

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Intrinsic stationarity Weak stationarity Strong stationarity

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 27

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Intrinsic stationarity Weak stationarity Strong stationarity

Isotropy

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 28

Point-referenced spatial data

Let Y (s) for s ∈ D ⊆ Rd be a spatial process. Let E[Y (s)] = 0 for all s ∈ D because we will model the mean separately. Assumptions: Stationarity

Intrinsic stationarity Weak stationarity Strong stationarity

Isotropy Gaussian process

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 6 / 51

SLIDE 29

Point-referenced spatial data

Example spatial process

3 4 5 6 7 400 600 800 1000 1400 MgCl2 dNTP 7 8 9 10 11

8 9 10 11

log of DNA amplification rate (KCL=29.77, KPO4=32.13)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 7 / 51

SLIDE 30

Point-referenced spatial data Assumptions

Intrinsic stationarity

Definition A process Y (s) is intrinsically stationary if (E[Y (s + h) − Y (s)] = 0 and) E[(Y (s + h) − Y (s))2] = V ar[Y (s + h) − Y (s)] = 2γ(h) when s, s + h ∈ D.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 8 / 51

SLIDE 31

Point-referenced spatial data Assumptions

Intrinsic stationarity

Definition A process Y (s) is intrinsically stationary if (E[Y (s + h) − Y (s)] = 0 and) E[(Y (s + h) − Y (s))2] = V ar[Y (s + h) − Y (s)] = 2γ(h) when s, s + h ∈ D. We call 2γ(h) the variogram and γ(h) the semivariogram.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 8 / 51

SLIDE 32

Point-referenced spatial data Assumptions

Intrinsic stationarity

Definition A process Y (s) is intrinsically stationary if (E[Y (s + h) − Y (s)] = 0 and) E[(Y (s + h) − Y (s))2] = V ar[Y (s + h) − Y (s)] = 2γ(h) when s, s + h ∈ D. We call 2γ(h) the variogram and γ(h) the semivariogram. Definition A process Y (s) is isotropic if the semivariogram function depends only on ||h||, the length of the separation vector. Otherwise the process is anisotropic.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 8 / 51

SLIDE 33

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 34

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 35

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram. Since γ(h) = C(0) − C(h), a weakly stationary process is also intrinsicly stationary.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 36

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram. Since γ(h) = C(0) − C(h), a weakly stationary process is also intrinsicly stationary. If the spatial process is ergodic, then C(h) → 0 as ||h|| → ∞

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 37

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram. Since γ(h) = C(0) − C(h), a weakly stationary process is also intrinsicly stationary. If the spatial process is ergodic, then C(h) → 0 as ||h|| → ∞ and lim||h||→∞ γ(h) = C(0).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 38

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram. Since γ(h) = C(0) − C(h), a weakly stationary process is also intrinsicly stationary. If the spatial process is ergodic, then C(h) → 0 as ||h|| → ∞ and lim||h||→∞ γ(h) = C(0). Thus C(h) = C(0) − γ(h) = lim

||u||→∞ γ(u) − γ(h).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 39

Point-referenced spatial data Assumptions

Weak stationarity

Definition A process Y (s) has weak stationarity if (E[Y (s)] = µ and) Cov[Y (s), Y (s + h)] = C(h) when s, s + h ∈ D. We call C(h) the covariance function or covariogram. Since γ(h) = C(0) − C(h), a weakly stationary process is also intrinsicly stationary. If the spatial process is ergodic, then C(h) → 0 as ||h|| → ∞ and lim||h||→∞ γ(h) = C(0). Thus C(h) = C(0) − γ(h) = lim

||u||→∞ γ(u) − γ(h).

Thus, if the process is ergodic, an intrinsicly stationary process is also weakly stationary.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51

SLIDE 40

Point-referenced spatial data Assumptions

Covariance functions for isotropic models

Model Covariance function, C(t) Semivariogram, γ(t) Linear C(t) does not exist γ(t) =

τ2 + σ2t

if t > 0

therwise

Spherical C(t) =      σ2 1 − 3

2 φt + 1 2 (φt)3

τ2 + σ2 γ(t) =      τ2 + σ2 if τ2 + σ2

3 2 φt − 1 2 (φt)3

therwise

Exponential C(t) =

σ2 exp(−φt)

τ2 + σ2 γ(t) =

τ2 + σ2 [1 − exp(−φt)]

t >

therwise

Powered exponential C(t) =

σ2 exp(−|φt|p)

τ2 + σ2 γ(t) =

τ2 + σ2 [1 − exp(−|φt|p)]

t

therwise

Gaussian C(t) =

σ2 exp(−φ2t2)

τ2 + σ2 γ(t) =

τ2 + σ2

1 − exp(−φ2t2)

Rational quadratic

C(t) =    σ2

1 −

t2 (1+φ2)

τ2 + σ2

γ(t) =

τ2 + σ2

t2 (1+φ2)

t > 0

therwise

Wave C(t) =

σ2 sin(φt)

φt

τ2 + σ2 γ(t) =

τ2 + σ2

1 − sin(φt)

φt

t > 0
therwise

Power law C(t) does not exist γ(t) =

τ2 + σ2tλ

t > 0

therwise

Mat´ ern C(t) =   

σ2 2ν−1Γ(ν) (2√vtφ)νKν(2√νtφ)

τ2 + σ2 γ(t) =    τ2 + σ2

1 − (2√νtφ)ν

2ν−1Γ(ν) (2√v

Mat´ ern (ν = 3/2) C(t) =

σ2(1 + φt) exp(−φt)

τ2 + σ2 γ(t) =

τ2 + σ2 [1 − (1 + φt) exp(−φt

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 10 / 51

SLIDE 41

Point-referenced spatial data Assumptions

Spherical semivariogram

0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2

Spherical semivariogram

t γ(t) nugget (τ2) partial sill (σ2) sill (τ2 + σ2) range (1 φ)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 11 / 51

SLIDE 42

Point-referenced spatial data Assumptions

Mat´ ern

Perhaps the most important isotropic process is the Mat´ ern process with covariance C(t) =

σ2

2ν−1Γ(ν)(2√νtφ)νKν(2√νtφ)

t > 0 τ 2 + σ2 t = 0 and variogram γ(t) =

τ 2 + σ2

1 − (2√νtφ)ν

2ν−1Γ(ν)Kν(2√νtφ)

t > 0
therwise

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 12 / 51

SLIDE 43

Point-referenced spatial data Assumptions

Mat´ ern

Perhaps the most important isotropic process is the Mat´ ern process with covariance C(t) =

σ2

2ν−1Γ(ν)(2√νtφ)νKν(2√νtφ)

t > 0 τ 2 + σ2 t = 0 and variogram γ(t) =

τ 2 + σ2

1 − (2√νtφ)ν

2ν−1Γ(ν)Kν(2√νtφ)

t > 0
therwise

where ν controls the smoothness of the spatial process (⌊ν⌋ number of times process realizations are mean square differentiable) while φ is a spatial scale parameter.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 12 / 51

SLIDE 44

Point-referenced spatial data Assumptions

Mat´ ern

Perhaps the most important isotropic process is the Mat´ ern process with covariance C(t) =

σ2

2ν−1Γ(ν)(2√νtφ)νKν(2√νtφ)

t > 0 τ 2 + σ2 t = 0 and variogram γ(t) =

τ 2 + σ2

1 − (2√νtφ)ν

2ν−1Γ(ν)Kν(2√νtφ)

t > 0
therwise

where ν controls the smoothness of the spatial process (⌊ν⌋ number of times process realizations are mean square differentiable) while φ is a spatial scale parameter. Special cases are the exponential (ν = 1/2) and Gaussian (ν → ∞).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 12 / 51

SLIDE 45

Point-referenced spatial data Assumptions

Strong stationarity

Definition A process Y (s) is strongly (or strictly) stationary if, for any set of n ≥ 1 sites {s1, . . . , sn} and any h ∈ Rd, (Y (s1), . . . , Y (sn))⊤ d = (Y (s1 + h), . . . , Y (sn + h))⊤ where d = means equal in distribution.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 13 / 51

SLIDE 46

Point-referenced spatial data Assumptions

Strong stationarity

Definition A process Y (s) is strongly (or strictly) stationary if, for any set of n ≥ 1 sites {s1, . . . , sn} and any h ∈ Rd, (Y (s1), . . . , Y (sn))⊤ d = (Y (s1 + h), . . . , Y (sn + h))⊤ where d = means equal in distribution. If we assume all variances exist, then strong stationarity implies weak stationarity.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 13 / 51

SLIDE 47

Point-referenced spatial data Assumptions

Strong stationarity

Definition A process Y (s) is strongly (or strictly) stationary if, for any set of n ≥ 1 sites {s1, . . . , sn} and any h ∈ Rd, (Y (s1), . . . , Y (sn))⊤ d = (Y (s1 + h), . . . , Y (sn + h))⊤ where d = means equal in distribution. If we assume all variances exist, then strong stationarity implies weak stationarity. The reverse is not necessarily true.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 13 / 51

SLIDE 48

Point-referenced spatial data Gaussian process

Gaussian process

Definition Y (s) is a Gaussian process if, for any n ≥ 1 and any set of sites {s1, . . . , sn}, Y = (Y (s1), . . . , Y (sn))⊤ has a multivariate normal distribution.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 14 / 51

SLIDE 49

Point-referenced spatial data Gaussian process

Gaussian process

Definition Y (s) is a Gaussian process if, for any n ≥ 1 and any set of sites {s1, . . . , sn}, Y = (Y (s1), . . . , Y (sn))⊤ has a multivariate normal distribution. For a Gaussian process, weak stationarity and strong stationarity are equivalent.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 14 / 51

SLIDE 50

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 51

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 52

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 53

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process and according to a particular covariance function

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 54

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process and according to a particular covariance function. Collectively refer to the parameters as θ, then our objective is p(θ|y) ∝ p(y|θ)p(θ).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 55

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process and according to a particular covariance function. Collectively refer to the parameters as θ, then our objective is p(θ|y) ∝ p(y|θ)p(θ). Suppose we assume the Mat´ ern covariance function and a common mean µ so that θ = (µ, ν, φ, τ 2, σ2).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 56

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process and according to a particular covariance function. Collectively refer to the parameters as θ, then our objective is p(θ|y) ∝ p(y|θ)p(θ). Suppose we assume the Mat´ ern covariance function and a common mean µ so that θ = (µ, ν, φ, τ 2, σ2). Then we have p(µ, ν, φ, τ 2, σ2|y) ∝ N(y; µ, Σ)p(µ, ν, φ, τ 2, σ2)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 57

Point-referenced spatial data Estimation

Bayesian estimation of Gaussian process parameters

Suppose we observe data at some locations s1, . . . , sn. Collectively, we have y = (y(s1), . . . , y(sn)). Let’s assume the data arise from a Gaussian process and according to a particular covariance function. Collectively refer to the parameters as θ, then our objective is p(θ|y) ∝ p(y|θ)p(θ). Suppose we assume the Mat´ ern covariance function and a common mean µ so that θ = (µ, ν, φ, τ 2, σ2). Then we have p(µ, ν, φ, τ 2, σ2|y) ∝ N(y; µ, Σ)p(µ, ν, φ, τ 2, σ2) where Σ is constructed from the parameters ν, φ, τ 2, and σ2 and the distances between locations, e.g. ||s1 − s2||.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51

SLIDE 58

Point-referenced spatial data Estimation

Consider point-referenced data at spatial locations s1, . . . , sn, model this data as Y (s) = µ(s) + w(s) + ǫ(s) If we constrain ourselves to isotropic models, the Mat´ ern class is suggested as a general tool (Banerjee pg. 37).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51

SLIDE 59

Point-referenced spatial data Estimation

Consider point-referenced data at spatial locations s1, . . . , sn, model this data as Y (s) = µ(s) + w(s) + ǫ(s) If we constrain ourselves to isotropic models, the Mat´ ern class is suggested as a general tool (Banerjee pg. 37). If w = (w(s1), . . . , w(sn))⊤ and ǫ = (ǫ(s1), . . . , ǫ(sn))⊤, then a general model is

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51

SLIDE 60

Point-referenced spatial data Estimation

Consider point-referenced data at spatial locations s1, . . . , sn, model this data as Y (s) = µ(s) + w(s) + ǫ(s) If we constrain ourselves to isotropic models, the Mat´ ern class is suggested as a general tool (Banerjee pg. 37). If w = (w(s1), . . . , w(sn))⊤ and ǫ = (ǫ(s1), . . . , ǫ(sn))⊤, then a general model is V ar[w] = σ2H(φ) V ar[ǫ] = τ 2I

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51

SLIDE 61

Point-referenced spatial data Estimation

Consider point-referenced data at spatial locations s1, . . . , sn, model this data as Y (s) = µ(s) + w(s) + ǫ(s) If we constrain ourselves to isotropic models, the Mat´ ern class is suggested as a general tool (Banerjee pg. 37). If w = (w(s1), . . . , w(sn))⊤ and ǫ = (ǫ(s1), . . . , ǫ(sn))⊤, then a general model is V ar[w] = σ2H(φ) V ar[ǫ] = τ 2I where H is a correlation matrix with Hij = ρ(si − sj; φ) and ρ is a valid isotropic correlation function on Rr, i.e. Mat´ ern: ρ(u; ν, φ) = (u/φ)νKν(u/φ) 2ν−1Γ(ν) as defined in geoR:matern.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51

SLIDE 62

Point-referenced spatial data Estimation

Consider point-referenced data at spatial locations s1, . . . , sn, model this data as Y (s) = µ(s) + w(s) + ǫ(s) If we constrain ourselves to isotropic models, the Mat´ ern class is suggested as a general tool (Banerjee pg. 37). If w = (w(s1), . . . , w(sn))⊤ and ǫ = (ǫ(s1), . . . , ǫ(sn))⊤, then a general model is V ar[w] = σ2H(φ) V ar[ǫ] = τ 2I where H is a correlation matrix with Hij = ρ(si − sj; φ) and ρ is a valid isotropic correlation function on Rr, i.e. Mat´ ern: ρ(u; ν, φ) = (u/φ)νKν(u/φ) 2ν−1Γ(ν) as defined in geoR:matern. The overall mean is modeled separately and uses covariates x(s) via µ(s) = x(s)⊤β.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51

SLIDE 63

Point-referenced spatial data Estimation

Bayesian estimation for spatial random effects

Let θ = (β, σ2, τ 2, φ), then parameter estimates may be obtained from the posterior distribution: p(θ|y) ∝ p(y|θ)p(θ)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 17 / 51

SLIDE 64

Point-referenced spatial data Estimation

Bayesian estimation for spatial random effects

Let θ = (β, σ2, τ 2, φ), then parameter estimates may be obtained from the posterior distribution: p(θ|y) ∝ p(y|θ)p(θ) where Y |θ ∼ N(Xβ, σ2H(φ) + τ 2I).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 17 / 51

SLIDE 65

Point-referenced spatial data Estimation

Bayesian estimation for spatial random effects

Let θ = (β, σ2, τ 2, φ), then parameter estimates may be obtained from the posterior distribution: p(θ|y) ∝ p(y|θ)p(θ) where Y |θ ∼ N(Xβ, σ2H(φ) + τ 2I). Typically, independent priors are chosen so that p(θ) = p(β)p(σ2)p(τ 2)p(φ). As a general rule, non-informative priors can be chosen for β, e.g. p(β) ∝ 1.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 17 / 51

SLIDE 66

Point-referenced spatial data Estimation

Bayesian estimation for spatial random effects

Let θ = (β, σ2, τ 2, φ), then parameter estimates may be obtained from the posterior distribution: p(θ|y) ∝ p(y|θ)p(θ) where Y |θ ∼ N(Xβ, σ2H(φ) + τ 2I). Typically, independent priors are chosen so that p(θ) = p(β)p(σ2)p(τ 2)p(φ). As a general rule, non-informative priors can be chosen for β, e.g. p(β) ∝ 1. However, improper (or vague proper) priors for the variance parameters can lead to improper (or computationally improper) posteriors.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 17 / 51

SLIDE 67

Point-referenced spatial data Example

Diameter at breast height (DBH) for an experimental forest

50 100 150 200 100 200 300

East_m North_m DBH_cm

50 100 150 200

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 18 / 51

SLIDE 68

Point-referenced spatial data Example

Silver Fir Western Hemlock Douglas Fir Grand Fir Noble Fir 100 200 300 100 200 300 100 200 300 50 100 150 200 50 100 150 200

East_m North_m DBH_cm

50 100 150 200

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 19 / 51

SLIDE 69

Point-referenced spatial data Example

Interpolation of mean DBH (ignoring species)

50 100 150 200 250 300 350 50 100 150 200 Easting (m) Northing (m) 20 30 40 50 60 70 80 90

20 20 2 30 3 30 30 30 3 30 30 3 40 40 40 40 40 4 40 40 4 40 4 4 50 5 50 5 50 50 50 60 6 60 60 60 6 70 7 7 8

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 20 / 51

SLIDE 70

Point-referenced spatial data Example

Regression

## variog: computing omnidirectional variogram ## variofit: covariance model used is exponential ## variofit: weights used: equal ## variofit: minimisation function used: nls ## ## Call: ## lm(formula = DBH_cm ~ Species, data = d) ## ## Residuals: ## Min 1Q Median 3Q Max ## -78.423

9.969
3.561

10.924 118.277 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 89.423 1.303 68.629 <2e-16 *** ## SpeciesGrand Fir

51.598

4.133 -12.483 <2e-16 *** ## SpeciesNoble Fir

5.873

15.744

0.373

0.709 ## SpeciesSilver Fir

68.347

1.461 -46.784 <2e-16 *** ## SpeciesWestern Hemlock

48.062

1.636 -29.377 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 22.19 on 1950 degrees of freedom ## Multiple R-squared: 0.5332,Adjusted R-squared: 0.5323 ## F-statistic: 556.9 on 4 and 1950 DF, p-value: < 2.2e-16 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 21 / 51

SLIDE 71

Point-referenced spatial data Example

Variogram (exponential model)

20 40 60 80 200 400 600 800 1000 1200

DBH

distance semivariance 20 40 60 80 200 250 300 350 400 450 500

DBH residuals

distance semivariance

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 22 / 51

SLIDE 72

Point-referenced spatial data Example

Isotropy?

distance semivariance

100 200 300 400 500 20 40 60 80

45 90

20 40 60 80 100 200 300 400 500

135

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 23 / 51

SLIDE 73

Point-referenced spatial data Example

spBayes

p = nlevels(d$Species) r = spLM(DBH_cm ~ Species, data = d, coords = as.matrix(d[c('East_m','North_m')]), knots = c(6,6,.1), # for spatial prediction cov.model = 'exponential', starting = list(tau.sq = fit.DBH.resid$nugget, sigma.sq = fit.DBH.resid$cov.pars[1], phi = fit.DBH.resid$cov.pars[2]), tuning = list(tau.sq = 0.015, sigma.sq = 0.015, phi = 0.015), priors = list(beta.Norm = list(rep(0,p), diag(1000,p)), phi.Unif = c(3/1,3/0.1), sigma.sq.IG = c(2,200), tau.sq.IG = c(3,300)), n.samples = 10000, n.report = 200, verbose=TRUE) Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 24 / 51

SLIDE 74

Point-referenced spatial data Example ## ---------------------------------------- ## General model description ## ---------------------------------------- ## Model fit with 1955 observations. ## ## Number of covariates 5 (including intercept if specified). ## ## Using the exponential spatial correlation model. ## ## Using modified predictive process with 36 knots. ## ## Number of MCMC samples 10000. ## ## Priors and hyperpriors: ## beta normal: ## mu: 0.000 0.000 0.000 0.000 0.000 ## cov: ## 1000.000 0.000 0.000 0.000 0.000 ## 0.000 1000.000 0.000 0.000 0.000 ## 0.000 0.000 1000.000 0.000 0.000 ## 0.000 0.000 0.000 1000.000 0.000 ## 0.000 0.000 0.000 0.000 1000.000 ## ## sigma.sq IG hyperpriors shape=2.00000 and scale=200.00000 ## tau.sq IG hyperpriors shape=3.00000 and scale=300.00000 ## phi Unif hyperpriors a=3.00000 and b=30.00000 ## ------------------------------------------------- ## Sampling ## ------------------------------------------------- ## Sampled: 200 of 10000, 2.00% ## Report interval Metrop. Acceptance rate: 36.50% ## Overall Metrop. Acceptance rate: 36.50% ## ------------------------------------------------- ## Sampled: 400 of 10000, 4.00% ## Report interval Metrop. Acceptance rate: 36.50% ## Overall Metrop. Acceptance rate: 36.50% ## ------------------------------------------------- Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 25 / 51

SLIDE 75

Point-referenced spatial data Example

Traceplots

plot(r$p.theta.samples, density=FALSE)

2000 4000 6000 8000 100 300 500 Iterations

Trace of sigma.sq

2000 4000 6000 8000 100 300 500 Iterations

Trace of tau.sq

2000 4000 6000 8000 5 15 25 Iterations

Trace of phi

burnin = 500 nreps = dim(r$p.theta.samples)[1] Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 26 / 51

SLIDE 76

Point-referenced spatial data Example r$p.theta.samples[burnin:nreps,] %>% as.data.frame %>% GGally::ggpairs() + theme_bw() Corr: −0.992 Corr: 0.179 Corr: −0.175

sigma.sq tau.sq phi sigma.sq tau.sq phi 100 200 300 400 500 100 200 300 400 500 10 20 0.000 0.001 0.002 0.003 100 200 300 400 500 10 20

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 27 / 51

SLIDE 77

Point-referenced spatial data Example

Traceplot 2s

plot(r$p.beta.samples, density=FALSE)

2000 4000 6000 8000 10000 86 92 Iterations

Trace of (Intercept)

2000 4000 6000 8000 10000 −65 −40 Iterations

Trace of SpeciesGrand Fir

2000 4000 6000 8000 10000 −60 Iterations

Trace of SpeciesNoble Fir

2000 4000 6000 8000 10000 −72 −64 Iterations

Trace of SpeciesSilver Fir

2000 4000 6000 8000 10000 −54 −44 Iterations

Trace of SpeciesWestern Hemlock

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 28 / 51

SLIDE 78

Point-referenced spatial data Example

Summary statistics

summary(r$p.theta.samples) ## ## Iterations = 1:10000 ## Thinning interval = 1 ## Number of chains = 1 ## Sample size per chain = 10000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## sigma.sq 246.59 127.324 1.27324 31.224 ## tau.sq 244.94 126.725 1.26725 30.490 ## phi 17.24 7.204 0.07204 2.743 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## sigma.sq 43.200 133.07 238.13 366.61 446.7 ## tau.sq 51.208 124.40 252.81 355.24 449.0 ## phi 4.593 11.46 17.55 23.58 27.9 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 29 / 51

SLIDE 79

Point-referenced spatial data Example

Summary statistics 2

summary(r$p.beta.samples) ## ## Iterations = 1:10000 ## Thinning interval = 1 ## Number of chains = 1 ## Sample size per chain = 10000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## (Intercept) 89.011 1.291 0.01291 0.01291 ## SpeciesGrand Fir

50.361

4.125 0.04125 0.04125 ## SpeciesNoble Fir

4.406 14.034

0.14034 0.14034 ## SpeciesSilver Fir

67.923

1.458 0.01458 0.01458 ## SpeciesWestern Hemlock -47.605 1.631 0.01631 0.01631 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## (Intercept) 86.48 88.13 89.006 89.873 91.53 ## SpeciesGrand Fir

58.49 -53.11 -50.306 -47.617 -42.21

## SpeciesNoble Fir

32.44 -13.91
4.382

5.187 22.83 ## SpeciesSilver Fir

70.79 -68.90 -67.917 -66.944 -65.09

## SpeciesWestern Hemlock -50.79 -48.69 -47.615 -46.506 -44.44 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 30 / 51

SLIDE 80

Point-referenced spatial data Example

Spatial surface

If interest resides in w, draws can be obtained using the following relationship p(w|y) =

p(w|σ2, φ, y)p(σ2, φ|y)dσ2dφ

which suggests the following strategy:

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 31 / 51

SLIDE 81

Point-referenced spatial data Example

Spatial surface

If interest resides in w, draws can be obtained using the following relationship p(w|y) =

p(w|σ2, φ, y)p(σ2, φ|y)dσ2dφ

which suggests the following strategy:

1. Run the MCMC sampler to obtain draws (σ2, φ)(g) ∼ p(σ2, φ|y)
2. After burn-in and for g = 1, . . . , G, sample w(g) ∼ p(w|(σ2, φ)(g), y).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 31 / 51

SLIDE 82

Point-referenced spatial data Example

Prediction

For prediction at points s01, . . . , s0m and denoting Y0 = (Y (s01), . . . , Y (s0m))⊤ and design matrix X0 having rows x(s0j)⊤,

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51

SLIDE 83

Point-referenced spatial data Example

Prediction

For prediction at points s01, . . . , s0m and denoting Y0 = (Y (s01), . . . , Y (s0m))⊤ and design matrix X0 having rows x(s0j)⊤, we have the following relationship p(y0|y, X, X0) =

p(y0|y, θ, X0)p(θ|y, X)dθ

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51

SLIDE 84

Point-referenced spatial data Example

Prediction

For prediction at points s01, . . . , s0m and denoting Y0 = (Y (s01), . . . , Y (s0m))⊤ and design matrix X0 having rows x(s0j)⊤, we have the following relationship p(y0|y, X, X0) =

p(y0|y, θ, X0)p(θ|y, X)dθ ≈ 1

G

g=1

p(y0|y, θ(g), X0).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51

SLIDE 85

Point-referenced spatial data Example

Prediction

For prediction at points s01, . . . , s0m and denoting Y0 = (Y (s01), . . . , Y (s0m))⊤ and design matrix X0 having rows x(s0j)⊤, we have the following relationship p(y0|y, X, X0) =

p(y0|y, θ, X0)p(θ|y, X)dθ ≈ 1

G

g=1

p(y0|y, θ(g), X0). It is more common to take draws y(g) ∼ p(y0|y, θ(g), X0) and estimate the predictive distribution using p(y0|y, X, X0) ≈ 1 G

G

g=1

δy(g)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51

SLIDE 86

Point-referenced spatial data Example

Prediction

For prediction at points s01, . . . , s0m and denoting Y0 = (Y (s01), . . . , Y (s0m))⊤ and design matrix X0 having rows x(s0j)⊤, we have the following relationship p(y0|y, X, X0) =

p(y0|y, θ, X0)p(θ|y, X)dθ ≈ 1

G

g=1

p(y0|y, θ(g), X0). It is more common to take draws y(g) ∼ p(y0|y, θ(g), X0) and estimate the predictive distribution using p(y0|y, X, X0) ≈ 1 G

G

g=1

δy(g) where p(y0|y, θ, X0) has a conditional normal distribution.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51

SLIDE 87

Point-referenced spatial data Example

Predictions are not conditionally independent

Consider the joint distribution for y and y0 = y(s0) (a scalar for simplicity),

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51

SLIDE 88

Point-referenced spatial data Example

Predictions are not conditionally independent

Consider the joint distribution for y and y0 = y(s0) (a scalar for simplicity), then y y0

∼ N

Xβ X0β

,

Ω11 Ω12 Ω21 Ω22

Jarad Niemi (STAT615@ISU)

Bayesian Spatial Analysis November 9, 2017 33 / 51

SLIDE 89

Point-referenced spatial data Example

Predictions are not conditionally independent

Consider the joint distribution for y and y0 = y(s0) (a scalar for simplicity), then y y0

∼ N

Xβ X0β

,

Ω11 Ω12 Ω21 Ω22

where

Ω11 = σ2H(φ) + τ 2I Ω22 = σ2 + τ 2 Ω⊤

12

= σ2[ρ(d01; φ), . . . , ρ(d0n; φ)] and dij = ||si − sj||.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51

SLIDE 90

Point-referenced spatial data Example

Predictions are not conditionally independent

Consider the joint distribution for y and y0 = y(s0) (a scalar for simplicity), then y y0

∼ N

Xβ X0β

,

Ω11 Ω12 Ω21 Ω22

where

Ω11 = σ2H(φ) + τ 2I Ω22 = σ2 + τ 2 Ω⊤

12

= σ2[ρ(d01; φ), . . . , ρ(d0n; φ)] and dij = ||si − sj||. Thus y0|y, θ, X, X0 is normal with E[Y (s0)|y, θ, X, X0] = x⊤

0 β + Ω⊤ 12Ω−1 22 (y − Xβ)

V ar[Y (s0)|y, θ, X, X0] = σ2 + τ 2 − Ω⊤

12Ω−1 22 Ω12

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51

SLIDE 91

Point-referenced spatial data Example

Generalized linear spatial modeling

Let Y (s) be the response of interest with E[Y (s)] = g−1(x(s)⊤β + w(s)) where w(s) is our spatial random effect.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 34 / 51

SLIDE 92

Point-referenced spatial data Example

Generalized linear spatial modeling

Let Y (s) be the response of interest with E[Y (s)] = g−1(x(s)⊤β + w(s)) where w(s) is our spatial random effect. For example, Poisson regression Y (s) ∼ Po(ex(s)⊤β+w(s)).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 34 / 51

SLIDE 93

Point-referenced spatial data Example

Generalized linear spatial modeling

Let Y (s) be the response of interest with E[Y (s)] = g−1(x(s)⊤β + w(s)) where w(s) is our spatial random effect. For example, Poisson regression Y (s) ∼ Po(ex(s)⊤β+w(s)). For GLMs (other than linear models), w(s) cannot be integrated out

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 34 / 51

SLIDE 94

Point-referenced spatial data Example

Generalized linear spatial modeling

Let Y (s) be the response of interest with E[Y (s)] = g−1(x(s)⊤β + w(s)) where w(s) is our spatial random effect. For example, Poisson regression Y (s) ∼ Po(ex(s)⊤β+w(s)). For GLMs (other than linear models), w(s) cannot be integrated out and therefore a common MCMC strategy is

1. Sample β| . . ..
2. Sample w| . . ..
3. Sample θ| . . . (the spatial parameters [no nugget]).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 34 / 51

SLIDE 95

Areal-referenced data

Choropleth

<0.7 0.7−1 1−1.2 >1.2

MN Lung Cancer SMR

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 35 / 51

SLIDE 96

Areal-referenced data

Modeling areal units

Let Yi represent the SMR for lung cancer in MN county i.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51

SLIDE 97

Areal-referenced data

Modeling areal units

Let Yi represent the SMR for lung cancer in MN county i. Consider the model defined by conditional distributions: Yi|y−i ∼ N  

j∈ni

yj/mi, τ 2/mi   where ni indicates the neighbors of i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51

SLIDE 98

Areal-referenced data

Modeling areal units

Let Yi represent the SMR for lung cancer in MN county i. Consider the model defined by conditional distributions: Yi|y−i ∼ N  

j∈ni

yj/mi, τ 2/mi   where ni indicates the neighbors of i mi indicates the number of neighbors for i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51

SLIDE 99

Areal-referenced data

Modeling areal units

Let Yi represent the SMR for lung cancer in MN county i. Consider the model defined by conditional distributions: Yi|y−i ∼ N  

j∈ni

yj/mi, τ 2/mi   where ni indicates the neighbors of i mi indicates the number of neighbors for i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51

SLIDE 100

Areal-referenced data

Modeling areal units

Let Yi represent the SMR for lung cancer in MN county i. Consider the model defined by conditional distributions: Yi|y−i ∼ N  

j∈ni

yj/mi, τ 2/mi   where ni indicates the neighbors of i mi indicates the number of neighbors for i This defines a Markov Random Field.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51

SLIDE 101

Areal-referenced data Brook’s Lemma

Brook’s Lemma

It is clear that given p(y1, . . . , yn), the full conditionals, i.e. p(yi|y−i), are determined.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 37 / 51

SLIDE 102

Areal-referenced data Brook’s Lemma

Brook’s Lemma

It is clear that given p(y1, . . . , yn), the full conditionals, i.e. p(yi|y−i), are determined. Definition Brook’s Lemma states that p(y1, . . . , yn) p(y′

1, . . . , y′ n) = p(y1|y2, . . . , yn)

p(y′

1|y2, . . . , yn)·p(y2|y′ 1, y3, . . . , yn)

p(y′

2|y′ 1, y3, . . . , yn) · · · p(yn|y′ 1, . . . , y′ n−1)

p(y′

n|y′ 1, . . . , y′ n−1)

for all (y′

1, . . . , y′ n).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 37 / 51

SLIDE 103

Areal-referenced data Brook’s Lemma

Brook’s Lemma

It is clear that given p(y1, . . . , yn), the full conditionals, i.e. p(yi|y−i), are determined. Definition Brook’s Lemma states that p(y1, . . . , yn) p(y′

1, . . . , y′ n) = p(y1|y2, . . . , yn)

p(y′

1|y2, . . . , yn)·p(y2|y′ 1, y3, . . . , yn)

p(y′

2|y′ 1, y3, . . . , yn) · · · p(yn|y′ 1, . . . , y′ n−1)

p(y′

n|y′ 1, . . . , y′ n−1)

for all (y′

1, . . . , y′ n).

If

p(y′

1, . . . , y′ n) =

p(y1|y2, . . . , yn)

p(y′

1|y2, . . . , yn)

· p(y2|y′

1, y3, . . . , yn)

p(y′

2|y′ 1, y3, . . . , yn)

· · · p(yn|y′

1, . . . , y′ n−1)

p(y′

n|y′ 1, . . . , y′ n−1)

dy1, . . . , dyn < ∞

then p(y1, . . . , yn) is a proper joint distribution.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 37 / 51

SLIDE 104

Areal-referenced data Conditionally autoregressive models

Conditionally autoregressive models

More generally, we can consider Yi|y−i ∼ N  

j=i

bijyj, τ 2

i

  Through Brook’s Lemma, we have p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

where

B has elements bij D is diagonal with elements τ 2

i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 38 / 51

SLIDE 105

Areal-referenced data Conditionally autoregressive models

Conditionally autoregressive models

More generally, we can consider Yi|y−i ∼ N  

j=i

bijyj, τ 2

i

  Through Brook’s Lemma, we have p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

where

B has elements bij D is diagonal with elements τ 2

i

In order for D−1[I − B] to be symmetric, we need bij

τ 2

i = bji

τ 2

j for all i, j. Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 38 / 51

SLIDE 106

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 107

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 108

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 109

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 110

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 111

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 112

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point neighbors defined by those who are within distance δ

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 113

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point neighbors defined by those who are within distance δ K-nearest neighbors

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 114

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point neighbors defined by those who are within distance δ K-nearest neighbors

“distance”

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 115

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point neighbors defined by those who are within distance δ K-nearest neighbors

“distance”

inverse intercentroidal distance

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 116

Areal-referenced data Proximity matrix

Proximity matrix

Definition A proximity matrix is a an n × n matrix, W, with elements wii = 0 and wij representing the “distance” between unit i and unit j Common choices for wij are 1 if i is a neighbor of j and 0 otherwise

neighbors defined by those who share an edge neighbors defined by those who share a point neighbors defined by those who are within distance δ K-nearest neighbors

“distance”

inverse intercentroidal distance inverse minimum distance plus c

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 39 / 51

SLIDE 117

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

Jarad Niemi (STAT615@ISU)

Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 118

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+,

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 119

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

Jarad Niemi (STAT615@ISU)

Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 120

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

where

W is our proximity matrix and

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 121

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

where

W is our proximity matrix and Dw has diagonal elements wi+

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 122

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

where

W is our proximity matrix and Dw has diagonal elements wi+

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 123

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

where

W is our proximity matrix and Dw has diagonal elements wi+ This can be rewritten as p(y1, . . . , yn) ∝ exp  − 1 2τ 2

i=j

wij(yi − yj)2  

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 124

Areal-referenced data Proximity matrix

Intrinsicially autoregressive model

Recall p(y1, . . . , yn) ∝ exp

−1

2y⊤D−1[I − B]y

if we set wi+ = n

j=1 wij, bij = wij/wi+, and τ 2 i = τ 2/wi+, we have

p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − W]y

where

W is our proximity matrix and Dw has diagonal elements wi+ This can be rewritten as p(y1, . . . , yn) ∝ exp  − 1 2τ 2

i=j

wij(yi − yj)2   This is called the intrinsically autoregressive model.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 40 / 51

SLIDE 125

Areal-referenced data Proper CAR models

Proper CAR models

To make this proper, p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − ρW]y

with

ρ ∈ (1/λ(n), 1/λ(1)) where

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 41 / 51

SLIDE 126

Areal-referenced data Proper CAR models

Proper CAR models

To make this proper, p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − ρW]y

with

ρ ∈ (1/λ(n), 1/λ(1)) where λ(1) < · · · < λ(n) are the ordered eigenvalues of D−1/2

w

WD−1/2

w

.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 41 / 51

SLIDE 127

Areal-referenced data Proper CAR models

Proper CAR models

To make this proper, p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − ρW]y

with

ρ ∈ (1/λ(n), 1/λ(1)) where λ(1) < · · · < λ(n) are the ordered eigenvalues of D−1/2

w

WD−1/2

w

.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 41 / 51

SLIDE 128

Areal-referenced data Proper CAR models

Proper CAR models

To make this proper, p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − ρW]y

with

ρ ∈ (1/λ(n), 1/λ(1)) where λ(1) < · · · < λ(n) are the ordered eigenvalues of D−1/2

w

WD−1/2

w

. The full conditionals are Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+  

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 41 / 51

SLIDE 129

Areal-referenced data Proper CAR models

Proper CAR models

To make this proper, p(y1, . . . , yn) ∝ exp

− 1

2τ 2 y⊤[Dw − ρW]y

with

ρ ∈ (1/λ(n), 1/λ(1)) where λ(1) < · · · < λ(n) are the ordered eigenvalues of D−1/2

w

WD−1/2

w

. The full conditionals are Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   a prior for ρ that induces a reasonable amount of spatial association should put most of its mass near 1.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 41 / 51

SLIDE 130

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model:

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 131

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 132

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association. ρ is a proportional reaction to the weighted average of its neighbors.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 133

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association. ρ is a proportional reaction to the weighted average of its neighbors. If ρ < 1, then the expected value of the current location is less than the weighted average of its neighbors.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 134

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association. ρ is a proportional reaction to the weighted average of its neighbors. If ρ < 1, then the expected value of the current location is less than the weighted average of its neighbors. If ρ = 0, then we have conditional independence.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 135

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association. ρ is a proportional reaction to the weighted average of its neighbors. If ρ < 1, then the expected value of the current location is less than the weighted average of its neighbors. If ρ = 0, then we have conditional independence. But the variance decreases with the number of neighbors which is perplexing.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 136

Areal-referenced data Proper CAR models

Issues with the proper CAR

The full condition for the proper CAR, i.e. Yi|y−i ∼ N  ρ

j=i

wijyj/wi+, τ 2/wi+   indicates some issues with this model: τ 2 does not play a role in spatial association. ρ is a proportional reaction to the weighted average of its neighbors. If ρ < 1, then the expected value of the current location is less than the weighted average of its neighbors. If ρ = 0, then we have conditional independence. But the variance decreases with the number of neighbors which is perplexing. ρ needs to be very close to 1 to obtain a consequential amount of spatial association.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 42 / 51

SLIDE 137

Areal-referenced data Proper CAR models

Dealing with ρ in the proper CAR

Choose ρ so the CAR model is proper

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 43 / 51

SLIDE 138

Areal-referenced data Proper CAR models

Dealing with ρ in the proper CAR

Choose ρ so the CAR model is proper Choose ρ = 1 (improper IAR model) and constrain n

i=1 Yi = 0

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 43 / 51

SLIDE 139

Areal-referenced data Proper CAR models

Dealing with ρ in the proper CAR

Choose ρ so the CAR model is proper Choose ρ = 1 (improper IAR model) and constrain n

i=1 Yi = 0

Choose ρ = 1 and estimate a mean (remove mean from the fixed effect)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 43 / 51

SLIDE 140

Areal-referenced data Proper CAR models

Dealing with ρ in the proper CAR

Choose ρ so the CAR model is proper Choose ρ = 1 (improper IAR model) and constrain n

i=1 Yi = 0

Choose ρ = 1 and estimate a mean (remove mean from the fixed effect) Let ρ ∼ Be(18, 2) (Banerjee pg 164) and estimate it.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 43 / 51

SLIDE 141

Areal-referenced data Leroux CAR

Leroux CAR

The Leroux et al. (1999) CAR tries to ameliorate these issues.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 44 / 51

SLIDE 142

Areal-referenced data Leroux CAR

Leroux CAR

The Leroux et al. (1999) CAR tries to ameliorate these issues. The joint distribution is Y ∼ N(0, τ 2[ρ(Dw − W) + (1 − rho)I]−1)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 44 / 51

SLIDE 143

Areal-referenced data Leroux CAR

Leroux CAR

The Leroux et al. (1999) CAR tries to ameliorate these issues. The joint distribution is Y ∼ N(0, τ 2[ρ(Dw − W) + (1 − rho)I]−1) and the conditional distributions are Yi|y−i ∼ N

ρ

j=i wijyj

rho

j=i wijyj + 1 − ρ,

τ 2 rho

j=i wijyj + 1 − ρ

.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 44 / 51

SLIDE 144

Areal-referenced data Leroux CAR

Leroux CAR

The Leroux et al. (1999) CAR tries to ameliorate these issues. The joint distribution is Y ∼ N(0, τ 2[ρ(Dw − W) + (1 − rho)I]−1) and the conditional distributions are Yi|y−i ∼ N

ρ

j=i wijyj

rho

j=i wijyj + 1 − ρ,

τ 2 rho

j=i wijyj + 1 − ρ

.

This distribution is proper so long as 0 ≤ ρ < 1.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 44 / 51

SLIDE 145

Areal-referenced data Leroux CAR

Leroux CAR

The Leroux et al. (1999) CAR tries to ameliorate these issues. The joint distribution is Y ∼ N(0, τ 2[ρ(Dw − W) + (1 − rho)I]−1) and the conditional distributions are Yi|y−i ∼ N

ρ

j=i wijyj

rho

j=i wijyj + 1 − ρ,

τ 2 rho

j=i wijyj + 1 − ρ

.

This distribution is proper so long as 0 ≤ ρ < 1. Lee (2011) argued that this CAR should be preferred for a variety of reasons.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 44 / 51

SLIDE 146

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 147

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 148

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 149

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 150

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i then a possible model is Yi = X⊤

i β + ωs[i] + ǫi

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 151

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i then a possible model is Yi = X⊤

i β + ωs[i] + ǫi

where ǫi

ind

∼ N(0, σ2) is noise

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 152

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i then a possible model is Yi = X⊤

i β + ωs[i] + ǫi

where ǫi

ind

∼ N(0, σ2) is noise ωs is the spatial random effect associated with areal unit s, e.g.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 153

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i then a possible model is Yi = X⊤

i β + ωs[i] + ǫi

where ǫi

ind

∼ N(0, σ2) is noise ωs is the spatial random effect associated with areal unit s, e.g.

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 154

Areal-referenced data Random effect model

CAR as a model for random effects

Let Yi represent the (continuous) response for observation i Xi represent explanatory variables for observation i s[i] represent the areal unit for observation i then a possible model is Yi = X⊤

i β + ωs[i] + ǫi

where ǫi

ind

∼ N(0, σ2) is noise ωs is the spatial random effect associated with areal unit s, e.g. p(ω1, . . . , ωS) ∝ exp

− 1

2τ 2 ω⊤[Dw − ρW]ω

Jarad Niemi (STAT615@ISU)

Bayesian Spatial Analysis November 9, 2017 45 / 51

SLIDE 155

Areal-referenced data Housing price in Glasgow

650000 660000 670000 680000 220000 230000 240000 250000 260000 270000

0 5000 m

50 100 150 200 250 300 350

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 46 / 51

SLIDE 156

Areal-referenced data Housing price in Glasgow

Housing price model

Let Yi be the logarithm of the median home price in each Intermediate Geography (IG) to the north of the river Clude in the Greater Glasgow and Clyde health board,

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 47 / 51

SLIDE 157

Areal-referenced data Housing price in Glasgow

Housing price model

Let Yi be the logarithm of the median home price in each Intermediate Geography (IG) to the north of the river Clude in the Greater Glasgow and Clyde health board, use explanatory variables

crime: crime rate (number of crimes per 10,000 people) in each IG (logged), rooms: median number of rooms in a property in each IG, type: predominant property type in each IG with levels: detached, flat, semi, terrace, sales: percentage of properties that sold in each IG in a year, and driveshop: average time taken to drive to a shopping centre in minutes (logged).

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 47 / 51

SLIDE 158

Areal-referenced data Housing price in Glasgow

Housing price model

Assume Yi

ind

∼ N(Xiβ + ωi, ν2)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 48 / 51

SLIDE 159

Areal-referenced data Housing price in Glasgow

Housing price model

Assume Yi

ind

∼ N(Xiβ + ωi, ν2)

r, alternatively,

Yi = Xβ + ωi + ǫi, ǫi

ind

∼ N(0, ν2)

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 48 / 51

SLIDE 160

Areal-referenced data Housing price in Glasgow

Housing price model

Assume Yi

ind

∼ N(Xiβ + ωi, ν2)

r, alternatively,

Yi = Xβ + ωi + ǫi, ǫi

ind

∼ N(0, ν2) where β are the regression parameters and

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 48 / 51

SLIDE 161

Areal-referenced data Housing price in Glasgow

Housing price model

Assume Yi

ind

∼ N(Xiβ + ωi, ν2)

r, alternatively,

Yi = Xβ + ωi + ǫi, ǫi

ind

∼ N(0, ν2) where β are the regression parameters and ωi are assumed to come from an intrinsic CAR model with proximity matrix indicating those regions that share a border

Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 48 / 51

SLIDE 162

Areal-referenced data Housing price in Glasgow propertydata.spatial@data$logprice <- log(propertydata.spatial@data$price) propertydata.spatial@data$logdriveshop <- log(propertydata.spatial@data$driveshop) ################################################### ### code chunk number 9: CARBayes.Rnw:495-498 ################################################### library(splines) form <- logprice~ns(crime,3)+rooms+sales+factor(type) + logdriveshop model <- lm(formula=form, data=propertydata.spatial@data) ################################################### ### code chunk number 10: CARBayes.Rnw:505-510 ################################################### library(spdep) W.nb <- poly2nb(propertydata.spatial, row.names = rownames(propertydata.spatial@data)) W.list <- nb2listw(W.nb, style="B") resid.model <- residuals(model) moran.mc(x=resid.model, listw=W.list, nsim=1000) ## ## Monte-Carlo simulation of Moran I ## ## data: resid.model ## weights: W.list ## number of simulations + 1: 1001 ## ## statistic = 0.2733, observed rank = 1001, p-value = 0.000999 ## alternative hypothesis: greater Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 49 / 51

SLIDE 163

Areal-referenced data Housing price in Glasgow ## ## Call: ## lm(formula = form, data = propertydata.spatial@data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.91319 -0.15992 0.00136 0.15647 0.81675 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.436135 0.157971 28.082 < 2e-16 *** ## ns(crime, 3)1

0.358967

0.089006

4.033 7.24e-05 ***

## ns(crime, 3)2

0.617084

0.165152

3.736 0.000229 ***

## ns(crime, 3)3

0.299454

0.126516

2.367 0.018670 *

## rooms 0.193827 0.029268 6.623 2.02e-10 *** ## sales 0.002034 0.000362 5.619 4.93e-08 *** ## factor(type)flat

0.215967

0.066412

3.252 0.001298 **

## factor(type)semi

0.153610

0.057750

2.660 0.008301 **

## factor(type)terrace -0.280023 0.072634

3.855 0.000146 ***

## logdriveshop

0.089084

0.025588

3.482 0.000585 ***

## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.2243 on 260 degrees of freedom ## Multiple R-squared: 0.6206,Adjusted R-squared: 0.6075 ## F-statistic: 47.26 on 9 and 260 DF, p-value: < 2.2e-16 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 50 / 51

SLIDE 164

Areal-referenced data Housing price in Glasgow

Bayesian analysis using Leroux CAR

W <- nb2mat(W.nb, style = "B") model.spatial <- S.CARleroux(formula = form, data = propertydata.spatial@data, family = "gaussian", W = W, burnin = 20000, n.sample = 120000, thin = 10) ## Setting up the model ## Error: the formula inputted contains an error, e.g the variables may be different lengths. Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 51 / 51