Lecture 4 bayesian analysis of mapped data
Download
1 / 54

Lecture 4: Bayesian analysis of mapped data - PowerPoint PPT Presentation


  • 95 Views
  • Uploaded on

Lecture #4: Bayesian analysis of mapped data. Spatial statistics in practice Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden. Topics for today’s lecture. Frequentist versus Bayesian perspectives. Implementing random effects models in GeoBUGS.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Lecture 4: Bayesian analysis of mapped data' - gur


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Lecture 4 bayesian analysis of mapped data

Lecture #4:Bayesian analysis of mapped data

Spatial statistics in practice

Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden


Topics for today s lecture
Topics for today’s lecture

  • Frequentist versus Bayesian perspectives.

  • Implementing random effects models in GeoBUGS.

  • Spatially structured and unstructured random effects: the CAR, the ICAR, and the spatial filter specifications


Lmm glmm move in the direction of bayesian modeling
LMM & GLMM move in the direction of Bayesian modeling

  • Fixed effects are in keeping with a frequentist viewpoint―individual unknown parameters

  • Random effects are distributions of parameters, and are in keeping with a Bayesian viewpoint

  • The model specifications tend to be the same, with estimation methods tending to differ


Sampling distribution based inference the reported p values
Sampling distribution-based inference: the reported p values

  • Frequentism: assigns probabilities to random events only according to their relative frequencies of occurrence, rejecting degree-of-belief interpretations of mathematical probability theory.

  • The frequentist approach considers  to be an unknown constant. Allowing for the possibility that  can take on a range of values, frequentist inference is based on a hypothetical repeated sampling principle that obtains desirable and (often) physically interpretable properties (e.g., CIs).

  • Frequentist statistics typically is limited to posing questions in terms of a null hypothesis (i.e., H0) that a parameter takes on a single value.


What is bayesian analysis
What is Bayesian analysis?

  • Bayesianism: contends that mathematical probability theory pertains to degree of plausibility/belief; when used with Bayes theorem (which casts analysis in conditional terms), it becomes Bayesian inference.

  • The Bayesian approach considers  as the realized value of a random variable Θ with a probability density (mass) function called the prior distribution (a marginal probability distribution that is interpreted as a description of what is known about a variable in the absence of empirical/theoretical evidence).

  • Bayesian statistics furnishes a generic tool for inferring model parameter probability distribution functions from data: a parameter has a distribution, not a single value.


  • The Bayesian interpretation of probability allows (proper prior) probabilities to be assigned subjectively to random events, in accordance with a researcher’s beliefs.

  • Proper prior: a probability distribution that is normalized to one, and asymptotically dominated by the likelihood function

    • improper prior: a handy analytic form (such as a constant or logarithmic distribution) that, when integrated over a parameter space, tends to 1 and so is not normalizable.

  • Noninformative prior: because no prior information is available, the selected prior contains vague/general information about a parameter, having minimal influence on inferences.


  • conjugate prior: a prior probability distribution naturally connected with the likelihood function that, when combined with the likelihood and normalized, produces a posterior probability distribution that is of the same type as the prior, making the analytical mathematics of integration easy. (NOTE: MCMC techniques relaxes the need for this type of specification)

  • The hierarchical Bayes model is called “hierarchical” because it has two levels. At the higher level, hyper-parameter distributions are described by multivariate priors. Such distributions are characterized by vectors of means and covariance matrices; spatial autocorrelation is captured here. At the lower level, individuals’ behavior is described by probabilities of achieving some outcome that are governed by a particular model specification.


  • posterior distribution: a conditional probability distribution for a parameter taking empirical evidence into account  computed with Bayes’ theorem as a normalized product of the likelihood function and the prior distributions  that supports Bayesian inference about the parameter

  • Bayesian statistics assumes that the probability distribution in question is known, and hence involves integration to get the normalizing constant. This integration might be tricky, and in many cases there is no analytical solution. With the advent of computers and various integration techniques, this problem can partially be overcome. In many Bayesian statistics applications a prior is tabulated and then sophisticated numerical integration techniques (e.g., MCMC) are used to derive posterior distributions.


The controversy
The controversy

  • A frequentist generally has no objection to the use of the Bayes formula, but when prior probabilities are lacking s/he deplores the Bayesians’ tendency to make them up out of thin air.

  • Does a parameter have one or many values?

  • For many (but not all) of the simpler problems where frequentist methodology seems to give satisfying answers, the Bayesian approach yields basically the same answers, provided noninformative proper priors are used.


