Spatially structured random effects

by Daniel A. Griffith Ashbel Smith Professor of Geospatial Information Sciences

ABSTRACT

Researchers increasingly are employing

random effects modeling

to analyze data. When data are georeferenced, a random effect term needs to be

spatially structured

in order to account for spatial autocorrelation. Spatial structuring can be achieved in various ways, including the use of

semivariogram

,

spatial autoregressive

models. GeoBUGS as with the SAS , and

spatial filter

models. SAS implements the semivariogram option for linear mixed implements the spatial autoregressive option for either linear or generalized linear mixed modeling. Recently developed spatial filtering methodology can be used in either case, as well generalized linear mix model procedure, and furnishes one means of estimating space-time mixed models. This presentation summarizes comparisons of these three forms of spatial structuring, illustrating implementations with selected ecological data for the municipalities of Puerto Rico .

From Legendre

Spatial structures in communities indicate that some process hasbeen at work to create them. Two families of mechanisms • cangenerate spatial structures in communities: Autocorrelation model : the spatial structures are generated • by thespecies assemblage themselves (response variables). Induced spatial dependence model : forcing (explanatory) variablesare responsible for the spatial structures found in the speciesassemblage. They represent environmental or biotic control of thespecies assemblages, or historical dynamics To understand the mechanisms that generate these structures, we needto explicitly incorporate the spatial community structures, at all scales, into the statistical model.

Spatial autocorrelation (SA) is technically defined as the dependence, due to geographic proximity, present in the residuals of a [regression-type] model of a response variable y whicht akes into account all deterministic effects due to forcing variables.

Spatial autocorrelation can be interpreted in different ways

As a spatial process mechanism

– the cartoon

As a diagnostic tool –

the Cliff-Ord Eire example (the model specification should be nonlinear)

As a nuisance parameter –

eliminating spatial dependency to avoid statistical complications

As a spatial spillover effect –

georeferencing of pediatric lead poisoning cases in Syracuse, NY

As an outcome of areal unit demarcation –

the modifiable areal unit problem (MAUP)

As redundant information

– spatial sampling; map interpolation

As map pattern

– spatial filtering (to be discussed in this course)

As a missing variables indicator/surrogate

– a possible implication of spatial filtering

As self-correlation

– what is discussed next

The magic box is a physical model of spatial autocorrelation

The permutation perspective

The

SA SIM

game http://www.nku.edu/~longa/cgi-bin/cgi-tcl-examples/generic/SA/SA.cgi

Measures of spatial autocorrelation

MC : Moran Coefficient

;

GR : Geary Ratio

;

semivariogram n MC  i  1 n n  j  1 c ij i n   1 (y i  y ) n  j  1 c ij (y j  y ) i n   1 (y i  y ) 2 GR n 1  2 i  1 n n  j  1 c ij i  1 j n n   1 c ij (y i  y j ) 2 i n   1 (y i  y ) 2

γ( d

k

)

 i n   k 1

(y

i 

y

j

)

2

2 * n

k

f( d

k

), k



1,2,...,



K

Spherical Exponential Bessel function (1 st order, 2 nd kind)

Georeferenced data scatterplots

• The horizontal axis is the measurement scale for some attribute variable • The vertical axis is the measurement scale for neighboring values (

topological distance

-based ) of the same attribute variable

OR

• The horizontal axis is (usually)

Euclidean distance

between geocoded locations • The vertical axis is the measurement scale for geographic variability

Describing a scatterplot trend

positive relationship: H

igh Y with

H

igh X &

Medium Y with Medium X

&

L

ow Y with

L

ow X

negative relationship: H

igh Y with

L

ow X &

Medium Y with Medium X

&

L

ow Y with

H

igh X

Description of the Moran scatterplot

Positive spatial autocorrelation

- high values tend to be 2002 population density surrounded by nearby high values - intermediate values tend to be surrounded by nearby intermediate values - low values tend to be surrounded by nearby low values MC = 0.49

GR = 0.58

Description of the Moran scatterplot

Negative spatial autocorrelation

- high values tend to be competition for space surrounded by nearby low values - intermediate values tend to be surrounded by nearby intermediate values - low values tend to be surrounded by nearby high values MC = -0.16

