Part A

library(geoR)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------
library(gstat)
as.dat= dget("C:\\Users\\Carly\\Downloads\\as-1.data")
as.geor=as.geodata(as.dat)



variogram<-variog(as.geor)
## variog: computing omnidirectional variogram

h is the uvec vector in the variogram:

(h<-variogram$uvec)
##  [1] 0.05185866 0.15557597 0.25929328 0.36301059 0.46672790 0.57044521
##  [7] 0.67416252 0.77787983 0.88159714 0.98531445 1.08903177 1.19274908
## [13] 1.29646639

delta is the width of one bin:

((delta<-variogram$uvec[4]-variogram$uvec[3]))
## [1] 0.1037173

Variogram plot:

plot(variogram)

Part B

fit_exp<-variofit(variogram, ini=, fix.nugget = F, weights = "npairs", nugget=1, cov.model="exponential")
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(variogram, ini = , fix.nugget = F, weights =
## "npairs", : initial values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.39"  "0.83" "0.05" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 132.525519702429
fit<-summary(fit_exp)
fit$estimated.pars
##        tausq      sigmasq          phi 
## 6.520157e-02 2.210211e+03 7.206871e+03

The estimates of theta are \(\theta_1 = 2210.211\) , $_2 = 7206.9 $ , and the estimated nugget effect is 0.0652.

fit_power<-variofit(variogram, ini=, fix.nugget = F, weights = "npairs", cov.model="power")
## variofit: covariance model used is power 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
## Warning in variofit(variogram, ini = , fix.nugget = F, weights =
## "npairs", : initial values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "0.26"  "1.04" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 1138.14769446855
fit<-summary(fit_power)
fit$estimated.pars
##     tausq   sigmasq       phi 
## 0.1171672 0.2814871 1.6507010

The estimates of theta are \(\theta_1 = 0.281\) , \(\theta_2 = 1.651\) , and the estimated nugget effect is 0.117.

Part C

fit_exp_c<-variofit(variogram, ini=, fix.nugget = F, weights = "cressie", nugget=1, cov.model="exponential")
## variofit: covariance model used is exponential 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
## Warning in variofit(variogram, ini = , fix.nugget = F, weights =
## "cressie", : initial values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.39"  "0.83" "0.05" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 1745.16875744673
fit<-summary(fit_exp_c)
fit$estimated.pars
##        tausq      sigmasq          phi 
## 7.278071e-02 6.847142e+02 2.311238e+03

The estimates of theta are \(\theta_1 = 684.7\) , \(\theta_2 = 2311.2\) , and the estimated nugget effect is 0.0728.

fit_power_c<-variofit(variogram, ini=, fix.nugget = F, weights = "cressie", cov.model="power")
## variofit: covariance model used is power 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
## Warning in variofit(variogram, ini = , fix.nugget = F, weights =
## "cressie", : initial values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.26"  "1.04" "0.05" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 18990.7190944408
fit<-summary(fit_power_c)
fit$estimated.pars
##     tausq   sigmasq       phi 
## 0.0726233 0.2966663 1.0000000

The estimates of theta are \(\theta_1 = 0.297\) , \(\theta_2 = 1\) , and the estimated nugget effect is 0.0726.

To qualitatively judge which model is providing a better fit to the data, we can look at which model has the lowest minimised weighted sum of squares:

fit_exp$value
## [1] 66.96165
fit_power$value
## [1] 51.24559
fit_exp_c$value
## [1] 1360.536
fit_power_c$value
## [1] 1360.487

The “npairs” weighted power variogram model appears to provide the best fit to the data.

Standard Errors

If we wanted to find the standard error of our estimates, we could find the information matrix (I think one would have to use numerical methods we do not yet know), then the diagonal of 1 / Information matrix would provide us with an estimate of the variance. We can do this because we know that the MLE is approximately normal.

Also, after researching this online I found that there are cross-validation methods that can be used to provide this estimation. The standard deviation of the cross validation residuals will be of use when finding the standard error of our estimates.