# ======================================
# Fitting Generalized Linear Models
# ======================================

# Poisson distribution:
# p(x) = lambda^x exp(-lambda)/x!
# E(X) = Var(X) = lambda

rp  = rpois(1000,3)
summary(rp)
plot(rp)
hist(rp)

# -------------------------------------------------
# Spatial data
# n = 50 points in an specific area (unit square)
# -------------------------------------------------
n = 50
y  = rpois(n,3)      # random generation for Poisson (mean = lambda) 
summary(y)

xcoor = runif(n)
ycoor = runif(n)
plot(xcoor,ycoor,lwd=4)
text(xcoor,ycoor,1:n,cex=1,pos=2,col=4)

# Observations (map):                                  
col1 = heat.colors(max(y)+1)[(max(y)+1):1]
col2 = col1[y+1]                                        # colors for each value
aux  = max(y)/15
plot(0:1,0:1,type="n",xlab="x coord",ylab="y coord")    # square
points(xcoor,ycoor,bg=col2,pch=21,cex=(y+0.5)/aux)      # points with colors
text(xcoor,ycoor,y,cex=1.2)                             # values

# ----------------------------------
# Model:
#       y_i ~ Poisson(mu_i)
# log(mu_i) = b_0 + b_1 x_i
# ----------------------------------

# Simulated data:
n  = 50
x  = rnorm(n)
mu = exp(2 + 0.5*x)
y  = rpois(n, mu)
summary(y)
  
# Observations in a map:                                  
col1 = heat.colors(max(y)+1)[(max(y)+1):1]
col2 = col1[y+1]
aux  = max(y)/7
plot(0:1,0:1,type="n",xlab="x coord",ylab="y coord")
points(xcoor,ycoor,bg=col2,pch=21,cex=(y+0.5)/aux)
text(xcoor,ycoor,y,cex=1.2)



# Fit generalized linear model
fit = glm(y ~ x, family=poisson)
summary(fit)

# Obtains predictions:
p1 = predict(fit)
plot(y,exp(p1),ylab="observations",xlab="predictions",main="Model I")
abline(a=0,b=1,col=2,lty=2)


# ----------------------------------
# Model:
#    y_i ~ Gamma with mean equal to mu_i
# 1/mu_i = b_0 + b_1 x_i
# ----------------------------------

# Data
x = log(c(5,10,15,20,30,40,60,80,100))
y1 = c(118,58,42,35,27,25,21,19,18)
y2 = c(69,35,26,21,18,16,13,12,12)

# Modelo 1: 
fit1 = glm(y1 ~ x, family=Gamma)
summary(fit1)

# Obtains predictions:
p1 = predict(fit1)
plot(y1,1/p1,ylab="observations",xlab="predictions",main="Model I")
abline(a=0,b=1,col=2,lty=2)

# -----------------------------

# Modelo 2: 
fit2 = glm(y2 ~ x, family=Gamma)
summary(fit2)

# Obtains predictions:
p2 = predict(fit2)
plot(y2,1/p2,ylab="observations",xlab="predictions",main="Model II")
abline(a=0,b=1,col=2,lty=2)

# ========================================================
# Spatial model: 
# Gaussian distribution with spatial correlation
# ========================================================

install.packages("scatterplot3d")
install.packages("geoR")

# package: GeoR
library(geoR)
library(scatterplot3d)

# Simulated data 
s100
plot(s100)

summary(s100$data)
boxplot(s100$data)


# Observations in a map: 
xcoor = s100$coords[,1]
ycoor = s100$coords[,2]
y = s100$data

par(mfrow=c(1,2))
# plot 3D
scatterplot3d(xcoor, ycoor, y, type="h", 
              col.axis="blue", col.grid="lightblue",
              main="Spatial observations", pch=20,angle=50)

# another plot
points(s100, cex.min=.3, cex.max=4,col=heat.colors(15)) 

# Defining a prediction grid
loci <- expand.grid(seq(0,1,l=21), seq(0,1,l=21))
plot(loci,lwd=2)
points(xcoor,ycoor,col=4,lwd=2)
text(xcoor,ycoor,1:21,col=4,lwd=1,pos=2)

# Variogram models with the same "practical" range:
v.f <- function(x, ...){cov.spatial(x, ...)}
curve(v.f(x, cov.pars=c(1, 0.25)), from = 0, to = 1,
      xlab = "distance", ylab = expression(gamma(h)),
      main = "distance vs. correlation function")
      
# predicting by ordinary kriging
kc <- krige.conv(s100, loc=loci, krige=krige.control(cov.pars=c(1, .25)))

# -----------------------------------------------
# Plots: spatial observations and predictions 
# -----------------------------------------------

# Plots: "spatial observations" and "predictions" 
par(mfrow=c(1,2))
points(s100, cex.min=.3, cex.max=4,col=heat.colors(15))
image(kc, col=heat.colors(15))
title("kriging estimates")

# Plots:  "point estimates" and "standard errors"
par(mfrow=c(1,2))
image(kc, main="kriging estimates",col=heat.colors(15)); title("point estimates")
image(kc, val=sqrt(kc$krige.var),col=heat.colors(15)); title("kriging std. errors")

