Dear Thierry, Thank you. You're right. I have small data sets, which are problematic. I need to use the duplicate pair data, but the default weighting in gstat for variograms doesn't work well with closely spaced data.
I am looking now at REML instead. Best wishes, John John H. Carson Jr., PhD Senior Statistician Applied Sciences & Engineering Shaw Environmental & Infrastructure 16406 US Rte 224 East Findlay, OH 45840 Phone 419-425-6156 Fax 419-425-6085 john.car...@shawgrp.com http://www.shawgrp.com/ Shaw(tm) a world of Solutions(tm) -----Original Message----- From: ONKELINX, Thierry [mailto:thierry.onkel...@inbo.be] Sent: Monday, October 12, 2009 10:23 AM To: Carson, John; r-sig-geo@stat.math.ethz.ch Subject: RE: [R-sig-Geo] Problem with gstat variogram estimation Dear John, Isn't your problem rather a result from the unstable empirical variograms? Mainly due to the very low number of pairs per bin. Furthermore have a look at the weights that each bin gets in the fit.variogram() function. The default is N/h^2. In your case the first bin gets a weight of 2 whereas the other bins have weights ranging from 1e-3 tot 1e-4! Hence no surprise that the nugget is nearly identical to the semivariance of the first bin. In this case I would not trust the results for fit.variogram(). Not because of bugs in gstat, but because I don't trust empirical variogram with that low number of pairs per bin. HTH, Thierry ---------------------------------------------------------------------------- ir. Thierry Onkelinx Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology and quality assurance Gaverstraat 4 9500 Geraardsbergen Belgium tel. + 32 54/436 185 thierry.onkel...@inbo.be www.inbo.be To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger Brinner The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey -----Oorspronkelijk bericht----- Van: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] Namens Carson, John Verzonden: maandag 12 oktober 2009 15:35 Aan: r-sig-geo@stat.math.ethz.ch Onderwerp: [R-sig-Geo] Problem with gstat variogram estimation I have found anomalous behavior in gstat's variogram estimation. I have listed 3 example variograms below for small data sets. In order to better estimate the nugget effect, I slightly perturbed the locations (by 1 foot increments) of duplicate results. The empirical variograms are given below. Before I did this (I averaged duplicate values initially), a Gaussian model with 0 nugget was selected for the second variogram and pure nugget models for the first and third. I am using the candidate model list ('Nug', 'Exp', 'Sph', 'Gau', 'Mat', 'Cir', 'Lin', 'Bes') and selecting the model based on SSErr for preliminary testing purposes. Afterward, the pure nugget models had the lowest SSErr and were selected. Note that the variogram fits are completely controlled by the short range variance, because even the original pure nugget models are substantially different in the estimate of the nugget. The fitted models are listed below. Just by inspection, based on the numbers of pairs in these examples, a pure nugget model should be about halfway between the empirical semivariance of the last lag and the average of the other lags. However, the fitted nuggets are almost identical to the semivariance of the first last (dist = 1.4). It seems to me that this must be due to a bug in the GSTAT code. I pointed this out to Edzer Pebesma, and he asked me to post it here. The variograms are tmp.vgm [[1]] np dist gamma dir.hor dir.ver id 1 4 1.414214 0.14174537 0 0 PC1 2 2 44.742603 6.70989788 0 0 PC1 3 2 57.707880 1.76351594 0 0 PC1 4 4 59.987678 1.52197310 0 0 PC1 5 3 71.512518 1.21348268 0 0 PC1 6 1 84.852877 0.05381849 0 0 PC1 7 1 97.266495 1.21827622 0 0 PC1 8 3 112.237133 5.07947925 0 0 PC1 9 18 121.478856 1.93707676 0 0 PC1 [[2]] np dist gamma dir.hor dir.ver id 1 4 1.414214 0.09725079 0 0 PC2 2 2 44.742603 0.33598072 0 0 PC2 3 2 57.707880 0.39088727 0 0 PC2 4 4 59.987678 0.87315735 0 0 PC2 5 3 71.512518 0.14944845 0 0 PC2 6 1 84.852877 0.19809863 0 0 PC2 7 1 97.266495 0.63557814 0 0 PC2 8 3 112.237133 1.92063948 0 0 PC2 9 18 121.478856 0.65468693 0 0 PC2 [[3]] np dist gamma dir.hor dir.ver id 1 4 1.414214 0.035250817 0 0 PC3 2 2 44.742603 0.105299796 0 0 PC3 3 2 57.707880 0.020245674 0 0 PC3 4 4 59.987678 0.124159836 0 0 PC3 5 3 71.512518 0.008112554 0 0 PC3 6 1 84.852877 0.034337591 0 0 PC3 7 1 97.266495 0.053879459 0 0 PC3 8 3 112.237133 0.021922987 0 0 PC3 9 18 121.478856 0.085270969 0 0 PC3 But the fitted models are: tmp.vgm.fit [[1]] model psill range 1 Nug 0.1483120 0 [[2]] model psill range 1 Nug 0.09849419 0 [[3]] model psill range 1 Nug 0.03535234 0 John H. Carson Jr., PhD Senior Statistician Applied Sciences & Engineering Shaw Environmental & Infrastructure 16406 US Rte 224 East Findlay, OH 45840 Phone 419-425-6156 Fax 419-425-6085 john.car...@shawgrp.com http://www.shawgrp.com/ Shaw(tm) a world of Solutions(tm) ****Internet Email Confidentiality Footer**** Privileged/Confidential Information may be contained in this message. If you are not the addressee indicated in this message (or responsible for delivery of the message to such person), you may not copy or deliver this message to anyone. In such case, you should destroy this message and notify the sender by reply email. Please advise immediately if you or your employer do not consent to Internet email for messages of this kind. Opinions, conclusions and other information in this message that do not relate to the official business of The Shaw Group Inc. or its subsidiaries shall be understood as neither given nor endorsed by it. ______________________________________ The Shaw Group Inc. http://www.shawgrp.com _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo Druk dit bericht a.u.b. niet onnodig af. Please do not print this message unnecessarily. Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document. The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document. _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