Impact of sample size
Impact of sample size

prior

distribution

As the sample size increases, a prior distri-bution has less and less impact on results; BUT

likelihood

distribution

effective

sample size

for spatially

autocorrelated

data


What is bugs
What is BUGS?

Bayesian inference Using Gibbs Sampling

  • is a piece of computer software for the Bayesian analysis of complex statistical models using Markov chain Monte Carlo (MCMC) 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.


Classic BUGS

WinBUGS (Windows Version)

BUGS

GeoBUGS (spatial models)

PKBUGS (pharmokinetic modelling)

  • 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
What is WinBUGS?

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

  • GeoBUGS, an add-on to WinBUGS that fits spatial models and produces a range of maps as output.

  • PKBUGS, an efficient and user-friendly interface for specifying complex population pharmacokinetic and pharmacodynamic (PK/PD) models within the WinBUGS software.


What is geobugs
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.


Choice of priors
Choice of Priors

  • Can effect

    • model convergence

    • execution speed

    • computational overflow errors

  • Consider

    • Conjugate prior

    • Sensible prior for the range of a parameter

      • e.g. Poisson mean must be positive; hence, a normal distribution is not a good prior.

    • Truncated Prior

      • useful if it is known that sensible or useful parameter values are within a particular range


Parameter initialization
Parameter Initialization

  • Initialize important parameters.

    • Parameters distributed according to a non-informative prior must be initialized (otherwise an overflow error is very likely).

GeoBUGs

  • Regions are numbered sequentially from 1 to n

  • Each region is defined as a polygon in a map file

  • Each region is associated with a unique index

  • BUGs can import map files from Arcinfo, Epimap, SPLUS.


Review what is mcmc

Review: What is MCMC?

MCMC is used to simulate from some distribution p known only up to a constant factor, C:

pi = Cqi

where qi 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 joint (posterior) distribution.


Geobugs spatial poisson with a map
GeoBugs:spatial Poisson with a map

model;

{

for (i in 1:M){

for (j in 1:N){

#Poisson with a constant but unknown mean

y[i,j] ~ dpois(l)

#GeoBUGS numbers regions in a sequence

#This statement turns a 2-D array into a 1-D array for GeoBUGS

mapy[j + (i-1)N] <- y[I,j]

}

}

#priors

#Noninformative conjugate priors for Poisson mean.

l~dgamma(0.001,0.001)

}


Displaying a map with geobugs
Displaying a map with GeoBugs

  • Compile model, load data, load initial conditions

  • Set sample monitor on desired variables

  • Set trace

  • Set sample monitor on map variable OR set summary monitor on map variable

  • Run chain

  • Activate map tool, load appropriate map

  • Set cut points, colour spectrum as desired.


Poisson mixture on a spatial lattice
Poisson mixture on a spatial lattice

Poisson mean varies

x[i,j] in (0,1)

x[i,j] independent binomials

Map the hidden lattice

and the Poisson means

E(mu) = a/b; var(mu) = a/b2

a b


Introducing a covariate
Introducing a covariate

Poisson mean varies with location

Log link function

Normal prior for covariate parameter


Adding a random effect
Adding a random effect

Random effect


Geobugs example
GeoBUGS Example

NOTE: Possible mappings for the Scottish lip cancer data include

b – area-specific random effects term

mu – area-specific means

RR – area-specific relative risks

O – observed values

E – expected values


The scottish lip cancer example
The Scottish lip cancer example

Open a window in WinBUGS and insert the following items:

- model

- data

- initial

values

#MODEL

model

{

for (i in 1 : N) {

O[i] ~ dpois(mu[i])

log(mu[i]) <- log(E[i]) + alpha0 + alpha1 * X[i]/10 + b[i]

# Area-specific relative risk (for maps)

RR[i] <- exp(alpha0 + alpha1 * X[i]/10 + b[i])}

# CAR prior distribution for random effects:

b[1:N] ~ car.normal(adj[], weights[], num[], tau)

for(k in 1:sumNumNeigh) {

weights[k] <- 1

}

# Other priors:

alpha0 ~ dflat()

alpha1 ~ dnorm(0.0, 1.0E-5)

tau ~ dgamma(0.5, 0.0005) # prior on precision

sigma <- sqrt(1 / tau) # standard deviation

}


#DATA

