######################################################################
# Author: Michael Höhle <http://www.stat.uni-muenchen.de/~hoehle>
# Course: Space-Time Models in Epidemiology, SoSe 2009
#         Department of Statistics, University of Munich, Germany
#
#
# License:
#  Copyright 2009 Michael Hoehle
#  Distributed under GNU General Public License, v3
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Description:
#  R file for Spatial Statistics in Epidemiology course at the IBE
######################################################################

library("geoR")



#######
#(a)
#######

#Parameters of Z
mu <- 3
sigmasq <- 10

#Population
s1 <- c(0,0)
s2 <- c(1,1)
u <-  sqrt(sum((s1 - s2)^2))   
#Alternative:
#as.numeric(dist(rbind(s1,s2)))

#Matern covariance function
phi <- 1
(cu <- cov.spatial(u, cov.model="matern", cov.pars=c(sigmasq=10,phi=phi),kappa=0.5))
(rhou <- cu / sqrt(sigmasq * sigmasq))


#######
#(b)
#######

Sigma <- matrix( c(sigmasq, cu, cu, sigmasq), byrow=TRUE, 2,2)


#######
#(c)
#######

library("MASS")
z <- mvrnorm(n=1000, mu=rep(mu,2), Sigma=Sigma)
plot(z[,1],z[,2],asp=1)


#######
#(d)
#######

z.grf <- grf(n=2, grid= rbind(s1,s2), cov.model="matern", cov.pars=c(sigmasq,phi),kappa=0.5,nsim=1000,mean=mu)
plot(z.grf$data[1,],z.grf$data[2,],asp=1)
#plot(z.grf$data[1,]+3,z.grf$data[2,]+3,asp=1)
points(z[,1],z[,2],asp=1,col=2,cex=0.2)


#######
#(f)
#######

z2.grf <- grf(n=100, grid= "reg", cov.model="matern", cov.pars=c(sigmasq,phi),kappa=0.5,nsim=1,mean=mu)
image(z2.grf)


tcdd <- read.table('tcdd.txt',header=T, skip=9) 
borders <- rbind(c(min(tcdd$x),min(tcdd$y)),
                 c(min(tcdd$x),max(tcdd$y)),
                 c(max(tcdd$x),max(tcdd$y)),
                 c(max(tcdd$x),min(tcdd$y)))

tcdd <- as.geodata(tcdd, coords.col=1:2, data.col=3)
tcdd$borders <- borders


plot(0,type="n", xlim=range(tcdd$coord[,1]),ylim=range(tcdd$coord[,2]))
points(tcdd, col = gray(seq(1, 0.1, l = 100)),add=TRUE,cex.max=10)


par(mfcol=c(1,2))
hist(tcdd$data)
plot(tcdd$coord[,2],tcdd$data,xlab="y coordinate",ylab="log(conc)")


par(mfcol=c(1,2))
plot(variog(tcdd,option="cloud"))
plot(variog(tcdd,option="bin"))



tcdd$y <- matrix(tcdd$coord[,2],ncol=1)
m <- likfit(tcdd, trend = ~ I(y^2), cov.model="exp", ini.cov.pars=c(1,1))
#No spatial dependence
#trend <- lm( data ~I(y^2),data=tcdd)

m$beta[2] + c(-1,1) * 1.96 * m$beta.var[2,2]

par(mfcol=c(1,2))
plot(variog(tcdd,trend = ~ I(y^2),option="bin",direction=0),type="b",main="Angle 0")
plot(variog(tcdd,trend = ~ I(y^2),option="bin",direction=pi/2),type="b",main="Angle pi/2")


plot(tcdd$y,tcdd$data,xlab="y coordinate",ylab="log(conc)")
mu <- function(y) m$beta[1] + m$beta[2] * y^2
curve( mu, from=-20,35, n=1000, add=TRUE)


#http://www.math.umt.edu/graham/math544/empvariog.pdf
pred <- pred_grid(tcdd$borders, by=5)

KC <- krige.control(obj.model = m, trend.d = ~ I(y^2), trend.l = ~ I(pred[,2]^2))
k <- krige.conv( tcdd, loc=pred, krige=KC)
image(k,x.leg=c(-200,-100),y.leg=c(-100,-75),col=gray(seq(1,0.1,length=21)))
points(tcdd, col = gray(seq(1, 0.1, l = 100)),add=TRUE,cex.max=2)