#purpose of this exercise: empirical variograms, model fitting using Least squares methods

#setwd("C://stat406")#fixes the path leading to a certain directory so you don't have to specify it every time you need to import or export something from this R session

#read in data  from a web link
Lemp.table<-read.table("http://www.public.iastate.edu/~pcaragea/S40608/Compute/Emp_example.txt")
Lemp<-as.geodata(Lemp.table)

#plot observed locations
points(Lemp)

#and get summaries:
summary(Lemp)

#compute empirical variogram
# I will use 30 bins
vg.deft<-variog(Lemp,uvec=30)
#n is the number of pairs going into each class
vg.deft$n#nr. of pairs in each bin
vg.deft$uvec #center for each bin

#I will use only distances up to 240:
vg.deft<-variog(Lemp,uvec=25,max.dist=240)
vg.deft$n
vg.deft$uvec

#and you can plot it using the 'plot' function:
plot(vg.deft,pch=20)
#to plot the number of pairs in each bin on the variogram plot:
text(pos=1,vg.deft$u, vg.deft$v,labels=as.character(vg.deft$n))#not recommended for presentations, only helping us see better the nr. of pairs

#fit a model to the empirical variogram using LEAST SQUARES
#need to specify the sample variogram in vario
#need to specify starting values in ini.cov.pars for partial sill and range: I usually use the empirical variogram to guess these initial values. Always should try to assess sensitivity to starting values, try several and hope to see not much difference. 
#need to specify a model to try in cov.model (many models available in R, most useful gaussian, exponential, spherical, matern, wave, pure.nugget)
#need to specify weights (3 weights available in R: npairs=nr. pairs in each bin, cressie=closest to GLS, and equal=OLS)
#you could specify max distance here, and then only the bins up to that distance would be used to fit the desired theoretical model

exp.vg<-variofit(vario=vg.deft,ini.cov.pars=c(2,30),cov.model="exponential",weights="npairs")
#to see the resulting objects in the fitted model, type 
names(exp.vg)
#or, to see certain components in more detail, type
exp.vg$nugget #gives the estimate of the nugget
exp.vg$cov.pars #gives the estimates of the partial sill(sigma2) and the range
#to get the actual sill, you could add the nugget to sigma2 as:
sill.exp=exp.vg$nugget+exp.vg$cov.pars[1]

#to plot the theoretical model on top of the empirical variogram I would type
plot(vg.deft, pch=20, xlab="DISTANCE",ylab="GAMMA")#plots the empirical variogram. Change labels to something a little bit more meaningful, or at least better looking
lines(exp.vg)#plots the fitted model on top of the empirical variogram

#try other weights
#cressie
exp.vg.cr<-variofit(vario=vg.deft,ini.cov.pars=c(2,30),cov.model="exponential",weights="cressie")
exp.vg.cr$nugget 
exp.vg.cr$cov.pars 
lines(exp.vg.cr,col="red")#plots the new fitted model on top of the empirical variogram in RED


#equal (OLS)
exp.vg.ols<-variofit(vario=vg.deft,ini.cov.pars=c(2,30),cov.model="exponential",weights="equal")
exp.vg.ols$nugget 
exp.vg.ols$cov.pars 
lines(exp.vg.ols,col="blue")#plots the new fitted model on top of the empirical variogram in BLUE


#try other models
sph.vg<-variofit(vario=vg.deft,ini.cov.pars=c(2,30),cov.model="spherical",weights="npairs")
sph.vg.cr<-variofit(vario=vg.deft,ini.cov.pars=c(2,30),cov.model="spherical",weights="cressie")
plot(vg.deft, pch=20, xlab="DISTANCE",ylab="GAMMA")
lines(exp.vg)
lines(sph.vg,col="red")
lines(sph.vg.cr,lty=2,col="red")#change the line type as well as color