data("WEF.dat")
cleandata <-subset(WEF.dat, East_m != "NA")
coordinates(cleandata)<- c("East_m","North_m")
table(cleandata@data$Species)
## 
##   DF   GF   NF   SF  UNK   WH 
##  360   40    5 1336   34  539
douglasfirs.dat <- subset(cleandata, Species =="DF")

plot(douglasfirs.dat$East_m, douglasfirs.dat$North_m, xlab= "East", ylab="North", main="Location of Douglas Firs")

douglas.var <-variogram(douglasfirs.dat$DBH_cm~1, douglasfirs.dat,cutoff=300)
plot(douglas.var, main="Empirical Variogram for DBH")

douglas.var$np
##  [1] 2705 4763 5945 6894 7363 6862 5399 4491 3756 3412 2891 2508 2092 1828 1659

Empirical variogram seems to reach the partial sill at around 1200 and then begins to pick up a trend, there is a substantial nugget. Variogram suggests maybe a spherical or exponential model.

plot(douglasfirs.dat$ELEV_m, douglasfirs.dat$DBH_cm, xlab="Elevation", ylab="Diamater of Trunk (cm)", main="Elevation vs Diamater")  

douglasfirs.dat$ELEVcent <- douglasfirs.dat$ELEV_m - mean(douglasfirs.dat$ELEV_m)
elevmodel <- lm(DBH_cm ~ ELEVcent, data=douglasfirs.dat)
douglasresid <- add_residuals(douglasfirs.dat,elevmodel, var="resid")

douglasresid.var <- variogram(douglasresid$resid~1, douglasresid, cutoff=200)
plot(douglasresid.var, main="Detrended Empirical Variogram")

Empirical variogram reaches the partial sill at around 1000, there is a very large nugget. Variogram suggests exponential or spherical.

douglasmodel <- lgm(DBH_cm~ELEVcent, douglasfirs.dat, grid=60, boxcox = 1, fixBoxcox = T, aniso=F,fixShape=F)
summary(douglasmodel)
##                     estimate    stdErr      ci0.005      ci0.995      ci0.025
## (Intercept)        93.384683 0.4061381 92.338540448 9.443083e+01 92.588666815
## ELEVcent            3.757761 0.0385900  3.658359770 3.857162e+00  3.682126015
## sdNugget           20.662866        NA  0.027491230 1.553055e+04  0.133914508
## sdSpatial          26.692301        NA  0.274100305 2.599337e+03  0.819110151
## range             196.824018        NA  0.000216557 1.788891e+08  0.005757252
## shape               0.100000        NA -1.298331514 1.498332e+00 -0.963998846
## anisoRatio          1.000000        NA           NA           NA           NA
## anisoAngleRadians   0.000000        NA           NA           NA           NA
## anisoAngleDegrees   0.000000        NA           NA           NA           NA
## boxcox              1.000000        NA           NA           NA           NA
##                        ci0.975      ci0.05       ci0.95      ci0.1        ci0.9
## (Intercept)       9.418070e+01 92.71664514 9.405272e+01 92.8641959 9.390517e+01
## ELEVcent          3.833396e+00  3.69428613 3.821236e+00  3.7083060 3.807216e+00
## sdNugget          3.188258e+03  0.30106377 1.418151e+03  0.7661203 5.572937e+02
## sdSpatial         8.698207e+02  1.43416694 4.967894e+02  2.7356461 2.604427e+02
## range             6.728851e+06  0.03084185 1.256075e+06  0.2135721 1.813893e+05
## shape             1.163999e+00 -0.79293598 9.929360e-01 -0.5957115 7.957115e-01
## anisoRatio                  NA          NA           NA         NA           NA
## anisoAngleRadians           NA          NA           NA         NA           NA
## anisoAngleDegrees           NA          NA           NA         NA           NA
## boxcox                      NA          NA           NA         NA           NA
##                   pval Estimated
## (Intercept)          0      TRUE
## ELEVcent             0      TRUE
## sdNugget            NA      TRUE
## sdSpatial           NA      TRUE
## range               NA      TRUE
## shape               NA      TRUE
## anisoRatio          NA     FALSE
## anisoAngleRadians   NA     FALSE
## anisoAngleDegrees   NA     FALSE
## boxcox              NA     FALSE

Confidence intervals are calculated using MLE and \(\chi ^2\) distributions, 90% confidence intervals are estremely large for range shape and nugget.

douglasfitted <- fit.variogram(douglasresid.var, vgm(model="Mat", range=196, kap=.1, nugget=441, psill=729), fit.ranges = F, fit.kappa = F, fit.sills = c(F,F))


plot(douglasresid.var, douglasfitted, main="Fitted Variogram")

par(mfrow = c(1,2))
plot(douglasmodel$predict[["predict"]], main="Predicted Diameter")
plot(douglasmodel$predict[["krigeSd"]], main = "Kriging Standard Error")