list(N = 56,

O = c( 9, 39, 11, 9, 15, 8, 26, 7, 6, 20,

13, 5, 3, 8, 17, 9, 2, 7, 9, 7,

16, 31, 11, 7, 19, 15, 7, 10, 16, 11,

5, 3, 7, 8, 11, 9, 11, 8, 6, 4,

10, 8, 2, 6, 19, 3, 2, 3, 28, 6,

1, 1, 1, 1, 0, 0),

E = c( 1.4, 8.7, 3.0, 2.5, 4.3, 2.4, 8.1, 2.3, 2.0, 6.6,

4.4, 1.8, 1.1, 3.3, 7.8, 4.6, 1.1, 4.2, 5.5, 4.4,

10.5, 22.7, 8.8, 5.6, 15.5, 12.5, 6.0, 9.0, 14.4, 10.2,

4.8, 2.9, 7.0, 8.5, 12.3, 10.1, 12.7, 9.4, 7.2, 5.3,

18.8, 15.8, 4.3, 14.6, 50.7, 8.2, 5.6, 9.3, 88.7, 19.6,

3.4, 3.6, 5.7, 7.0, 4.2, 1.8),

X = c(16, 16, 10, 24, 10, 24, 10, 7, 7, 16,

7, 16, 10, 24, 7, 16, 10, 7, 7, 10,

7, 16, 10, 7, 1, 1, 7, 7, 10, 10,

7, 24, 10, 7, 7, 0, 10, 1, 16, 0,

1, 16, 16, 0, 1, 7, 1, 1, 0, 1,

1, 0, 1, 1, 16, 10),


43, 29, 25,

56, 32, 31, 24,

45, 33, 18, 4,

50, 43, 34, 26, 25, 23, 21, 17, 16, 15, 9,

55, 45, 44, 42, 38, 24,

47, 46, 35, 32, 27, 24, 14,

31, 27, 14,

55, 45, 28, 18,

54, 52, 51, 43, 42, 40, 39, 29, 23,

46, 37, 31, 14,

41, 37,

46, 41, 36, 35,

54, 51, 49, 44, 42, 30,

40, 34, 23,

52, 49, 39, 34,

53, 49, 46, 37, 36,

51, 43, 38, 34, 30,

42, 34, 29, 26,

49, 48, 38, 30, 24,

55, 33, 30, 28,

53, 47, 41, 37, 35, 31,

53, 49, 48, 46, 31, 24,

49, 47, 44, 24,

54, 53, 52, 48, 47, 44, 41, 40, 38,

29, 21,

54, 42, 38, 34,

54, 49, 40, 34,

49, 47, 46, 41,

52, 51, 49, 38, 34,

56, 45, 33, 30, 24, 18,

55, 27, 24, 20, 18

),

sumNumNeigh = 234)

num = c(3, 2, 1, 3, 3, 0, 5, 0, 5, 4,

0, 2, 3, 3, 2, 6, 6, 6, 5, 3,

3, 2, 4, 8, 3, 3, 4, 4, 11, 6,

7, 3, 4, 9, 4, 2, 4, 6, 3, 4,

5, 5, 4, 5, 4, 6, 6, 4, 9, 2,

4, 4, 4, 5, 6, 5

),