s MC = 0.075

GR = 1.04

Graphical portrayals of spatial autocorrelation latent in transformed

As

data

Constructing eigenfunctions for filtering spatial autocorrelation out of georeferenced variables:

Moran Coefficient = (n/

1

T

C1

)

x

Y

T (

I

–

11

T /n)

C

(

I

–

11

T /n)

Y

/

Y

T

(

I

–

11

T

/n)

Y

the eigenfunctions come from

(

I

–

11

T /n)

C

(

I

–

11

T /n)

Random effects model

Y

 f(

Xβ



ξ

,

ε

) ξ ε is a random observation effect (differences among individual observational units) is a time-varying residual error (links to change over time) The

composite error

term is the sum of the two.

Random effects model: normally distributed intercept term

• ξ σ 2 ~ N(0, ) and uncorrelated with covariates • supports inference beyond the nonrandom sample analyzed • simplest is where intercept is allowed to vary across areal units (repeated observations are individual time series) • The random effect variable is integrated out (with numerical methods) of the likelihood fcn • accounts for missing variables & within unit correlation ( commonality across time periods )

Spatial structuring of random effects

• • •

CAR

: conditional autoregressive model

ICAR

: improper conditional autoregressive model (spatial autocorrelation set to 1,and a spatially structured and a spatially unstructured variance component is estimated) ─should be specified as a convolution prior ( spatially structured & unstructured random effects )

SF

: spatial filter identified with a frequentist GLM

Definition of probability Point estimate Frequentist Bayesian

Long-run expected frequency in repeated (actual or hypothetical) experiments (Law of LN) Relative degree of belief in the state of the world Maximum likelihood estimate Mean, mode or median of the posterior probability distribution

Uncertainty intervals for parameters

“ confidence intervals ” based on the Likelihood Ratio Test (LRT) i.e., the expected probability distribution of the maximum likelihood estimate over many experiments “ credible intervals ” based on the posterior probability distribution

Uncertainty intervals of non parameters

Based on likelihood profile/LRT, or by resampling from the sampling distribution of the parameter

Model selection Difficulties

Discard terms that are not significantly different from a nested (null) model at a previously set confidence level Confidence intervals are confusing (range that will contain the true value in a proportion α of repeated experiments); rejection of model terms for “non significance” Calculated directly from the distribution of parameters Retain terms in models, on the argument that processes are not absent simply because they are not statistically significant Subjectivity; need to specify priors

Impact of sample size

•prior •distribution •likelihood •distribution As the sample size increases, a prior distri bution has less and less impact on results;

BUT

•effective •sample size •for spatially •autocorrelated •data

What is

BUGS

?

B

ayesian inference

U

sing

G

ibbs

S

ampling • is a piece of computer software for the Bayesian analysis of complex statistical models (MCMC) using Markov chain Monte Carlo methods. • It grew from a statistical research project at the MRC BIOSTATISTICAL UNIT in Cambridge, but now is developed jointly with the Imperial College School of Medicine at St Mary’s, London.

•BUGS • • •Classic BUGS •WinBUGS (Windows Version) GeoBUGS (spatial models) PKBUGS ( pharmokinetic modeling) • The

Classic BUGS

program uses text based model description and a command line interface, and versions are available for major computer platforms (e.g., Sparc, Dos). However, it is not being further developed.

What is

WinBUGS

?

• • WinBUGS and , a windows program with an option of a graphical user interface, the standard ‘point click’ windows interface, and on-line monitoring and convergence diagnostics. It also supports Batch-mode running (version 1.4).

GeoBUGS

spatial models and produces a range of maps as output. , an add-on to WinBUGS that fits • PKBUGS, an efficient and user-friendly interface for specifying complex population pharmacokinetic and pharmacodynamic (PK/PD) models within the WinBUGS software.

What is

GeoBUGS

?

• Available via http://www.mrc-bsu.cam.ac.uk/ bugs/winbugs/geobugs.shtml

• Bayesian inference is used to spatially smooth the standardized incidence ratios using Markov chain Monte Carlo (MCMC) methods.

GeoBUGS

