tests/Attempt_New_Model.R

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

library(SpTMixture)
library(raster)
library(rgeos)
library(maptools)

data(GoyderData)
subWaterQ.df = subset(WaterQ.df, param == levels(WaterQ.df$param)[6],
                      select = c(site, date, value))
subWaterQ.df = subset(subWaterQ.df, date < "2009-01-01", select = c(site, date, value))
subWaterQ.df = subset(subWaterQ.df, date > "2000-01-01", select = c(site, date, value))
ID_keep = names(which(table(subWaterQ.df$site)>100))
subWaterQ.df = subset(subWaterQ.df, site %in% ID_keep )
names(subWaterQ.df) <- c("ID" ,"date", "obs")
subWaterQ.df$ID <- as.factor(as.character(subWaterQ.df$ID))
subWaterQ.df$obs[subWaterQ.df$obs==0] <- NA
subWaterQ.df$obs <- log(subWaterQ.df$obs)

# LU = area2comp(LandUse.df, compCol = "DES_PRI_V6", IDcol = "SiteID")
# subLandUse.df = subset(LU, ID %in% unique(subWaterQ.df$ID))
# rownames(subLandUse.df) <- subLandUse.df$ID
# subLandUse.df$ID <- NULL
# keep_lu = names(which(apply(subLandUse.df,2,max) > .487))
# # Merging GMP and OMU as they seem to behave identically
# subLandUse.df = subset(subLandUse.df, select = keep_lu)
# subWaterQ.df = subset(subWaterQ.df, ID %in% rownames(subLandUse.df))
# subWaterQ.df$ID <- as.factor(as.character(subWaterQ.df$ID))
# subCoords.df = unique(subset(Coords.df, site %in% rownames(subLandUse.df)))
# names(subCoords.df)[1] <- c("ID")

# bbox = t(apply(subCoords.df[,2:3], 2,function(x) c(min(x), max(x))))
# rownames(bbox) <- c("lon","lat")
# colnames(bbox) <- c('min','max')
# polybox = asSpatialPolygonsBbox(bbox, zoom = 0.85)
# proj4string(polybox) <- proj4string(CatchmentMap)
# CatchmentMapC <- crop(CatchmentMap,polybox, byid = TRUE, drop_lower_td = TRUE)
# CatchmentMapC@data$id = CatchmentMapC@data$NewCatchId
# Catchment.points = fortify(CatchmentMapC, region="id")
# Catchment.df = join(Catchment.points,CatchmentMapC@data, by="id")
# Catchment.df$SiteID = FALSE
# Catchment.df$SiteID[Catchment.df$Site_Code %in% levels(subWaterQ.df$ID)] = TRUE
# Catchment.df$SiteID = as.factor(as.character(Catchment.df$SiteID))

# p.old = subLandUse.df[,which(colSums(subLandUse.df) > 0)]
# p.old[is.na(p.old)] = 0
# p.old= as.matrix((p.old)/rowSums(p.old))
df = SPTMData(subWaterQ.df, agg = "week", seasonality = c("month"), debug=F,method = "splines", degFree = 7)

# p.old = p.old[match(levels(df$obs.data$ID),rownames(p.old)),]
y = df$obs.data

bF0 = df$basis[,df$list.idx[[2]]]
bF1 = df$basis[,df$list.idx[[3]]]

theta = c(runif(2 + 2*ncol(bF0) + ncol(bF1),0,1),0)

par(mfrow = c(1,1), mar = c(2,2,2,2))
for(i in 1:8){
Y = y[y$ID == unique(y$ID)[i],'obs']
X = y[y$ID == unique(y$ID)[i],'date']
n.obs = length(which(!is.na(Y)))
if(n.obs > 50){
conv = 1
thetaI = theta
while(conv == 1){
res = optim(par = thetaI, fn = ll, y = Y, bF0 = bF0, bF1=bF1, control = list(maxit = 20000))
thetaI = res$par
conv = res$convergence
}
Yp = mF(theta=res$par, bF0 = bF0, bF1 = bF1)
CI = (cbind(Yp$full -1.96*(exp(res$par[length(res$par)])), Yp$full +1.96*(exp(res$par[length(res$par)]))))
ylim = range((CI))
if(is.na(max(CI))){ylim = range(Y, na.rm=T)}
plot(X , (Y), ylim = ylim, xlim = range(X[which(!is.na(Y))]), cex = 0.5)
polygon(c(X, rev(X)), c(CI[,1], rev(CI[,2])),border = NA, col = "lightgrey", density = 70)
points(X,(Y), cex=.5)
points(X,(Yp$full), col = "blue", type="l")
points(X,(Yp$trend), col = "orange", type="l", lwd = 2)
# points(Yp$seasonality, type= "l")
# points(Yp$trendS, type= "l")
hist(exp(Y) - exp(Yp$full),20, main = "")
}
}
####
f11 = function(x){atan(x)*2/pi}
f01 = function(x){atan(x)/pi + 0.5}

mF = function(theta,site = 1, cluster = 1, nsites =1, ncluster=1, bF0, bF1){

    idx_theta_Site = (2*site-1):(2*site)
    idx_theta_cluster0 = 2*nsites + (cluster-1)*(ncol(bF0) + ncol(bF1)) + 1:ncol(bF0)
    idx_theta_cluster00 = max(idx_theta_cluster0) + 1:ncol(bF0)
    idx_theta_cluster1 = max(idx_theta_cluster00) + 1:ncol(bF1)

    beta0 = theta[idx_theta_Site[1]]
    beta1 = theta[idx_theta_Site[2]]
    nu = matrix(f11(theta[idx_theta_cluster0]), ncol = 1)
    rho = matrix(f11(theta[idx_theta_cluster1]), ncol = 1)
    mu = matrix(f01(theta[idx_theta_cluster00]), ncol = 1)

    season =  (bF1 %*% rho)
    trendS = (1 + bF0 %*% mu)
    trend = beta0 * (1 + bF0 %*% nu)
    yT = trend + beta1 * season * trendS
    return(list(full = yT, seasonality = season, trend = trend, trendS = trendS ))
}

ll = function(theta, y,site = 1, cluster = 1, nsites =1, ncluster=1, bF0, bF1){
  -sum(dnorm(y - mF(theta,site, cluster, nsites, ncluster, bF0, bF1)$full,0,exp(theta[length(theta)]), log = T), na.rm=T)
}