adj = c(

19, 9, 5,

10, 7,

12,

28, 20, 18,

19, 12, 1,

17, 16, 13, 10, 2,

29, 23, 19, 17, 1,

22, 16, 7, 2,

5, 3,

19, 17, 7,

35, 32, 31,

29, 25,

29, 22, 21, 17, 10, 7,

29, 19, 16, 13, 9, 7,

56, 55, 33, 28, 20, 4,

17, 13, 9, 5, 1,

56, 18, 4,

50, 29, 16,

16, 10,

39, 34, 29, 9,

56, 55, 48, 47, 44, 31, 30, 27,

29, 26, 15,


#INITIAL VALUES

list(tau = 1, alpha0 = 0, alpha1 = 0,

b=c(0,0,0,0,0,NA,0,NA,0,0,

NA,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))


Geobugs execution instructions

GeoBUGS Execution Instructions

Step 1: open a new window in WinBUGS (this will be referred to as the user window)

Step 2: enter the model syntax, data, and initial values, using the WinBUGS format, in the user window

Step 3: select “Specification” in the “Model” pull down window

Step 4: highlight “model” in the user window, appearing at the beginning of the model syntax, and click once on the “check model” button in the “Specification Tool” window

NOTE: feedback from the program appears in the lower left-hand corner of the WinBUGS program window,

and should be monitored


Step 5: highlight “list” in the user window, appearing at the beginning of the data, and click once on the “load data” button in the “Specification Tool” window

Step 6: insert the number of chains to be run (the default number is 1) in the “Specification Tool” window

Step 7: click once on the “compile” button in the “Specification Tool” window

Step 8: highlight “list” in the user window, appearing at the beginning of the initial values, and click once on the “load inits” button in the “Specification Tool” window (one set is needed for each chain to be run; clicking the “gen inits” button can be dangerous for sound analysis)


Step 9: close the “Specification Tool” window

Step 10: select “Samples” in the “Inference” pull down menu

Step 11: in the “Sample Monitor Tool” window: type in the name (as it appears in the model syntax) of the 1st parameter to be monitored, and click once on the “set” button; type in the name of the 2nd parameter to be monitored, and click once on the “set” button; …; type in the name of the pth parameter to be monitored, and click once on the “set” button

Step 12: close the “Sample Monitor Tool” window


Step 13: select “Update” in the “Model” pull down menu

Step 14: the default value in the “updates” box is 1000this value most likely needs to be substantially increased (say, to 50,000; once any desired changes are made, click once on the “updates” button

Step 15: once the number appearing in the “iteration” box equals the number in the “updates” box, close the “Update Tool” window

Step 16: select “Samples” in the “Inference” pull down menu

Step 17: click on the down arrow to the right of the “node” box, and select the parameter to be monitored


Step 18: click once of the “history” button to view the time-series plot for all iterations; click once on the “trace” plot to view the time-series plot for the last iterations; click once on the “auto cor” button to view the time-series correlogram (you may wish to enlarge the graphic in this window); click once on the “stats” button to view a parameter estimate’s statistics

Step 19: close the “Sample Monitor Tool” window

Step 20: select “Mapping Tool” from the “Map” pull down menu

Step 21: select the appropriate map from the list appearing for the “Map” box when the down arrow to its right is clicked once (hint: your map that you imported from an Arc shapefile should appear here)

Step 22: type the name of the variable (exactly as it appears in the model syntax) to be mapped in the “variable” box, and click once on the “plot” box


Puerto rico binomial with spatial autocorrelation example
Puerto Rico Binomial with Spatial Autocorrelation Example

#MODEL

model

{

for (i in 1 : N) {

U[i] ~ dbin(p[i], T[i])

p[i] <- exp(alpha0 + b[i])/(1+exp(alpha0 + b[i]))

}

# CAR prior distribution for random effects:

b[1:N] ~ car.normal(adj[], weights[], num[], tau)

for(k in 1:sumNumNeigh) {

weights[k] <- 1

}

# Other priors:

alpha0 ~ dflat()

tau ~ dgamma(0.5, 0.0005) # prior on precision

sigma <- sqrt(1 / tau) # standard deviation

}


#DATA

list(N = 76,

U = c(42527, 11048, 46236, 35859, 21330,

25584, 3792, 32126, 32389, 18664,

41997, 2839, 11787, 95880, 37713,

27605, 21499, 29709, 17200, 14688,

23829, 178792, 24196, 14767, 50242,

23848, 28462, 34650, 434374, 35130,

38583, 17412, 63929, 94085, 75728,

23852, 36971, 59572, 23364, 37238,

40919, 11062, 42042, 64685, 27305,

23077, 25387, 91593, 18346, 22105,

27850, 224044, 40875, 139445, 30886,

42467, 185703, 30071, 43707, 16671,

14262, 40457, 29802, 16800, 35270,

33421, 39958, 10176, 20682, 40395,

21087, 99850, 35476, 36201, 16472,

58848

),

T = c(44444, 17318, 50531, 36452, 26261,

34415, 11061, 34485, 32537, 19817,

45409, 6449, 12741, 98434, 39697,

29965, 23753, 29709, 23844, 20152,

26719, 186475, 25450, 14767, 52362,

25935, 31113, 37105, 434374, 40997,

44204, 21665, 63929, 94085, 75728,

35336, 37910, 61929, 27913, 39246,

46384, 19143, 42042, 64685, 29032,

26493, 28348, 100131, 19117, 22322,

28909, 224044, 46911, 140502, 35244,

43335, 186076, 30071, 47370, 18004,

19811, 42753, 37597, 20002, 36867,

34017, 40712, 12367, 21888, 44301,

23072, 100053, 36743, 38925, 16614,

59035

),

num = c(4, 3, 5, 3, 4, 4, 3, 3, 5, 5,

3, 4, 5, 5, 3, 3, 6, 4, 7, 6,

4, 5, 4, 6, 7, 3, 2, 4, 6, 2,

7, 6, 8, 7, 5, 6, 7, 5, 6, 5,

4, 4, 8, 5, 5, 7, 6, 5, 6, 7,

6, 7, 6, 3, 5, 3, 7, 6, 5, 6,

7, 4, 3, 6, 5, 3, 3, 5, 6, 4,

5, 4, 2, 3, 2, 3

),


9, 11, 12, 19, 32, 35,

1, 6, 7, 14, 33, 36, 44,

20, 27, 41,

26, 41,

13, 17, 37, 38,

4, 9, 10, 31, 32, 43,

21, 36,

3, 10, 23, 29, 34, 40, 43,

9, 24, 29, 35, 43, 50,

5, 6, 25, 34, 44, 49, 53, 60,

3, 5, 31, 33, 40, 49, 59,

19, 24, 32, 45, 50,

14, 21, 25, 30, 44, 47,

13, 28, 38, 39, 51, 52, 61,

17, 28, 37, 48, 52,

13, 18, 19, 37, 45, 51,

31, 34, 43, 59, 64,

20, 26, 27, 42,

20, 41, 46, 54,

29, 31, 32, 40, 50, 56, 57, 64,

25, 33, 36, 47, 53,

19, 35, 39, 50, 51,

20, 22, 42, 48, 52, 54, 68,

36, 44, 53, 58, 62, 63,

17, 22, 38, 46, 52,

33, 34, 59, 60, 66, 67,

32, 35, 43, 45, 51, 55, 57,

37, 39, 45, 50, 55, 61,

37, 38, 46, 48, 61, 68, 69,

33, 44, 47, 58, 60, 65,

42, 46, 68,

50, 51, 57, 61, 71,

43, 57, 64,

43, 50, 55, 56, 64, 71, 74,

47, 53, 62, 63, 65, 72,

34, 40, 49, 64, 67,

33, 49, 53, 65, 66, 76,

37, 51, 52, 55, 69, 70, 71,

47, 58, 63, 72,

47, 58, 62,

40, 43, 56, 57, 59, 74,

53, 58, 60, 72, 76,

49, 60, 67,

49, 59, 66,

46, 52, 54, 69, 73,

52, 61, 68, 70, 73, 75,

61, 69, 71, 75,

55, 57, 61, 70, 74,

58, 62, 65, 76,

68, 69,

57, 64, 71,

69, 70,

60, 65, 72

),

sumNumNeigh = 366)

#INITIAL VALUES

list(tau = 1, alpha0 = 0,

b=c(0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0)

)

queen’s

connectivity

definition

adj = c(

2, 7, 14, 25,

1, 14, 21,

5, 8, 23, 31, 34,

9, 10, 29,

3, 6, 33, 34,

5, 7, 25, 33,

1, 6, 25,

3, 10, 23,

4, 11, 24, 29, 32,

4, 8, 23, 29, 31,

9, 12, 24,

11, 16, 19, 24,

17, 18, 28, 37, 39,

1, 2, 21, 25, 36,

17, 20, 22,

12, 18, 19,

13, 15, 22, 28, 38, 48,

13, 16, 19, 39,

12, 16, 18, 24, 35, 39, 45,

15, 22, 26, 41, 42, 46,

2, 14, 30, 36,

15, 17, 20, 46, 48,

3, 8, 10, 31,

9, 11, 12, 19, 32, 35,

1, 6, 7, 14, 33, 36, 44,

20, 27, 41,

26, 41,

13, 17, 37, 38,

4, 9, 10, 31, 32, 43,

21, 36,

3, 10, 23, 29, 34, 40, 43,

9, 24, 29, 35, 43, 50,

5, 6, 25, 34, 44, 49, 53, 60,

3, 5, 31, 33, 40, 49, 59,

19, 24, 32, 45, 50,

14, 21, 25, 30, 44, 47,

13, 28, 38, 39, 51, 52, 61,

17, 28, 37, 48, 52,

13, 18, 19, 37, 45, 51,

31, 34, 43, 59, 64,

20, 26, 27, 42,

20, 41, 46, 54,

29, 31, 32, 40, 50, 56, 57, 64,

25, 33, 36, 47, 53,

19, 35, 39, 50, 51,

20, 22, 42, 48, 52, 54, 68,

36, 44, 53, 58, 62, 63,

17, 22, 38, 46, 52,

33, 34, 59, 60, 66, 67,

32, 35, 43, 45, 51, 55, 57,

37, 39, 45, 50, 55, 61,

37, 38, 46, 48, 61, 68, 69,

33, 44, 47, 58, 60, 65,

42, 46, 68,

50, 51, 57, 61, 71,

43, 57, 64,

43, 50, 55, 56, 64, 71, 74,

47, 53, 62, 63, 65, 72,

34, 40, 49, 64, 67,

33, 49, 53, 65, 66, 76,

37, 51, 52, 55, 69, 70, 71,

47, 58, 63, 72,

47, 58, 62,

40, 43, 56, 57, 59, 74,

53, 58, 60, 72, 76,

49, 60, 67,

49, 59, 66,

46, 52, 54, 69, 73,

52, 61, 68, 70, 73, 75,

61, 69, 71, 75,

55, 57, 61, 70, 74,

58, 62, 65, 76,

68, 69,

57, 64, 71,

69, 70,

60, 65, 72, 2, 7, 14, 25,

1, 14, 21,

5, 8, 23, 31, 34,

9, 10, 29,

3, 6, 33, 34,

5, 7, 25, 33,

1, 6, 25,

3, 10, 23,

4, 11, 24, 29, 32,

4, 8, 23, 29, 31,

9, 12, 24,

11, 16, 19, 24,

17, 18, 28, 37, 39,

1, 2, 21, 25, 36,

17, 20, 22,

12, 18, 19,

13, 15, 22, 28, 38, 48,

13, 16, 19, 39,

12, 16, 18, 24, 35, 39, 45,

15, 22, 26, 41, 42, 46,

2, 14, 30, 36,

15, 17, 20, 46, 48,

3, 8, 10, 31,


Model {

for (i in 1:N) {m[i] <- 1/num[i]}

cumsum[1] <- 0

for(i in 2:(N+1)) {cumsum[i] <- sum(num[1:(i-1)])}

for(k in 1 : sumNumNeigh)

{for(i in 1:N){pick[k,i] <- step(k - cumsum[i] - epsilon)*step(cumsum[i+1] - k)

# pick[k,i] = 1 if cumsum[i] < k <= cumsum[i=1];

otherwise, pick[k,i] = 0}

C[k] <- sqrt(num[adj[k]]/inprod(num[], pick[k,])) # weight for each pair of neighbours}

epsilon <- 0.0001

for (i in 1 : N) {

U[i] ~ dbin(p[i], T[i])

p[i] <- exp(S[i])/(1+exp(S[i]))

theta[i] <- alpha}

# proper CAR prior distribution for random effects:

S[1:N] ~ car.proper(theta[],C[],adj[],num[],m[],prec,rho)

for(k in 1:sumNumNeigh) {weights[k] <- 1}

# Other priors:

alpha ~ dnorm(0,0.0001)

prec ~ dgamma(0.5,0.0005) # prior on precision

sigma <- sqrt(1/prec) # standard deviation

rho.min <- min.bound(C[],adj[],num[],m[])

rho.max <- max.bound(C[],adj[],num[],m[])

rho ~ dunif(rho.min,rho.max)

}


node mean sd MC error 2.5% median 97.5% start samplerhoCAR 0.9378 0.05255 5.995E-4 0.8047 0.9514 0.9973 50001 10000

MCMC iteration time series plot

MCMC iteration correlogram


Geobugs demonstration urban population in puerto rico no random effect
GeoBUGS demonstration: % urban population in Puerto Rico – no random effect

#MODEL

model

{

for (i in 1 : N) {

U[i] ~ dbin(p[i], T[i])

p[i] <- exp(alpha)/(1+exp(alpha ))

}

# priors:

alpha ~ dflat()

}


#DATA

list(N = 73,

U = c(11062, 42042, 64685, 27305, 23077, 25387, 91593, 18346, 32281, 27850, 254115,

40875, 139445, 30886, 185703, 43707, 16671, 14262, 40457, 29802, 16800, 35270,

33421, 39958, 20682, 40395, 21087, 99850, 35476, 36201, 16472, 58848, 42527, 11048,

46236, 35859, 21330, 25584, 3792, 32126, 74856, 18664, 41997, 2839, 11787, 95880,

37713, 27605, 21499, 29709, 17200, 14688, 23829, 178792, 24196, 14767, 50242,

23848, 28462, 34650, 434374, 35130, 38583, 17412, 63929, 94085, 75728, 23852,

36971, 59572, 23364, 37238, 40919

),

T = c(19143, 42042, 64685, 29032, 26493, 28348, 100131, 19117, 34689, 28909, 254115,

46911, 140502, 35244, 186076, 47370, 18004, 19811, 42753, 37597, 20002, 36867,

34017, 40712, 21888, 44301, 23072, 100053, 36743, 38925, 16614, 59035, 44444,

17318, 50531, 36452, 26261, 34415, 11061, 34485, 75872, 19817, 45409, 6449, 12741,

98434, 39697, 29965, 23753, 29709, 23844, 20152, 26719, 186475, 25450, 14767,

52362, 25935, 31113, 37105, 434374, 40997, 44204, 21665, 63929, 94085, 75728,

35336, 37910, 61929, 27913, 39246, 46384

),

)

#INITIAL VALUES

list(alpha=-3)


The German labor market revisited:

frequentist and Bayesian

random effects models


Frequentist random intercept term
Frequentist random intercept term

The spatial filter contains 27 (of 98) eigenvectors, with R2 = 0.4542, P(S-Wresiduals) < 0.0001.


Mcmc results from geobugs
MCMC results from GeoBUGS

1000 weeded

the

proper

CAR

model

25,000

burn-in


Bayesian model estimation
Bayesian model estimation

computer

execution

time

required

to

estimate

is

prohibitive


The icar alternative
The ICAR alternative

  • #MODEL

  • model

  • {for (i in 1 : N) {

  • W[i] ~ dpois(mu[i])

  • log(mu[i]) <- alpha + S[i] + U[i] + log(A[i])

  • U[i] ~ dnorm(0,tau.U) }

  • # improper CAR prior distribution for random effects:

  • S[1:N] ~ car.normal(adj[],weights[],num[],tau.S)

  • for(k in 1:sumNumNeigh) {weights[k] <- 1}

  • # Other priors:

  • alpha ~ dflat()

  • tau.S ~ dgamma(0.5, 0.0005)

  • sigma2.S <- 1/tau.S

  • tau.U ~ dgamma(0.5, 0.0005)

  • sigma2.U <- 1/tau.U }

offset


attribute

variables

  • #DATA

  • list(N = 439,

  • W = c(

  • 25335, 74162, 64273, 25080, 39848, 59407, 49685, 60880, 101563, 39206, 83512,

  • 32999, 26534, 45267, 35232, 35788, 41961, 36129

  • ),

  • A = c(

  • 23.51, 43.81, 87.14, 27.53, 558.57, 512.17, 855.22, 578.86, 261.45, 447.18, 877.26,

  • 369.50, 382.46, 452.02, 350.84, 219.97

  • ),

  • num = c(1, 2, 4, 3, 4, 7, 2, 4, 4, 5,

  • 8, 7, 8, 4, 7, 7, 6, 8, 6

  • ),

  • adj = c(

  • 12,

  • 11, 10,

  • 359, 15, 8, 6,

  • 438, 400, 390, 375, 372, 368

  • ),

  • sumNumNeigh = 2314)

number of neighbors

This can be

generated by

GeoBUGS; the

queen’s definition

of adjacency is used

lists of neighbors


  • #INITIAL VALUES

  • list(tau.U=1, tau.S=1, alpha=0,

  • S=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

  • 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),

  • U=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

  • 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)

  • )


Intercept
intercept

Convergence has not been achieved (10,000 burn in; + 25,000/100)


Random effects spatially structured and unstructured
Random effects: spatially structured and unstructured

SA:

2.67/(30.39+2.67)

= 0.08

(8% of variance)



Bugs r functions for running winbugs from r
bugs.R: functions for running WinBUGS from R

The bugs() function takes data and starting values as input, calls a Bugs model, summarizes inferences and convergence in a table and graph, and saves the simulations in arrays for easy access in R.

1st: Download WinBugs1.4

2nd: Download OpenBUGS and extract the zip file into "c:/Program Files/“

3rd: Go into R and type install.packages("R2WinBUGS")


Some prediction comparisons
Some prediction comparisons

n = 56;

effective sample size = 24


What have we

learned today ?


ad