1 / 26

Lab #4: - using GeoBUGS

Lab #4: - using GeoBUGS. Spatial statistics in practice Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden. Converting Arc shapefiles to GeoBUGS maps. Open ArcMap Use Conversion Tools to covert the shapefile “ to coverage ”

destiny-clayton
Download Presentation

Lab #4: - using GeoBUGS

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. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Lab #4:- using GeoBUGS Spatial statistics in practice Center for Tropical Ecology and Biodiversity, Tunghai University & Fushan Botanical Garden

  2. Converting Arc shapefiles to GeoBUGS maps • Open ArcMap • Use Conversion Tools to covert the shapefile “to coverage” • Use Coverage tool “from coverage” to ungenerate (*** in version #9) • Choose “coverage generated” “feature type” • Choose “fixed” “numeric format” • Use the .txt extension for the output file

  3. not optional 4. Edit the ASCII output file by inserting • map:n number of areal units • Add Xscale: *** Yscale: *** • Consecutive integers 1, 2, …, n assignment to areal unit alphanumeric names (one areal unit per record line) • regions ArcMap number assignment to areal unit alphanumeric names (one areal unit per record line) END • Open a new window in WinBUGS & paste the edited output file from ArcMAP • Open the MAP pulldown menu and select “import ARCINFO”

  4. Puerto Rico example map n 1 -66.957328796537695 18.489151001160103 -66.952575683427725 18.477293014667296 -66.953041076730273 18.476812362760878 -66.962738037179861 18.473377228162935 -66.963653564340419 18.471937179979335 -66.958953857184653 18.465316772733544 -66.963493347318519 18.455747604562497 -66.966682434196983 18.449995041130034 -66.958244323698267 18.440557479468140 -66.949417114597665 18.438266754451419 -66.947563171530049 18.437808990450790 -66.947090148925568 18.436859131011197 -66.943260193118775 18.424034118939282 -66.941383362145373 18.421670913452495 -66.935333251734193 18.418399810947768 -66.934387207365916 18.417457580810105 -66.934379577223297 18.416503906299198 -66.934829711921992 18.415544510022020 -66.939422607169362 18.410736084262879 -66.941734314237607 18.409759521846183 -66.944046020465152 18.409736633099701 -66.945884704928901 18.408287048078517 -66.947257995249402 18.406368255524157 -66.948158263806107 18.403499603530204 -66.947212219437830 18.401601791345978 -66.942550658791703 18.399269103907603 -66.940231323262225 18.399291991813392 -66.931427002067423 18.398427963071374 -66.926795959469700 18.398473739723645 -66.924484253242156 18.398496627629434 -66.923553466637031 18.398506164467019 -66.922622680872607 18.398038864469498 -66.922142028966192 18.396614074889758 -66.922119140219706 18.393280029593246 -66.929054260583740 18.392734527359387 -66.930900573508723 18.391763687549659 -66.935035705596164 18.388387679831492 -66.935485839454174 18.387426376018659 -66.937316894616004 18.384548187187118 -66.937286376567585 18.381690978725711 -66.937278747265651 18.381214141890606 -66.936798095359237 18.379312515475760 -66.936775207453451 18.376930237154500 -66.945564270044400 18.375411987575255 -66.948814392178988 18.376808166642238 -66.959526062395611 18.381938933980852 -66.973922729160265 18.386077881139602 -66.979484558363112 18.386493682718580 -66.984558105516896 18.385011672954022 -66.985023498819459 18.385482788022856 -66.992004394154364 18.388744354530690 -66.997589111262997 18.391065597595823 -67.002243041766505 18.392923354894052 -67.005958557203655 18.394309997964143 -67.027252196982005 18.393598556101406 -67.050949096292442 18.443834304579834 -67.052383422709767 18.447631835643250 -67.053825378429025 18.452852248750748 -67.056228637961112 18.460924148553378 -67.056312561123022 18.468544006121654 -67.057701110888075 18.468528747517794 -67.058631896652500 18.468517303144555 -67.059181213116986 18.476133346482211 -67.062973022414127 18.484186172609672 -67.073211670026254 18.487878799300344 -67.079238891690949 18.487810134742283 -67.079696655691578 18.487802505440353 -67.082496643127470 18.489200592042991 -67.082984924335818 18.491100310922182 -67.082565307685542 18.495391845800910 -67.085388183867906 18.499166488958537 -67.091926574646749 18.503376007370314 -67.095703125340037 18.509048461871473 -67.097595214917291 18.512357711726541 -67.108276367085494 18.514137268469817 -67.101333618260227 18.515171051071331 -67.097625732125010 18.515689849487398 -67.090682983299743 18.516252517860771 -67.053573608102568 18.515251160002631 -67.045684814068323 18.514390945491229 -67.043823242539489 18.514408111630743 -67.036849975665803 18.512582779075888 -67.029884338094064 18.511707305960627 -67.027099609262024 18.511262893868896 -67.021530151597958 18.510374068844737 -67.013160705657299 18.508560180663125 -67.005241394415336 18.504835128388386 -66.977287292505196 18.491323470735875 -66.957328796537695 18.489151001160103 END . . . END regions 1 Isabela 2 Aguadilla 3 Arecibo 4 Vega_Baja 5 Hatillo 6 Camuy 7 Quebradillas 8 Barceloneta 9 Vega_Alta 10 Manati 11 Dorado . . . 64 Juana_Diaz 65 Sabana_Grande 66 Guayanilla 67 Penuelas 68 Yabucoa 69 Patillas 70 Guayama 71 Salinas 72 Lajas 73 Maunabo 74 Santa_Isabel 75 Arroyo 76 Guanica END map:76 Xscale: 69.1 Yscale: 65.7 1 Isabela 2 Aguadilla 3 Arecibo 4 Vega_Baja 5 Hatillo 6 Camuy 7 Quebradillas 8 Barceloneta 9 Vega_Alta 10 Manati 11 Dorado . . . 67 Penuelas 68 Yabucoa 69 Patillas 70 Guayama 71 Salinas 72 Lajas 73 Maunabo 74 Santa_Isabel 75 Arroyo 76 Guanica scale concatenated digitized boundaries ArcInfo sequence of integers with names consecutive integers with names

  5. File entries • The boundaries file generated by ArcGIS has areal units in the same order as appears in the *.shp attribute file (i.e., the *.dbf file) • The order of the data included in the GeoBUGS files must be the same as this order • The names should be in alphabetical order; WinBUGS checks that the pairs of alpha-numeric names are exactly the same • REMINDER: WinBUGS default is the queen’s definition of adjacency

  6. Data entries • Each variable needs to be extracted from the ArcView/ArcGIS *.dbf file, inserted into MSWord, and then converted to a sequence of values (with line breaks) that are comma-delimited. • The MSWord replacement (use global search and replace) command for paragraph is ^p, which needs to be replace with a space followed by a comma

  7. Model syntax • The model (e.g., binomial probability, Poisson log-mean) needs to be specified • The prior distribution for each parameter, including its variance, needs to be specified: dflat() should be used for an intercept term • The empirical variable values need to be included • A random effectsMAY be included • Initial values should be posited • Spatial autocorrelation can be included via a CAR model, an ICAR model, or a SF model

  8. 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

  9. 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

  10. 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 will need to be substantially increased (say, to 25,000 first for burn-in, and then 500,000 for MCMC; once any desired changes are made, click once on the “updates” button

  11. 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: the burn-in period (say, 10,000) should be specified, as should the weeding (say, 100); click once of the “history” button to view the time-series plot for all 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

  12. 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 • 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

  13. Sample model syntax: example #1 #MODEL model { for (i in 1 : N) { U[i] ~ dbin(p[i], T[i]) p[i] <- exp(alpha)/(1+exp(alpha )) } # Other 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)

  14. Sample model syntax: example #2 #MODEL model { for (i in 1 : N) { U[i] ~ dbin(p[i], T[i]) p[i] <- 1/(1+exp(alpha + b*SF[i] + S[i])) S[i] ~ dnorm(0,0.0001) } # Other priors: alpha ~ dflat() b ~ dnorm(1,precb) precb ~ dgamma(0.5,0.0005) # prior on precision sigmab <- sqrt(1/precb) # standard deviation } #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 ),

  15. SF = c(-1.137518, 1.491700, 1.069692, 1.662677, -0.748245, 0.262096, -0.402159, -0.017504, -0.164057, 0.409578, 2.613552, 0.304349, 0.341998, -1.310238, 0.466386, -0.957026, 0.362092, -0.889712, 0.164100, -1.232902, 1.451166, 1.630955, 1.673424, 0.299560, -0.322737, -0.484777, -0.209917, 1.505092, 0.001573, -0.953810, 0.230846, 0.340325, 0.054938, -1.335408, -1.410285, -0.394615, 0.086024, -2.688523, -1.730705, -0.183643, -0.139794, 0.280239, -0.018101, -1.971644, 0.284003, -0.599356, 0.799498, 0.479685, 0.310890, 2.101611, -0.919954, -0.322202, -0.072132, -0.955641, -0.713250, 0.981736, 0.121009, -0.832610, -0.958571, -0.350464, 1.379148, -0.657765, -1.397524, -0.706984, 2.837594, 1.746381, 0.747843, -1.527760, 1.852185, 0.608635, -0.818181, 0.103756, -1.520626 ), ) #INITIAL VALUES list(alpha=-3, b=1, precb=0.001, 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) )

  16. Executing GeoBUGS: summary First the map must be loaded • Step 1: check the model syntax • Step 2: compile the model • Step 3: load the data • Step 4: load the initial values • Step 5: set the parameters to be estimated • Step 6: execute R iterations • Step 7: remove burn-in and weed • Step 8: correlogram, history, estimates

  17. Model syntax for CAR model #MODEL 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) }

  18. #DATA list(N = 76, U = c(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, 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 ), T = c(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, 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 ), num = c(6, 4, 3, 6, 5, 6, 5, 2, 3, 7, 7, 4, 7, 4, 5, 6, 3, 7, 4, 7, 6, 7, 5, 6, 3, 2, 4, 3, 4, 3, 4, 5, 4, 3, 3, 4, 5, 6, 5, 4, 8, 5, 7, 3, 3, 5, 6, 2, 6, 5, 6, 4, 5, 8, 6, 3, 5, 3, 2, 6, 5, 5, 6, 5, 7, 7, 3, 6, 4, 4, 7, 5, 3, 3, 5, 6 ),

  19. 65, 43, 32, 16, 15, 63, 61, 28, 12, 76, 71, 66, 47, 42, 33, 14, 1, 66, 49, 47, 41, 6, 75, 65, 60, 52, 39, 35, 15, 60, 16, 15, 60, 26, 19, 73, 51, 27, 20, 9, 76, 63, 61, 49, 42, 41, 75, 55, 63, 47, 42, 34, 12, 6, 66, 36, 6, 3, 2, 73, 72, 54, 46, 24, 20, 60, 43, 35, 19, 68, 24, 23, 11, 10, 74, 51, 38, 37, 24, 22, 20, 10, 75, 65, 48, 29, 18, 8, 57, 30, 1, 71, 56, 38, 37, 1, 66, 36, 14, 6, 2, 52, 45, 44, 43, 19, 15, 76, 63, 47, 40, 28, 67, 29, 22, 18, 5, 61, 49, 47, 40, 34, 12, 70, 31, 16, 13, 4, 75, 55, 43, 39, 32, 18, 13, 58, 50, 42, 41, 36, 14, 6, 62, 38, 22, 72, 69, 53, 25, 24, 11, 68, 25, 17, 11, 64, 32, 16, 13, 57, 41, 37, 33, 20, 7, 1, 73, 68, 51, 25, 24, 72, 51, 46, 54, 38, 22, 65, 55, 48, 43, 35, 61, 47, 41, 30, 28, 1 ), sumNumNeigh = 366) adj = c( 76, 71, 57, 56, 41, 30, 59, 50, 6, 3, 50, 36, 2, 64, 31, 23, 21, 13, 11, 62, 22, 21, 18, 10, 66, 59, 50, 49, 42, 2, 71, 33, 27, 20, 9, 55, 29, 46, 27, 7, 54, 53, 24, 23, 22, 21, 5, 69, 68, 53, 31, 23, 17, 4, 63, 49, 40, 34, 70, 65, 64, 32, 21, 18, 4, 66, 58, 41, 33, 60, 44, 43, 39, 16, 70, 64, 44, 39, 32, 15, 69, 31, 11, 65, 62, 55, 29, 21, 13, 5, 60, 52, 45, 26, 71, 54, 51, 46, 37, 27, 7, 23, 18, 13, 10, 5, 4, 74, 67, 62, 54, 38, 10, 5, 53, 21, 11, 10, 4, 72, 68, 54, 53, 51, 10, 72, 69, 68, 45, 19, 46, 20, 9, 7, 76, 61, 40, 62, 55, 18, 8, 76, 56, 1, 64, 17, 11, 4, 70, 65, 39, 16, 13, 71, 41, 14, 7, 63, 49, 12, 75, 52, 43, 66, 58, 50, 3, 71, 57, 54, 38, 20, 74, 67, 57, 54, 37, 22,

  20. #INITIAL VALUES list(prec=1, alpha=3, rho=0.1, 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) )

  21. Density of Puerto Rico population: spatial filter model syntax #MODEL model { for (i in 1 : N) { P[i] ~ dpois(mu[i]) offset[i] <- log(A[i]) log(mu[i]) <- alpha+beta*SF[i]+U[i]+(offset[i]-mean(offset[])) U[i] ~ dnorm(0,tau.U) } # Other priors: alpha ~ dnorm(0.0,1.0E-6) beta ~ dnorm(0.0,1.0E-6) tau.U ~ dgamma(1.0E-3,1.0E-3) # prior on precision sigma2.U <- 1/(tau.U*tau.U) # standard deviation }

  22. #DATA list(N = 73, P = 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 ), SF = c( -0.87980, 0.18990, 0.15104, 0.87052, -0.43893, 0.09154, -0.68331, -0.04907, -0.50465, -0.39226, 1.14477, 0.21607, 0.63555, -0.40064, 0.90284, -0.32918, 0.55484, -0.95903, -0.14897, -0.22116, 0.24814, 0.09793, 0.83398, 0.41653, -0.37741, -0.15932, -0.44595, 0.86744, 0.63944, -0.66191, 0.23526, -0.07711, 0.06443, -0.48191, -0.04134, 0.16879, -0.09492, -1.13629, -0.33841, 0.04491, 0.40173, 0.51533, -0.50619, -0.62497, -0.16845, 0.07312, 0.12536, -0.40689, 0.24930, 0.43163, -0.52047, -0.34204, -0.27139, -0.38079, -0.05208, 0.14143, 0.48815, -0.49467, -0.15929, -0.02422, 1.14529, -0.07535, -0.28340, 0.12182, 1.11191, 0.96355, 0.93429, -1.04985, 0.46577, -0.09707, -0.00740, -0.30082, -0.95594 ), A = c( 17384.0, 8011.8, 9489.2, 7911.6, 8110.9, 10240.6, 33304.0, 3897.9, 8861.7, 8867.7, 12857.7, 18873.8, 15287.3, 12105.0, 12490.4, 13454.8, 7683.3, 17283.9, 9443.4, 20205.6, 7391.4, 11020.7, 6092.9, 8004.2, 9670.5, 17076.5, 11194.5, 7023.4, 7333.5, 10893.7, 2932.9, 11632.2, 14404.2, 11527.1, 15734.5, 6893.1, 15792.5, 15965.8, 12037.5, 8790.6, 13896.7, 6705.8, 12048.6, 9478.6, 5470.8, 14314.6, 13042.7, 10085.1, 13443.7, 7170.6, 16525.0, 12230.2, 11374.8, 30047.0, 5988.9, 3695.6, 15759.1, 9280.3, 17747.6, 14094.2, 12909.2, 13775.6, 18458.6, 8793.6, 7136.4, 6133.4, 5554.1, 29783.6, 7346.9, 12503.6, 9576.9, 14310.8, 17802.0 ) )

  23. #INITIAL VALUES list(alpha=10.6, beta=1.0, tau.U=0.2, 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) ) From SAS

  24. Comparison: SAS PROC NLMIXED & GeoBUGS, with a SF

  25. SF: mean = -0.0000, std = 0.5391 MC = 0.8488, GR = 0.2535 • WinBUGS Bayesian priors: intercept, beta ~ N(0.0,0.000001) sre ~ gamma(0.001,0.001) • WinBUGS/random effects (RE; 125 replicates): 525,000 iterations 25,000-iteration burn-in 400 weeding average RE

  26. What you should achieve with today’s lab: • An appreciation of GeoBUGS

More Related