implements models for data that are collected within discrete regions (not at the individual level), and smoothing is done based on Markov random field models for the neighborhood structure of the regions relative to each other.

What is

MCMC

?

MCMC

is used to simulate from some distribution

p

known only up to a constant factor, C:

p i

= Cq i where q i is known but C is unknown and too horrible to calculate.

MCMC

begins with

conditional

(marginal) distributions, and MCMC sampling outputs a sample of parameters drawn from their

posterior

(

joint

) distribution.

The geographic distribution of elevation across the island of Puerto Rico

From a USGS DEM containing 87,358,136 points. Darkness of gray scale is directly proportional to elevation.

SAS PROC MIXED summary results for a Semivario gram model Variance (nugget) Spatial correlation Residual b 0 b u 2 b uv b v b v 2 none -- -- 0.217

6.080*** -0.349*** -0.263*** -0.270*** -0.527*** spherical 0.0331

< 0.0001

0.169

6.080*** -0.349*** -0.263*** -0.270*** -0.527*** expo nential 0.2151

Gaussian 0.2210

power 0.2514

Bessel 0.2450

0.7643

0.5730

0.2702

< 0.0001

6.157*** 0.0253

6.201*** < 0.0001

6.157*** -0.463*** -0.486*** -0.463*** -0.254** -0.114

-0.255** -0.168

-0.254** -0.114

-0.529*** -0.561*** -0.529*** 0.6089

0.0057

6.202*** -0.486*** -0.257** -0.168

-0.569***

The average random effects term example MCMC chain from a WinBUGS run ICAR spatial filter (SF)

WinBUGS: geographic distributions of unstructured (left) and spatially structured (right) random effects

spatial filter (SF) WinBUGS: ICAR

Comparative parameter estimates for a LMM Param eter SAS semivariogram (Bessel) model estimate se SAS SF estimate se GeoBUGS-ICAR WinBUGS-SF (100 weeded replications) estimate se estimate se b 0 b u 2 b uv b v 2 var varure varssre 6.1906 0.287

-0.5048 0.122

-0.2229 0.123

-0.5314 0.125

0.0055 0.019

0.2856 0.091

0.7205 0.221

6.1101 0.055

-0.3881 0.031

-0.2939

-0.5193

0 0.0001

0.0282

0.030

0.032

-- -- -- 6.5175

0.168

6.1101

0.061

-0.7507 0.153 -0.3878 0.035

-0.2031 0.101 -0.2920 0.030

-0.5683 0.061 -0.5190 0.037

0.0049

0.0047

0.4854

0.007

0.007

0.093

0.0305

0.0318

0.0301

0.024

0.025

-- varure denotes the variance of the unstructured random effects varssre denotes the variance of the spatially structured random effects

binomial GLMM random effects

SAS SF WinBUGS SF WinBUGS ICAR

SF GLMM

statistic b 0 b elev  elev b

E

1 b

E

4 μˆ ξ σˆ 2 ξ σˆ ξ 2  SS P(S-W) MC ss GR ss MC ξˆ GR MC ξˆ ξˆ  SS GR r ξˆ , elev ξˆ  SS r ξˆ ,

E

1 r ξˆ ,

E

4 SAS NLMIXED (SF) estimate se -1.2867

-0.0111

3.0646

3.2182

0.0015

0.7045

0.9787

<0.0001

0.967

0.158

0.119

1.045

0.356

0.739

0.001

-0.001

0.001

0.2624

0.0013

1.1559

1.3433

WinBUGS (100 weeded replications) SF ICAR estimate se estimate se -1.3114

-0.0110

3.0600

3.0116

0.0054

0.7144

0.9783

< 0.0001

0.975

0.154

0.132

1.000

0.357

0.739

0.001

-0.009

0.022

0.2852

0.0014

1.3632

1.4256

-1.5340

-0.0100

*** *** 0.0066

0.3727

0.9583

< 0.0001

0.787

0.177

0.036

1.129

0.388

0.696

0.011

*** *** 0.2419

0.0013

Graphical diagnostics of residuals for the GLMM estimated with SAS

Scatterplot of the SAS and mean WinBUGS estimated spatially structured random effects terms

