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
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
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
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
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
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
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 ) = || u ||→∞ γ ( u ) − γ ( h ) . lim Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 9 / 51
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 ) = || u ||→∞ γ ( u ) − γ ( h ) . lim 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
Point-referenced spatial data Assumptions Covariance functions for isotropic models Model Covariance function, C ( t ) Semivariogram, γ ( t ) τ 2 + σ 2 t � if t > 0 Linear C ( t ) does not exist γ ( t ) = 0 otherwise τ 2 + σ 2 0 if σ 2 � 2 ( φt ) 3 � 1 − 3 2 φt + 1 τ 2 + σ 2 � 2 ( φt ) 3 � 3 2 φt − 1 Spherical C ( t ) = γ ( t ) = 0 τ 2 + σ 2 0 otherwise σ 2 exp( − φt ) τ 2 + σ 2 [1 − exp( − φt )] � � t > Exponential C ( t ) = γ ( t ) = τ 2 + σ 2 0 otherwise σ 2 exp( −| φt | p ) τ 2 + σ 2 [1 − exp( −| φt | p )] � � t Powered exponential C ( t ) = γ ( t ) = τ 2 + σ 2 0 otherwise σ 2 exp( − φ 2 t 2 ) � τ 2 + σ 2 � 1 − exp( − φ 2 t 2 ) � � Gaussian C ( t ) = γ ( t ) = τ 2 + σ 2 0 � t 2 � t 2 σ 2 � τ 2 + σ 2 1 − t > 0 (1+ φ 2) (1+ φ 2) Rational quadratic C ( t ) = γ ( t ) = τ 2 + σ 2 0 otherwise � σ 2 sin( φt ) � τ 2 + σ 2 � 1 − sin( φt ) � t > 0 Wave C ( t ) = φt γ ( t ) = φt τ 2 + σ 2 0 otherwise τ 2 + σ 2 t λ � t > 0 Power law C ( t ) does not exist γ ( t ) = 0 otherwise 2 ν − 1Γ( ν ) (2 √ vtφ ) ν K ν (2 √ νtφ ) 1 − ( 2 √ νtφ ) ν σ 2 � 2 ν − 1Γ( ν ) (2 √ v τ 2 + σ 2 Mat´ ern C ( t ) = γ ( t ) = τ 2 + σ 2 0 σ 2 (1 + φt ) exp( − φt ) τ 2 + σ 2 [1 − (1 + φt ) exp( − φt � � Mat´ ern ( ν = 3 / 2 ) C ( t ) = γ ( t ) = τ 2 + σ 2 0 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 10 / 51
Point-referenced spatial data Assumptions Spherical semivariogram Spherical semivariogram 1.2 1.0 partial sill ( σ 2 ) 0.8 sill ( τ 2 + σ 2 ) γ ( t ) 0.6 range (1 φ ) 0.4 0.2 nugget ( τ 2 ) 0.0 0.0 0.5 1.0 1.5 2.0 t Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 11 / 51
Point-referenced spatial data Assumptions Mat´ ern Perhaps the most important isotropic process is the Mat´ ern process with covariance 2 ν − 1 Γ( ν ) (2 √ νtφ ) ν K ν (2 √ νtφ ) � σ 2 t > 0 C ( t ) = τ 2 + σ 2 t = 0 and variogram 1 − (2 √ νtφ ) ν 2 ν − 1 Γ( ν ) K ν (2 √ νtφ ) � τ 2 + σ 2 � � t > 0 γ ( t ) = 0 otherwise Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 12 / 51
Point-referenced spatial data Assumptions Mat´ ern Perhaps the most important isotropic process is the Mat´ ern process with covariance 2 ν − 1 Γ( ν ) (2 √ νtφ ) ν K ν (2 √ νtφ ) � σ 2 t > 0 C ( t ) = τ 2 + σ 2 t = 0 and variogram 1 − (2 √ νtφ ) ν 2 ν − 1 Γ( ν ) K ν (2 √ νtφ ) � τ 2 + σ 2 � � t > 0 γ ( t ) = 0 otherwise 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
Point-referenced spatial data Assumptions Mat´ ern Perhaps the most important isotropic process is the Mat´ ern process with covariance 2 ν − 1 Γ( ν ) (2 √ νtφ ) ν K ν (2 √ νtφ ) � σ 2 t > 0 C ( t ) = τ 2 + σ 2 t = 0 and variogram 1 − (2 √ νtφ ) ν 2 ν − 1 Γ( ν ) K ν (2 √ νtφ ) � τ 2 + σ 2 � � t > 0 γ ( t ) = 0 otherwise 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
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 { s 1 , . . . , s n } and any h ∈ R d , ( Y ( s 1 ) , . . . , Y ( s n )) ⊤ d = ( Y ( s 1 + h ) , . . . , Y ( s n + h )) ⊤ where d = means equal in distribution. Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 13 / 51
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 { s 1 , . . . , s n } and any h ∈ R d , ( Y ( s 1 ) , . . . , Y ( s n )) ⊤ d = ( Y ( s 1 + h ) , . . . , Y ( s n + 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
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 { s 1 , . . . , s n } and any h ∈ R d , ( Y ( s 1 ) , . . . , Y ( s n )) ⊤ d = ( Y ( s 1 + h ) , . . . , Y ( s n + 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
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 { s 1 , . . . , s n } , Y = ( Y ( s 1 ) , . . . , Y ( s n )) ⊤ has a multivariate normal distribution. Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 14 / 51
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 { s 1 , . . . , s n } , Y = ( Y ( s 1 ) , . . . , Y ( s n )) ⊤ 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
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . Let’s assume the data arise from a Gaussian process Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . 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
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . 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
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . 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
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . 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
Point-referenced spatial data Estimation Bayesian estimation of Gaussian process parameters Suppose we observe data at some locations s 1 , . . . , s n . Collectively, we have y = ( y ( s 1 ) , . . . , y ( s n )) . 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. || s 1 − s 2 || . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 15 / 51
Point-referenced spatial data Estimation Consider point-referenced data at spatial locations s 1 , . . . , s n , 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
Point-referenced spatial data Estimation Consider point-referenced data at spatial locations s 1 , . . . , s n , 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 ( s 1 ) , . . . , w ( s n )) ⊤ and ǫ = ( ǫ ( s 1 ) , . . . , ǫ ( s n )) ⊤ , then a general model is Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51
Point-referenced spatial data Estimation Consider point-referenced data at spatial locations s 1 , . . . , s n , 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 ( s 1 ) , . . . , w ( s n )) ⊤ and ǫ = ( ǫ ( s 1 ) , . . . , ǫ ( s n )) ⊤ , then a general model is V ar [ w ] = σ 2 H ( φ ) V ar [ ǫ ] = τ 2 I Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 16 / 51
Point-referenced spatial data Estimation Consider point-referenced data at spatial locations s 1 , . . . , s n , 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 ( s 1 ) , . . . , w ( s n )) ⊤ and ǫ = ( ǫ ( s 1 ) , . . . , ǫ ( s n )) ⊤ , then a general model is V ar [ w ] = σ 2 H ( φ ) V ar [ ǫ ] = τ 2 I where H is a correlation matrix with H ij = ρ ( s i − s j ; φ ) and ρ is a valid isotropic correlation function on R r , 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
Point-referenced spatial data Estimation Consider point-referenced data at spatial locations s 1 , . . . , s n , 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 ( s 1 ) , . . . , w ( s n )) ⊤ and ǫ = ( ǫ ( s 1 ) , . . . , ǫ ( s n )) ⊤ , then a general model is V ar [ w ] = σ 2 H ( φ ) V ar [ ǫ ] = τ 2 I where H is a correlation matrix with H ij = ρ ( s i − s j ; φ ) and ρ is a valid isotropic correlation function on R r , 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
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
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β, σ 2 H ( φ ) + τ 2 I) . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 17 / 51
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β, σ 2 H ( φ ) + τ 2 I) . 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
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β, σ 2 H ( φ ) + τ 2 I) . 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
Point-referenced spatial data Example Diameter at breast height (DBH) for an experimental forest 200 150 DBH_cm North_m 50 100 150 100 200 50 0 0 100 200 300 East_m Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 18 / 51
Point-referenced spatial data Example Douglas Fir Grand Fir Noble Fir 200 150 100 50 DBH_cm 50 North_m 0 0 100 200 300 100 Silver Fir Western Hemlock 150 200 200 150 100 50 0 0 100 200 300 0 100 200 300 East_m Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 19 / 51
Point-referenced spatial data Example Interpolation of mean DBH (ignoring species) 90 3 0 40 200 5 0 80 50 4 0 60 50 60 70 40 3 0 30 150 40 Northing (m) 4 0 20 60 6 0 40 60 50 30 40 50 50 4 0 100 5 0 40 30 20 40 0 2 30 40 30 50 0 4 6 0 7 0 0 7 30 60 30 50 0 8 40 70 20 3 0 0 0 50 100 150 200 250 300 350 Easting (m) Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 20 / 51
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
Point-referenced spatial data Example Variogram (exponential model) DBH DBH residuals 1200 500 450 1000 400 semivariance semivariance 800 350 600 300 400 250 200 200 0 20 40 60 80 0 20 40 60 80 distance distance Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 22 / 51
Point-referenced spatial data Example Isotropy? 20 40 60 80 90 135 500 400 300 200 semivariance 100 0 45 500 400 300 200 100 20 40 60 80 distance Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 23 / 51
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
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 ## -------------------------------------------------
Point-referenced spatial data Example Traceplots plot(r$p.theta.samples, density=FALSE) Trace of sigma.sq Trace of tau.sq 500 500 300 300 100 100 0 2000 4000 6000 8000 0 2000 4000 6000 8000 Iterations Iterations Trace of phi 25 15 5 0 2000 4000 6000 8000 Iterations Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 26 / 51 burnin = 500 nreps = dim(r$p.theta.samples)[1]
Point-referenced spatial data Example r$p.theta.samples[burnin:nreps,] %>% as.data.frame %>% GGally::ggpairs() + theme_bw() sigma.sq tau.sq phi 0.003 sigma.sq 0.002 Corr: Corr: −0.992 0.179 0.001 0.000 500 400 300 Corr: tau.sq −0.175 200 100 20 phi 10 0 100 200 300 400 500 100 200 300 400 500 10 20 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 27 / 51
Point-referenced spatial data Example Traceplot 2s plot(r$p.beta.samples, density=FALSE) Trace of (Intercept) Trace of SpeciesGrand Fir −40 92 86 −65 0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000 Iterations Iterations Trace of SpeciesNoble Fir Trace of SpeciesSilver Fir −64 0 −72 −60 0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000 Iterations Iterations Trace of SpeciesWestern Hemlock −44 −54 0 2000 4000 6000 8000 10000 Iterations Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 28 / 51
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
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
Point-referenced spatial data Example Spatial surface If interest resides in w , draws can be obtained using the following relationship � p ( w | σ 2 , φ, y ) p ( σ 2 , φ | y ) dσ 2 dφ p ( w | y ) = which suggests the following strategy: Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 31 / 51
Point-referenced spatial data Example Spatial surface If interest resides in w , draws can be obtained using the following relationship � p ( w | σ 2 , φ, y ) p ( σ 2 , φ | y ) dσ 2 dφ p ( w | y ) = 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
Point-referenced spatial data Example Prediction For prediction at points s 01 , . . . , s 0 m and denoting Y 0 = ( Y ( s 01 ) , . . . , Y ( s 0 m )) ⊤ and design matrix X 0 having rows x ( s 0 j ) ⊤ , Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51
Point-referenced spatial data Example Prediction For prediction at points s 01 , . . . , s 0 m and denoting Y 0 = ( Y ( s 01 ) , . . . , Y ( s 0 m )) ⊤ and design matrix X 0 having rows x ( s 0 j ) ⊤ , we have the following relationship � p ( y 0 | y, X, X 0 ) = p ( y 0 | y, θ, X 0 ) p ( θ | y, X ) dθ Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51
Point-referenced spatial data Example Prediction For prediction at points s 01 , . . . , s 0 m and denoting Y 0 = ( Y ( s 01 ) , . . . , Y ( s 0 m )) ⊤ and design matrix X 0 having rows x ( s 0 j ) ⊤ , we have the following relationship G � p ( y 0 | y, θ, X 0 ) p ( θ | y, X ) dθ ≈ 1 � p ( y 0 | y, θ ( g ) , X 0 ) . p ( y 0 | y, X, X 0 ) = G g =1 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51
Point-referenced spatial data Example Prediction For prediction at points s 01 , . . . , s 0 m and denoting Y 0 = ( Y ( s 01 ) , . . . , Y ( s 0 m )) ⊤ and design matrix X 0 having rows x ( s 0 j ) ⊤ , we have the following relationship G � p ( y 0 | y, θ, X 0 ) p ( θ | y, X ) dθ ≈ 1 � p ( y 0 | y, θ ( g ) , X 0 ) . p ( y 0 | y, X, X 0 ) = G g =1 It is more common to take draws y ( g ) ∼ p ( y 0 | y, θ ( g ) , X 0 ) and estimate the 0 predictive distribution using G p ( y 0 | y, X, X 0 ) ≈ 1 � δ y ( g ) G 0 g =1 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51
Point-referenced spatial data Example Prediction For prediction at points s 01 , . . . , s 0 m and denoting Y 0 = ( Y ( s 01 ) , . . . , Y ( s 0 m )) ⊤ and design matrix X 0 having rows x ( s 0 j ) ⊤ , we have the following relationship G � p ( y 0 | y, θ, X 0 ) p ( θ | y, X ) dθ ≈ 1 � p ( y 0 | y, θ ( g ) , X 0 ) . p ( y 0 | y, X, X 0 ) = G g =1 It is more common to take draws y ( g ) ∼ p ( y 0 | y, θ ( g ) , X 0 ) and estimate the 0 predictive distribution using G p ( y 0 | y, X, X 0 ) ≈ 1 � δ y ( g ) G 0 g =1 where p ( y 0 | y, θ, X 0 ) has a conditional normal distribution. Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 32 / 51
Point-referenced spatial data Example Predictions are not conditionally independent Consider the joint distribution for y and y 0 = y ( s 0 ) (a scalar for simplicity), Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51
Point-referenced spatial data Example Predictions are not conditionally independent Consider the joint distribution for y and y 0 = y ( s 0 ) (a scalar for simplicity), then � y �� Xβ � Ω 11 � � �� Ω 12 ∼ N , y 0 X 0 β Ω 21 Ω 22 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51
Point-referenced spatial data Example Predictions are not conditionally independent Consider the joint distribution for y and y 0 = y ( s 0 ) (a scalar for simplicity), then � y �� Xβ � Ω 11 � � �� Ω 12 ∼ N , y 0 X 0 β Ω 21 Ω 22 where = σ 2 H ( φ ) + τ 2 I Ω 11 = σ 2 + τ 2 Ω 22 Ω ⊤ = σ 2 [ ρ ( d 01 ; φ ) , . . . , ρ ( d 0 n ; φ )] 12 and d ij = || s i − s j || . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51
Point-referenced spatial data Example Predictions are not conditionally independent Consider the joint distribution for y and y 0 = y ( s 0 ) (a scalar for simplicity), then � y �� Xβ � Ω 11 � � �� Ω 12 ∼ N , y 0 X 0 β Ω 21 Ω 22 where = σ 2 H ( φ ) + τ 2 I Ω 11 = σ 2 + τ 2 Ω 22 Ω ⊤ = σ 2 [ ρ ( d 01 ; φ ) , . . . , ρ ( d 0 n ; φ )] 12 and d ij = || s i − s j || . Thus y 0 | y, θ, X, X 0 is normal with = x ⊤ 0 β + Ω ⊤ 12 Ω − 1 E [ Y ( s 0 ) | y, θ, X, X 0 ] 22 ( y − Xβ ) = σ 2 + τ 2 − Ω ⊤ 12 Ω − 1 V ar [ Y ( s 0 ) | y, θ, X, X 0 ] 22 Ω 12 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 33 / 51
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
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 ( e x ( s ) ⊤ β + w ( s ) ) . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 34 / 51
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 ( e x ( 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
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 ( e x ( 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
Areal-referenced data Choropleth MN Lung Cancer SMR <0.7 0.7−1 1−1.2 >1.2 Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 35 / 51
Areal-referenced data Modeling areal units Let Y i represent the SMR for lung cancer in MN county i . Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51
Areal-referenced data Modeling areal units Let Y i represent the SMR for lung cancer in MN county i . Consider the model defined by conditional distributions: � y j /m i , τ 2 /m i Y i | y − i ∼ N j ∈ n i where n i indicates the neighbors of i Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51
Areal-referenced data Modeling areal units Let Y i represent the SMR for lung cancer in MN county i . Consider the model defined by conditional distributions: � y j /m i , τ 2 /m i Y i | y − i ∼ N j ∈ n i where n i indicates the neighbors of i m i indicates the number of neighbors for i Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51
Areal-referenced data Modeling areal units Let Y i represent the SMR for lung cancer in MN county i . Consider the model defined by conditional distributions: � y j /m i , τ 2 /m i Y i | y − i ∼ N j ∈ n i where n i indicates the neighbors of i m i indicates the number of neighbors for i Jarad Niemi (STAT615@ISU) Bayesian Spatial Analysis November 9, 2017 36 / 51
Areal-referenced data Modeling areal units Let Y i represent the SMR for lung cancer in MN county i . Consider the model defined by conditional distributions: � y j /m i , τ 2 /m i Y i | y − i ∼ N j ∈ n i where n i indicates the neighbors of i m i 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
Recommend
More recommend
