vignettes/SpTMixture_OtherData.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

## ---- message=FALSE, warning=FALSE---------------------------------------
library(SpTMixture, quiet = TRUE)
require("rgdal", quiet = TRUE) # requires sp, will use proj.4 if installed
require("maptools", quiet = TRUE)
require("ggplot2", quiet = TRUE)
require("plyr", quiet = TRUE)
library(grid, quiet = TRUE)
library(gridExtra, quiet = TRUE)
library(mytoolbox, quiet = TRUE)
library(raster, quiet=TRUE)
library(reshape)
library(maps)

## ------------------------------------------------------------------------
data(mesa.data.raw, package="SpatioTemporal")
AirQ.df = melt(mesa.data.raw$obs)
names(AirQ.df) <- c("date", "ID", "obs")
AirQ.df$date = as.Date(as.character(AirQ.df$date))

## ------------------------------------------------------------------------
ggplot(data = AirQ.df, aes(x = date, y = obs, color = ID, group = ID)) + 
  geom_point() + ylab("log(obs)") + geom_line() + facet_wrap(~ID)

## ------------------------------------------------------------------------
subCoords.df = mesa.data.raw$X[,c(1,2,3,4,5)]

## ------------------------------------------------------------------------
m <- map("county", "california", plot =F, fill = TRUE)
IDs <- sapply(strsplit(m$names, ":"), function(x) x[1])

Cal.poly <- map2SpatialPolygons(m, IDs=IDs, proj4string=CRS("+proj=longlat +ellps=GRS80 +no_defs"))

bbox = t(apply(subCoords.df[,4:5], 2,function(x) c(min(x), max(x))))
rownames(bbox) <- c("lon","lat")
colnames(bbox) <- c('min','max')
polybox = asSpatialPolygonsBbox(bbox, zoom = 1.2)
polybox <- spTransform(polybox, CRS(proj4string(Cal.poly)))
CalMap <- clipToExtent(Cal.poly, polybox)
Cal.df = fortify(CalMap, region="id")


## ------------------------------------------------------------------------
nLu = 5
p.old = matrix(runif(nLu*nrow(subCoords.df),0,1), ncol = nLu, nrow = nrow(subCoords.df))
p.old= as.matrix((p.old)/rowSums(p.old))
colnames(p.old) <- paste("Cl",1:nLu)

## ------------------------------------------------------------------------
df = SPTMData(AirQ.df, agg = "month", seasonality = c("quarter"), debug=F,
               method = "splines", degFree = 10)
df$basis<-NULL
df = updateSPTMData(df, "all")
y = df$obs.data

## ---- fig.height=6-------------------------------------------------------
season.sites=df$season
season.lu = updateSeason(season.sites, p.sites = p.old)
plot(season.lu, select = 1:nLu, cex.axis = 0.7)


## ------------------------------------------------------------------------
theta.mean = rnorm(ncol(p.old), mean = rep(mean(y$obs, na.rm=T),ncol(p.old)),
                  sd = rep(mean(abs(y$obs), na.rm=T),ncol(p.old))/10)
theta.season = rnorm(max(ncol(df$basis),2)*ncol(p.old), mean = 1, sd = 0.5)
theta.sigma = rnorm(ncol(p.old), mean = -0.5, sd = 0.05)
theta.0 = c(theta.mean,theta.season, theta.sigma)
subCoords.df = subCoords.df[match(levels(df$obs.data$ID),subCoords.df$ID),]
proj.xy = cbind(subCoords.df$x,subCoords.df$y)
subCoords.df$x = proj.xy[,1]
subCoords.df$y = proj.xy[,2]
sptmod = SPTModel(df, p.old, proj.xy)

## ------------------------------------------------------------------------
resEM = estimEM(sptmod, theta.init = NULL, print.res = T, debug=F, hessian = F)
theta.EM = resEM$Result$theta
p.EM = resEM$Result$latent

## ---- fig.height=8-------------------------------------------------------
plot(resEM, type = "residuals")
plot(resEM, type = "tempBasis")
# plot(resEM, type= "prediction")

## ------------------------------------------------------------------------
# colnames(p.EM)<-colnames(p.old)
# season.lu = updateSeason(season.sites, p.sites=p.EM)
# season.lu[is.na(season.lu)] <- NA
# dataS = data.frame(
#   Comp = as.factor(sapply(apply(p.EM, 1, which.max), function(x) which(sort(theta.EM[1:ncol(p.old)],index.return=T)$ix == x))),
#   Longitude = subCoords.df[,4], 
#   Latitude = subCoords.df[,5]
#   )
# 
# SpPt = SpatialPointsDataFrame(coords = as.matrix(subCoords.df[,4:5]), data = dataS, proj4string = CRS("+proj=longlat +ellps=GRS80 +no_defs"))

## ---- fig.height=8-------------------------------------------------------
# ggplot() + geom_polygon(data = Cal.df, aes(long,lat, group = group)) + 
#   geom_path(data = Cal.df, aes(long,lat, group = group),color="darkgrey", size = 1) +
#   coord_equal() + 
#   geom_point(data = SpPt@data, aes(x = Longitude, y = Latitude, col =Comp), size=3) + 
#   scale_color_manual(values = rev(brewer.pal(ncol(p.EM),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.EM)[sort(theta.EM[1:ncol(p.old)],index.return=T)$ix])+ 
#   theme(legend.position="bottom")