Individual GLMM estimation results for each Puerto Rican sugar cane crop year

Crop year Individ ual SF eigen vector #s Raw per centages MC GR

1965/ 66 1966/ 67 1967/ 68

0.484 0.458

1,4,6,24 0.490 0.454

0.474 0.434

# 0s 3 4 6 -a Point-in-time estimation b elev  elev μˆ ξ σˆ 2 ξ P(S W) 1.1959 0.0064 0.0020 1.0135 0.0055

1.3786 0.0060 0.0021 1.1138 0.0027

1.7354 0.0050 0.0018 1.0856 0.0017

Space-time data: preliminaries

random effect (re) re + ss: 1965/66 re + ss: 1966/67 re + ss: 1967/68 Random effects term is constant across time; spatial structuring changes over time

Space-time GLMM: Puerto Rican sugar cane crop years 1965/66-1967/68 when all fixed effects are year-specific statistic b 0 b elev  elev b

E

1 b

E

4 b

E

6 b

E

2 4 pseudo-R 2 MC ξˆ  SS GR ξˆ  SS MC residuals GR residuals crop year 1965/66 estimate se -1.2291

-0.0065

0.2336

0.0009

4.5040

4.9713

1.1684

1.2432

-4.5091

-4.0290

0.9950

0.449

1.1620

1.0994

0.580

0.019

0.916

crop year 1966/7 estimate se -1.3122

-0.0064

0.2336

0.0009

4.5785

5.3372

1.1684

1.2432

-4.8209

-3.9657

0.9976

0.467

1.1620

1.0994

0.562

0.042

0.839

crop year 1967/68 estimate se -1.4520

-0.0065

0.2336

0.0009

4.9226

5.8053

1.1684

1.2432

-5.0203

-4.0285

0.9929

0.493

1.1621

1.0995

0.537

0.009

0.808

Discussion & Implications

1. All three common specifications of spatial structuring —semivariogram, spatial autoregressive and SF models —for a random effect term in mixed statistical models perform in an equivalent fashion. 2. Matching Bayesian model priors with their implicit frequentist counterparts yields estimation results from both approaches that are essentially the same.

3. making use of spatially structured random effects tends to furnish an alternative to quasi likelihood estimation techniques for GLMMs

4. Semivariogram models offer a geostatistical theoretical basis and have been implemented in SAS for LMMs.

• A spatial statistics practitioner with the necessary computer programming skills can employ WinBUGS in order to utilize them with GLMMs.

5. Spatial autoregressive modeling offers a theoretical basis for spatial structuring, and is available in GeoBUGS.

• This would be very difficult to trick SAS into doing. 6. Spatial filtering, which can be derived from spatial • • • autoregressive model specifications, tends to be more exploratory in nature (being akin to principal components analysis) can be implemented in either SAS or WinBUGS for either LMMs or GLMMs, and can be easily extended to space-time datasets with either of these software packages.

7. Illustrative Puerto Rico sugar cane examples tend to have a random effect term that virtually equates to the corresponding LMM/GLMM residual variate.

• This is not always the case, as is highlighted by the extension of a GLMM specification to a space-time sugar cane dataset.

8. All of the estimated random effects terms for the various Puerto Rico examples tend to be non-normal.

9. once a random effect term has been estimated with a frequentist approach, using it when calculating a deviance statistic allows its number of degrees of freedom to be approximated for GLMMs.

• Although n values are estimated, because they are correlated, the resulting number of degrees of freedom is less than n.

• This particular finding should help spatial statistics practitioners better understand the cost of employing a statistical mixed model.

A df aside: future research

• Spiegelhalter et al. (2002) address the df problem for complex hierarchical models in which the number of parameters is not clearly defined because, for instance, of the presence of random effects.

• An information-theoretic argument is used to approximate the effective number of parameters in a model, equivalent to the trace of the product of the Fisher information and the posterior covariance matrices.

– this particular approximation is equivalent to the trace of the ‘hat’ matrix for linear models with a normally distributed error term.

k dfs for random effects k



n



(p



1)



deviance 1

binomial

k



n



(p



1)



deviance deviance

Poisson neative binomial