tests/Attempt_New_Model_Mixed.R

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

library(reshape)
library(maps)
library(spatstat)
data(mesa.data.raw, package="SpatioTemporal")
AirQ.df = melt(mesa.data.raw$obs)
Xcov = mesa.data.raw$X[,7:12]
names(AirQ.df) <- c("date", "ID", "obs")
AirQ.df$date = as.Date(as.character(AirQ.df$date))
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 = "quarter", seasonality = c("quarter"), debug=F,method = "splines", degFree = 10)
df$basis<-NULL
df = updateSPTMData(df, "all")
y = df$obs.data
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, cov = NULL)
season.lu = updateSeason(df$season, p.sites = p.old)
resEM = estimEM(sptmod, theta.init = NULL, print.res = T, debug=F, hessian = T)
theta.EM = resEM$Result$theta
p.EM = resEM$Result$latent
plot(resEM, type= "residuals")
plot(resEM)



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))
sptmod = SPTModel(df, p.old, proj.xy, cov = Xcov[,1:4])
set.seed(1)
resGibbs = estimGibbs(sptmod, method = 'spat', priors = priors, print.res = F,
                      N.run = 1500, 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

plot(resGibbs, type= "residuals")
plot(resGibbs, type = "tempBasis")
plot(resGibbs, type= "prediction")

diagPlot(resGibbs, par = "theta", breaks = 50, main = "")
diagPlot(resGibbs, par = "beta")
diagPlot(resGibbs, par = "z")
diagPlot(resGibbs, par = "spatial")


df = SPTMData(AirQ.df, agg = "month", seasonality = c("quarter"), debug=F,method = "splines", degFree = 10)
sptmod = SPTModel(df, p.old, proj.xy, cov = NULL)
priors = list(theta = list(m0 =rep(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))
resGibbsM = estimGibbs(sptmod, method = 'mixedeffect', priors = priors, print.res = F,
                   N.run = 5000, beta.range = c(0,3), meanFunc = meanFunc, debug=F)
p.GibbsM = resGibbsM$Result$latent
colnames(p.GibbsM)<-colnames(p.old)
theta.GibbsM = resGibbsM$Result$theta

plot(resGibbsM, type= "residuals")
plot(resGibbsM, type = "tempBasis")
plot(resGibbsM, type= "prediction")

diagPlot(resGibbsM, par = "theta", breaks = 50, main = "")
diagPlot(resGibbsM, par = "random-effect", breaks = 20, main = "")
diagPlot(resGibbsM, par = "covariates", breaks = 20, main = "")
diagPlot(resGibbsM, par = "beta")
diagPlot(resGibbsM, par = "z")
diagPlot(resGibbsM, par = "spatial")