We will discuss the use of SAS and R for more sophisticated analysis of spatial data. SAS is expensive software available for yearly rent through the Stat/Math center on campus or on campus computers for free. R is free software that can be found at:
The R Project for Statistical Computing
geoR is a contributed package to R and can loaded from the R program once you have installed R. A guide to geoR is available at:
http://www.est.ufpr.br/geoR/geoRdoc/geoRintro.html
Here are the commands we may use in class are as follows. They load the library, read the data, print and plot the data, create and plot a binned variogram, and determine the maximum likelihood estimates for a gaussian model.
library(geoR)
rain <- read.geodata('SwissRain.txt', header=FALSE, coords.col = 1:2, data.col=3)
summary(rain)
plot(rain)
bin <- variog(rain, option = "bin", uvec=seq(0,340000,10000))
plot(bin)
like <- likfit(rain, ini = c(10000, 31000), nugget=3000, fix.kappa=TRUE, kappa=2, cov.model="powered.exponential", control=list(parscale=c(10000,.1)))
like
/* Name the dataset "rain" */
DATA rain;
/* Tell SAS where to get the data */
INFILE 'C:\Documents and Settings\Elizabeth Housworth\My Documents\Math468\SwissRain.txt';
/* Tell SAS the variables names */
INPUT x y rainfall;
/* Execute the code above */
RUN;
/* Most things in SAS are procedures. We want to run the printing procedure
to make sure we loaded the data properly. */
Proc print Data=rain;
Run;
/* The variogram procedure computes the variogram. The input is out data set
rain. The output is a new data set called vario. We then tell later
procedures which data set to use - rain or vario. */
PROC VARIOGRAM DATA=rain OUTVAR=vario;
/* You control the binning process with this statement. LAGD is the distance
for the bin and MAXLAGS is the maximum number of bins. */
COMPUTE LAGD=8000 MAXLAGS=42;
/* COORDINATES tells SAS which of the variables in the DATA set are the
X-coordinates and which are the Y-coordinates. */
COORDINATES xc=x yc=y;
/* VAR tells SAS which variable contains the spatial responses. */
VAR rainfall;
RUN;
/* Check out what you created. Notice that SAS named the
variogram estimates VARIOG. */
PROC PRINT DATA=vario;
RUN;
/* Now plot it. */
PROC GPLOT DATA=vario;
/* These two statements define the sizes of the vertical and
horizontal axes, respectively. These are assigned to the
vertical and horizontal axes by the options in the PLOT
statement. Be careful to ensure that they include the range for your data. */
AXIS1 ORDER=0 TO 20000 BY 2000;
AXIS2 ORDER=0 TO 350000 BY 20000;
/* The PLOT command says to plot variog on the y axis versys distance on the x axis */
PLOT variog*distance / VAXIS=axis1 HAXIS=axis2;
RUN;
/* Now we want an initial fit of some model to the variogram. SAS
does not seem to be able to work with the Matern distribution easily.
It can handle spherical, exponential (powered exponential with k=1),
gaussian (powered exponential with k=2), the powered exponential generally,
and some other possibilities: linear, log linear, etc...
Let's fit a general powered exponential to the data: */
/* The METHOD tells SAS which optimization algorithm to use
to minimize the sum of squares. The options are: */
PROC NLIN METHOD=dud;
/* This statement defines the model parameters and gives
them initial values for SAS to use to start the optimization
algorithm. You need to start with reasonable parameters.
You could use VarioWin of just the variogram plot from
before to get some intial guesses. */
PARAMETERS tau2=2000 sigma2=15000.0 psi=50000.0 kappa=1;
/* This statement places bounds on those parameters to make
sure the results don't wander into invalid values. *.
BOUNDS c>0, a>0;
/* This statement defines the exponential model. */
poweredexp=tau2 + sigma2*(1-exp(-((distance/psi)**kappa)));
MODEL variog=poweredexp;
/* This gives the values of the weights to use. The
underscores around the variable name mean that it is a
special SAS variable name. In this case, the weight is the
inverse of the variance of the model. */
_WEIGHT_=count/(poweredexp);
/* This puts the estimated model values into a data set
called VARIOMOD*/
OUTPUT OUT=variomod PRED=evhat;
RUN;
/* Now plot out the two model fits against the sample
variogram. */
PROC GPLOT DATA=variomod;
AXIS1 ORDER=0 TO 20000 BY 2000;
AXIS2 ORDER=0 TO 350000 BY 20000;
/* This is one way to get several plots on the same axes. */
PLOT (variog evhat)*distance /
OVERLAY VAXIS=AXIS1 HAXIS=AXIS2;
/* These statements tell SAS how to plot the different data
sets. The first draws black bullets for the first variable,
VARIOG; the second draws a solid red line (L=1) through the
model estimates EVHAT; The line
control L can take on values up to 46, with values 2-46
giving a different kind of dashed line. */
SYMBOL1 VALUE=J FONT=special INTERPOL=none COLOR=black;
SYMBOL2 VALUE=J FONT=special INTERPOL=none COLOR=red;
RUN;
/* Now we have some decent parameter estimates:
Parameter Estimate Std Error Approximate 95% Confidence Limits
tau2 3269.4 753.0 1746.3 4792.6
sigma2 10382.9 798.9 8767.0 11998.8
psi 41361.9 2304.9 36699.9 46024.0
kappa 2.7628 0.6078 1.5334 3.9921
We will use thess in Procedure MIXED to come up with the maximum likelihood estimates
for all the parameters (including the overall mean) in the model. */
DATA rain;
/* Tell SAS where to get the data */
INFILE 'C:\Documents and Settings\Elizabeth Housworth\My Documents\Math468\SwissRain.txt';
/* Tell SAS the variables names */
INPUT x y rainfall;
/* Execute the code above */
RUN;
/* Proc Mixed handles general mixture models. Mixture refers to, among other things, that
our variance matrix for the data is a mixture (sum) of a variance times the correlation matrix
for the signal and another variance times the identity matrix. */
PROC MIXED METHOD=ml;
/* The model can be a regression model. We just want to estimate the mean so our model is empty.
That means the model is rainfall = intercept.
An intercept (mean) is assumed by SAS. */
MODEL rainfall= / SOLUTION;
/* Here is where we specify the form for the variance matrices.
The SUBJECT tells what is fitting the covariance model
LOCAL tells SAS to add tau^2 * I to the covariance structure
TYPE=SP(Gau)(x y) tells SAS to use a spatial model, specifically the
Gaussian one, with spatial coordinates given by the variables x and y */
REPEATED / SUBJECT=INTERCEPT LOCAL TYPE=SP(Gau)(x y);
/* These are the inital parameters in the model: sigma_squared, psi, and then tau_squared */
PARMS (10382) (41361) (3269);
RUN;
The Swiss rainfall dataset is given below. It is available in many places. I found it at
http://www.cas.umt.edu/math/graham/math595/math595.html
x y rainfall -159812 -39393 215 -159806 -68316 167 -158368 -42768 271 -156586 -62838 166 -153432 -29548 263 -153090 -46242 149 -148708 -24105 190 -148491 -46357 187 -146620 -64200 220 -145922 -35296 149 -143939 -16435 228 -140463 -30977 151 -136211 -12166 255 -129848 -72 210 -128276 -33474 210 -127684 -5680 234 -126536 -23500 110 -125353 -3505 217 -124685 -7968 100 -124038 -13542 120 -123184 -9111 79 -121276 9758 191 -120674 1960 194 -119464 -13635 133 -118084 17479 234 -117307 -20351 214 -116742 8555 183 -116780 -32596 334 -115480 -4816 107 -114528 -34864 371 -113301 -10418 230 -111242 -22691 290 -110314 27340 228 -109680 20655 296 -109192 -36075 394 -109172 -34963 394 -108729 -10503 300 -108027 -13852 278 -107305 26173 260 -104273 -17256 371 -103712 14987 225 -103460 29441 281 -101596 -39544 323 -101256 -65129 338 -96756 22655 290 -97055 -42955 297 -96147 -34073 383 -95716 -55210 255 -94512 -27426 432 -93956 -40780 415 -93540 -14096 324 -93350 -1866 301 -93332 -754 322 -92287 65956 105 -92261 19247 223 -93102 -85277 82 -89885 25884 119 -88961 36991 251 -88794 48109 266 -89164 -79776 99 -88936 -64210 214 -87683 -30866 154 -86804 -76474 83 -85629 5803 289 -84732 68067 135 -85061 -62042 214 -84906 -50923 213 -83936 18012 330 -83716 33578 96 -83690 -18690 585 -82882 39127 194 -83783 -80965 114 -82771 -63186 220 -81792 -47630 268 -80712 -24291 407 -79712 -6511 493 -79475 11279 342 -78818 3485 382 -79258 -88810 33 -78609 -38776 264 -78030 -53240 180 -75594 17901 334 -75197 49033 209 -74547 -19922 390 -73962 -97772 98 -72356 34541 180 -73017 -83327 57 -71263 62329 230 -70390 8940 347 -69838 -8859 415 -68588 33384 256 -69022 -71142 51 -68707 -43345 131 -68015 -50025 165 -67728 -24452 283 -67188 23360 300 -66815 -11118 405 -64428 67810 196 -64887 -44499 104 -63854 -92327 78 -63599 -67865 45 -62900 -1152 393 -62677 19974 398 -61617 48876 208 -61512 58883 176 -60921 -31197 121 -59304 -22317 174 -59026 -72360 98 -56696 11018 434 -54591 -7906 351 -52692 34334 349 -52257 -2366 406 -51974 -57968 99 -51393 -80213 39 -50660 6516 434 -50299 -41303 135 -49564 46539 288 -49480 56546 267 -49688 -60211 47 -49277 -10175 221 -48722 -35756 119 -47282 48744 258 -47061 -19089 198 -46443 -37998 141 -45799 45396 306 -45574 74306 136 -45522 80978 108 -45400 -1310 192 -45168 -72479 43 -44038 78742 123 -43485 52051 304 -43354 69841 102 -42618 67612 153 -42594 70948 151 -42545 77619 133 -42513 82067 107 -42749 -56929 55 -42159 26466 375 -42112 33137 328 -42112 -75836 65 -41935 -50262 52 -41642 -8010 200 -41067 74273 142 -41044 77608 148 -40170 -15803 143 -39542 77598 152 -39481 -26928 97 -39247 8654 517 -38077 72029 151 -38004 83148 152 -38015 -36945 121 -37133 -18047 142 -35931 54224 310 -35877 63119 320 -35257 40876 350 -34132 -26960 145 -33999 -3610 223 -33631 60882 400 -33561 73113 253 -33529 78673 161 -33071 26408 445 -32168 -89241 108 -31367 61982 386 -30591 66425 289 -29884 57526 334 -29097 64194 325 -29037 76425 260 -28596 11929 429 -28425 47511 327 -28312 70862 325 -27780 24157 378 -27723 36388 322 -27345 -47010 37 -26794 74190 315 -26758 81974 213 -26097 61956 356 -26043 74187 338 -25311 69736 331 -25141 -70371 0 -25096 -59252 17 -24768 19696 412 -24717 31927 400 -24666 44159 311 -24508 81964 238 -23978 27477 327 -23234 24138 426 -22628 -15895 82 -22047 -65935 13 -21560 67497 380 -21153 -30356 123 -20973 20793 365 -20275 3000 251 -20090 57484 301 -20059 66380 350 -19610 -25914 94 -18685 25234 314 -18597 53032 322 -17799 67485 335 -17218 8550 185 -17186 19669 239 -15935 -81522 58 -15792 -27037 118 -14797 65252 367 -14028 71922 330 -13567 -57065 30 -12630 22994 225 -11145 6311 145 -11083 39670 196 -11016 75252 324 -8908 -22603 141 -8209 -70417 98 -8127 -9261 68 -7332 17425 166 -6524 61902 254 -3518 53004 226 -2812 -49296 20 -2045 -48184 34 -2004 71906 266 -508 -4818 72 254 -1483 53 251 70794 286 1790 -51520 78 1778 -11490 73 1765 34100 115 1762 41883 212 1759 55227 216 4803 18534 142 4779 47445 165 5503 80804 225 6333 6304 57 6322 16312 71 8627 -3700 90 9412 -19266 140 10104 22989 87 10969 -35943 179 10774 70804 212 12486 -29268 148 12324 48568 178 13152 15211 32 13024 73032 183 14022 -32601 62 13971 -10362 71 16178 18554 59 16790 69706 132 18515 -1454 60 18311 64151 164 19006 81944 124 19861 50812 127 20697 27464 70 21414 38586 60 22095 58604 163 22045 71947 147 21952 96409 156 22714 93076 124 23427 101974 138 24830 -57029 354 24364 56389 124 25443 -20332 185 25738 86417 153 26705 38608 63 26674 45279 73 26579 66406 146 27305 71969 133 27284 76417 144 28253 30831 93 28158 50846 120 29878 8600 59 29728 38622 85 29582 67532 143 30369 60865 75 31263 34182 137 31216 43078 86 32301 -18075 129 31941 48641 124 31917 53089 77 33278 -56988 402 33150 -33638 345 33088 -22518 149 34145 -74775 441 33209 90903 113 34603 -20286 262 33933 95354 120 33874 105361 184 35222 4181 33 35084 27532 86 36429 -70314 444 35888 19752 86 37632 102049 121 39277 -34713 346 38768 43122 67 39590 33120 94 39426 57582 106 39239 85380 96 40111 67594 94 41878 -79174 278 41380 -6899 270 41899 27576 113 42166 95408 100 44129 33152 65 45012 16479 106 44793 45389 45 46026 -16873 209 45513 49842 66 45436 59849 115 45341 72080 104 46624 4259 171 46477 23162 91 46407 32057 91 46372 36505 58 45988 85429 107 48000 -73569 308 47974 -70233 312 47728 -39100 359 47261 19832 140 46630 98778 126 49386 -54654 278 48503 53202 58 49452 29858 72 50350 13186 162 50086 44319 48 51878 -76873 208 52863 -13481 171 52391 39891 54 53046 51017 120 54930 -73510 227 53719 59919 119 53526 81045 123 53434 91052 139 54545 52142 88 54225 86611 131 56563 -83503 90 55465 34359 82 55434 37694 64 56485 6568 190 56253 31030 103 55862 72171 141 55851 73283 131 59158 -33438 452 58337 49954 65 58991 59969 118 61279 -92353 74 59947 39962 86 59710 63313 89 61646 -52317 503 60382 71103 134 60291 79998 129 61816 5509 247 61111 73335 126 60902 93348 131 63769 -109008 0 63585 -91218 16 62727 -8938 287 63265 12196 199 63147 23316 155 64917 -71188 240 63172 91149 136 66081 -107871 1 64505 37786 117 64818 77822 109 67605 -105630 0 67556 -101183 0 65718 64488 130 65456 87838 118 67298 -8888 343 66647 48930 91 66256 83399 109 68438 24486 145 68585 76753 142 68402 92319 142 69848 33399 184 71226 -18851 254 70784 81227 140 72860 -27727 285 72550 -2154 221 71928 48992 152 71629 73454 128 74004 -58851 164 73312 59018 60 75610 -64391 196 74401 32342 115 76497 -13226 269 76199 10123 271 77388 -23223 348 76383 54608 95 78368 -39891 261 78296 -34332 283 76801 80191 90 77816 60187 63 77743 65746 61 79458 -6515 241 81551 -48745 257 80645 73570 68 81859 40225 184 81960 86933 127 84694 230 220 83475 85842 115 84689 53610 72 84561 62504 126 84305 80294 109 85540 46949 178 88193 -30859 276 87981 -16404 218 86885 58090 137 87887 41423 140 87738 51430 118 88559 46994 116 89751 18097 284 89094 61459 119 92088 62617 111 94583 378 298 93700 55970 185 93505 68200 125 95768 -26294 94 95450 -6281 313 94901 74895 144 97099 -14039 132 96720 9309 230 96467 24875 282 98187 -34040 129 97523 52695 231 97339 63814 140 98769 68285 116 100124 32720 235 99992 40503 205 99917 44950 311 101069 21615 200 101815 66113 149 103672 2753 170 105420 -53936 58 105188 -40595 65 103398 61692 392 103261 69475 303 106209 -11660 186 105735 15023 240 107070 -17206 190 107040 68430 249 109500 -24948 190 109815 -474 72 111743 18469 175 113592 -39330 96 113152 -15982 170 111636 64067 156 114546 -49321 131 115226 -4821 126 114798 17415 205 116224 -17036 184 120891 -21392 10 122870 -43597 99 122236 -12467 82 125007 937 92 126724 -44629 36 128096 -1222 67 131057 -32302 16 130476 -5620 100 132684 -36715 39 131778 4419 121 138143 -41041 40 142077 -45399 18 143715 -49810 30 145302 -51996 20 145757 -7492 8 146655 -13031 6 150921 -61870 55 154392 -17287 5 160524 -18238 29 161600 -1523 0 172891 -23457 15