## ---- fig.height=8-------------------------------------------------------
# Ypred = predict(resEM, transform = "none")
# Ypred$ClName = as.factor(colnames(p.EM)[as.numeric(as.character(Ypred$Cluster))])
# Ypred$ClName <- ordered(Ypred$ClName, levels = colnames(p.EM)[sort(theta.EM[1:ncol(p.old)],index.return=T)$ix])
# ggplot(data = Ypred) + geom_point(aes(x=date, y = obs), size = 0.4)+
#   geom_line(aes(x = date, y = pred, col = ClName), size=0.5)+
#   geom_ribbon(aes(x = date, ymin=CIlow,ymax=CIup, col = ClName),alpha=0.3)+
#     facet_wrap(~ID, scales = "free") + 
#     scale_color_manual(values = rev(brewer.pal(ncol(p.EM),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.EM)[sort(theta.EM[1:ncol(p.old)],index.return=T)$ix]) + 
#   scale_fill_manual(values = rev(brewer.pal(ncol(p.EM),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.EM)[sort(theta.EM[1:ncol(p.old)],index.return=T)$ix]) + 
#   theme(legend.position="bottom")


## ------------------------------------------------------------------------
# theta.EM[(length(theta.EM) - ncol(p.EM) + 1):length(theta.EM)] = exp(theta.EM[(length(theta.EM) - ncol(p.EM) + 1):length(theta.EM)])
# 
# fisher_info<-solve(-resEM$Result$Hessian)
# prop_theta<-sqrt((diag(fisher_info)))
# 
# idx.s = (ncol(season.lu)+1):length(theta.EM)
# priors = list(theta = list(m0 = theta.EM[1:(ncol(season.lu))], s0 = prop_theta[1:ncol(season.lu)]),
#               sigma2 = list(a0 = theta.EM[idx.s]^2 / ((prop_theta[idx.s])) + 2, b0 = (theta.EM[idx.s]^2 / ((prop_theta[idx.s])) + 1)*theta.EM[idx.s]),
#                z = list(a = p.EM)) 
# set.seed(1)
# resGibbs = estimGibbs(sptmod, method = 'simple', priors = priors, print.res = F,
#                     N.run = 250, beta.range = c(0,3), meanFunc = meanFunc, debug=F)
# p.Gibbs = resGibbs$Result$latent
# colnames(p.Gibbs)<-colnames(p.old)
# theta.Gibbs = resGibbs$Result$theta

## ---- fig.height=8-------------------------------------------------------
# plot(resGibbs, type = "residuals")
# plot(resGibbs, type = "tempBasis")


## ---- fig.height=8-------------------------------------------------------
# diagPlot(resGibbs, par = "theta",main = "", breaks = 200, plot.cluster = 1:5)
# diagPlot(resGibbs, par = "z",main = "")
# diagPlot(resGibbs, par = "paired-theta", plot.cluster = 1)

## ---- fig.height=8-------------------------------------------------------
# Ypred = predict(resGibbs)
# Ypred$ClName = as.factor(colnames(p.old)[as.numeric(as.character(Ypred$Cluster))])
# Ypred$ClName <- ordered(Ypred$ClName, levels = colnames(p.old)[sort(theta.Gibbs[1:ncol(p.old)],index.return=T)$ix])
# ggplot(data = Ypred) + geom_point(aes(x=date, y = obs), size = 0.4)+
#   geom_line(aes(x = date, y = pred, color = (ClName)))+
#   geom_ribbon(aes(x = date, ymin=CIlow,ymax=CIup, fill = ClName),alpha=0.3)+
#     facet_wrap(~ID, scales = "free") + 
#     scale_color_manual(values = rev(brewer.pal(ncol(p.Gibbs),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.old)[sort(theta.Gibbs[1:ncol(p.old)],index.return=T)$ix]) + 
#    scale_fill_manual(values = rev(brewer.pal(ncol(p.Gibbs),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.old)[sort(theta.Gibbs[1:ncol(p.old)],index.return=T)$ix]) + 
#   theme(legend.position="bottom")


## ------------------------------------------------------------------------
# colnames(p.Gibbs)<-colnames(p.old)
# season.lu = updateSeason(season.sites, p.sites=p.Gibbs)
# season.lu[is.na(season.lu)] <- NA
# dataS = data.frame(
#   Comp = as.factor(sapply(apply(p.Gibbs, 1, which.max), function(x) which(sort(theta.Gibbs[1:ncol(p.old)],index.return=T)$ix == x))),
#   Longitude = subCoords.df[,4], 
#   Latitude = subCoords.df[,5]
#   )
# 
# SpPt = SpatialPointsDataFrame(coords = as.matrix(subCoords.df[,4:5]), data = dataS, proj4string = CRS("+proj=longlat +ellps=GRS80 +no_defs"))

## ---- fig.height=8-------------------------------------------------------
# ggplot() + geom_polygon(data = Cal.df, aes(long,lat, group = group)) + 
#   geom_path(data = Cal.df, aes(long,lat, group = group),color="darkgrey", size = 1) +
#   coord_equal() + 
#   geom_point(data = SpPt@data, aes(x = Longitude, y = Latitude, col =Comp), size=3) + 
#   scale_color_manual(values = rev(brewer.pal(ncol(p.Gibbs),"RdYlGn")), name = "LU Composition",
#                      labels = colnames(p.Gibbs)[sort(theta.Gibbs[1:ncol(p.old)],index.return=T)$ix])+ 
#   theme(legend.position="bottom")