Nothing
predict.glmssn <- function(object, predpointsID, ...)
{
if(length(object$estimates$Warnlog) > 0 &&
length(grep("Algorithm diverging",object$estimates$Warnlog)) > 0)
stop("No predictions for diverging algorithm")
datao <- object$ssn.object@obspoints@SSNPoints[[1]]@point.data
ocoord <- object$ssn.object@obspoints@SSNPoints[[1]]@point.coords
theta <- object$estimates$theta
ind <- object$sampinfo$ind.obs
nobs <- length(ind)
CorModels <- object$args$CorModels
useTailDownWeight <-object$args$useTailDownWeight
REs <- object$sampinfo$REs
distord <- order(as.integer(as.character(datao[,"netID"])),
datao[,"pid"])
a.mat <- NULL
b.mat <- NULL
net.zero <- NULL
w.matrix <- NULL
dist.hydro <- NULL
xcoord <- NULL
ycoord <- NULL
xyobs <- NULL
xypred <- NULL
rnames <- NULL
if(predpointsID == "_MissingObs_")
{
if(length(grep("tail",CorModels)) > 0 | useTailDownWeight == TRUE){
if(length(grep("taildown",CorModels)) > 1)
stop("Cannot have more than 1 tailup model")
n.all <- length(datao[,1])
dist.junc <- matrix(0, nrow = nobs, ncol = nobs)
net.zero <- matrix(0, nrow = nobs, ncol = nobs)
nsofar <- 0
nIDs <- sort(as.integer(as.character(unique(datao[,"netID"]))))
for(i in nIDs){
workspace.name <- paste("dist.net", i, ".RData", sep="")
path <- file.path(object$ssn.object@path,
"distance", "obs", workspace.name)
file_handle <- file(path, "rb")
distmat <- unserialize(file_handle)
close(file_handle)
oi <- order(as.numeric(rownames(distmat)))
distmat <- distmat[oi, oi, drop = F]
rnames <- c(rnames,rownames(distmat))
ni <- length(distmat[1,])
dist.junc[(nsofar + 1):(nsofar + ni),
(nsofar + 1):(nsofar + ni)] <- distmat
net.zero[(nsofar + 1):(nsofar + ni),
(nsofar + 1):(nsofar + ni)] <- 1
nsofar <- nsofar + ni
}
rownames(dist.junc) <- rnames
# maximum distance to common junction between two sites
a.mat <- pmax(dist.junc,t(dist.junc))
# minimum distance to common junction between two sites
b.mat <- pmin(dist.junc,t(dist.junc))
# hydrological distance
dist.hydro <- as.matrix(dist.junc + t(dist.junc))*net.zero
flow.con.mat <- 1 - (b.mat > 0)*1
addfunccol <- object$args$addfunccol
w.matrix <- sqrt(pmin(outer(datao[distord,addfunccol],
rep(1, times = nobs)),
t(outer(datao[distord,addfunccol],rep(1, times = nobs)))) /
pmax(outer(datao[distord,addfunccol],rep(1, times = nobs)),
t(outer(datao[distord,addfunccol],rep(1, times = nobs))))) *
flow.con.mat*net.zero
if(length(grep("tailup",CorModels)) > 0){
if(length(grep("tailup",CorModels)) > 1)
stop("Cannot have more than 1 tailup model")
flow.con.mat <- 1 - (b.mat > 0)*1
}
}
xcoord <- ocoord[distord,1,drop = F]
ycoord <- ocoord[distord,2,drop = F]
REind <- which(names(datao) %in% CorModels)
if(length(REind)) {
REs <- list()
REnames <- sort(names(datao)[REind])
## model matrix for a RE factor
for(ii in 1:length(REind)) REs[[REnames[ii]]] <-
model.matrix(~datao[distord,REnames[ii]] - 1)
rownames(REs[[REnames[ii]]]) <- datao[distord,"pid"]
## corresponding block matrix
for(ii in 1:length(REind)) REs[[ii]] <- REs[[ii]] %*% t(REs[[ii]])
}
Vpred <- makeCovMat(theta = theta, dist.hydro = dist.hydro,
a.mat = a.mat, b.mat = b.mat, w.matrix = w.matrix,
net.zero = net.zero, x.row = xcoord, y.row = ycoord,
x.col = xcoord, y.col = ycoord,
CorModels, useTailDownWeight = useTailDownWeight,
use.nugget = object$args$use.nugget,
use.anisotropy = FALSE, REs)
for(i in 1:length(object$ssn.object@predpoints@SSNPoints))
if(object$ssn.object@predpoints@ID[i] == "_MissingObs_"){
datap <- object$ssn.object@predpoints@SSNPoints[[i]]@point.data
distordp <- order(as.integer(as.character(datap[,"netID"])),
datap[,"pid"])
datap <- datap[distordp, , drop = F]
pcoord <- object$ssn.object@predpoints@SSNPoints[[i]]@
point.coords
pcoord <- pcoord[distordp, , drop = F]
netpcoord <- object$ssn.object@predpoints@
SSNPoints[[i]]@network.point.coords
netpcoord <- netpcoord[distordp, , drop = F]
}
Vpred <- Vpred[!(datao[distord,"pid"] %in% datap[,"pid"]),
datao[distord,"pid"] %in% datap[,"pid"], drop = F]
}
if(predpointsID != "_MissingObs_")
{
#get the predicted data data.frame and coordinates from glmssn object
for(i in 1:length(object$ssn.object@predpoints@SSNPoints))
if(object$ssn.object@predpoints@ID[i] == predpointsID){
datap <- object$ssn.object@predpoints@SSNPoints[[i]]@point.data
pcoord <- object$ssn.object@predpoints@SSNPoints[[i]]@point.coords
netpcoord <- object$ssn.object@predpoints@SSNPoints[[i]]@network.point.coords
}
npred <- length(datap[,1])
distordp <- order(as.integer(as.character(datap[,"netID"])),
datap[,"pid"])
datap <- datap[distordp, , drop = F]
pcoord <- pcoord[distordp, , drop = F]
netpcoord <- netpcoord[distordp, , drop = F]
rnames <- NULL
cnames <- NULL
if(length(grep("tail",CorModels)) > 0){
if(length(grep("taildown",CorModels)) > 1)
stop("Cannot have more than 1 tailup model")
dist.junc.a <- matrix(0, nrow = nobs, ncol = npred)
dist.junc.b <- matrix(0, nrow = npred, ncol = nobs)
net.zero <- matrix(0, nrow = nobs, ncol = npred)
nSoFari <- 0
nSoFarj <- 0
nIDs <- sort(unique(c(as.integer(as.character(datao[,"netID"])),
as.integer(as.character(datap[,"netID"])))))
#loop through the networks to build covariance matrix
for(i in nIDs) {
ind.obs <- as.numeric(as.character(datao$netID)) ==
as.numeric(i)
site.no <- nrow(datao[ind.obs,])
ind.pred <- as.numeric(as.character(datap$netID)) ==
as.numeric(i)
pred.no <- nrow(datap[ind.pred,])
if(site.no*pred.no != 0) {
workspace.name.a <- paste("dist.net", i,
".a.RData", sep = "")
workspace.name.b <- paste("dist.net", i,
".b.RData", sep = "")
path.a <- file.path(object$ssn.object@path,
"distance", predpointsID, workspace.name.a)
path.b <- file.path(object$ssn.object@path,
"distance", predpointsID, workspace.name.b)
# distance matrix a
file_handle = file(path.a, open="rb")
distmata <- unserialize(file_handle)
close(file_handle)
# distance matrix b
file_handle = file(path.b, open="rb")
distmatb <- unserialize(file_handle)
close(file_handle)
# pid's for observed data in matrices above
ordoi <- order(as.numeric(rownames(distmata)))
ordpi <- order(as.numeric(rownames(distmatb)))
ni <- length(ordoi)
nj <- length(ordpi)
distmata <- distmata[ordoi, ordpi, drop = F]
distmatb <- distmatb[ordpi, ordoi, drop = F]
if(any(rownames(distmata)!=colnames(distmatb)))
stop("rownames of distmata do not match colnames of distmatb")
if(any(colnames(distmata)!=rownames(distmatb)))
stop("colnames of distmata do not match rownames of distmatb")
dist.junc.a[(nSoFari + 1):(nSoFari + ni),
(nSoFarj + 1):(nSoFarj + nj)] <- distmata
dist.junc.b[(nSoFarj + 1):(nSoFarj + nj),
(nSoFari + 1):(nSoFari + ni)] <- distmatb
net.zero[(nSoFari + 1):(nSoFari + ni),
(nSoFarj + 1):(nSoFarj + nj)] <- 1
}
else
{
ni <- site.no
nj <- pred.no
}
nSoFari <- nSoFari + ni
nSoFarj <- nSoFarj + nj
}
# rownames are for observed data, colnames are prediction data
rnames <- datao[distord,"pid"]
cnames <- datap[,"pid"]
rownames(dist.junc.a) <- rnames
colnames(dist.junc.a) <- cnames
# creat A matrix (longest distance to junction of two points)
a.mat <- pmax(dist.junc.a,t(dist.junc.b))
# creat B matrix (shorted distance to junction of two points)
b.mat <- pmin(dist.junc.a,t(dist.junc.b))
# get hydrologic distance
dist.hydro <- as.matrix(dist.junc.a + t(dist.junc.b))
# create indicator matrix of flow-connected
if(length(grep("tailup",CorModels)) > 0){
if(length(grep("tailup",CorModels)) > 1)
stop("Cannot have more than 1 tailup model")
flow.con.mat <- 1 - (b.mat > 0)*1
# weight matrix based on additive function
addfunccol <- object$args$addfunccol
w.matrix <- sqrt(pmin(outer(datao[distord,addfunccol],
rep(1, times = npred)),
t(outer(datap[,addfunccol],rep(1, times = nobs)))) /
pmax(outer(datao[distord,addfunccol],rep(1, times = npred)),
t(outer(datap[,addfunccol],rep(1, times = nobs))))) *
flow.con.mat*net.zero
if(any(rownames(w.matrix)!=rownames(a.mat)))
stop("rownames of w.matrix do not match rownames of a.mat")
if(any(colnames(w.matrix)!=colnames(a.mat)))
stop("colnames of w.matrix do not match colnames of a.mat")
}
}
# create Euclidean matrix among observed data
xyobs <- ocoord
x.samp <- ocoord[distord,1,drop = F]
y.samp <- ocoord[distord,2,drop = F]
xypred <- pcoord
x.pred <- pcoord[, 1, drop = F]
y.pred <- pcoord[, 2, drop = F]
if(any(rownames(x.samp)!=rownames(a.mat)))
stop("rownames of x.samp do not match rownames of a.mat")
if(any(rownames(x.pred)!=colnames(a.mat)))
stop("rownames of x.pred do not match colnames of a.mat")
REPs <- NULL
if(!is.null(REs)) {
REnames <- names(REs)
for(ii in 1:length(REnames)) if(any(is.na(datap[,REnames[ii]])))
stop("Cannot having missing values when creating random effects")
REOs <- list()
REPs <- list()
ObsSimDF <- datao
PredSimDF <- datap
## model matrix for a RE factor
for(ii in 1:length(REnames)){
#we'll add "o" to observed levels and "p" to prediction
# levels so create all possible levels
plevels <- unique(c(levels(PredSimDF[,REnames[[ii]]]),
paste("o",levels(ObsSimDF[,REnames[[ii]]]),sep = ""),
paste("p",levels(PredSimDF[,REnames[[ii]]]),sep = "")))
# sites with prediction levels same as observation levels
pino <- PredSimDF[,REnames[[ii]]] %in% ObsSimDF[,REnames[[ii]]]
#add "o" to observed levels
ObsSimDF[,REnames[[ii]]] <- paste("o",
ObsSimDF[,REnames[[ii]]], sep = "")
ObsSimDF[,REnames[[ii]]] <- as.factor(as.character(
ObsSimDF[,REnames[[ii]]]))
#add all possible levels to prediction data frame
levels(PredSimDF[,REnames[[ii]]]) <- plevels
# add "o" to prediction sites with observation levels
if(any(pino)) PredSimDF[pino,REnames[[ii]]] <- paste("o",
PredSimDF[pino,REnames[[ii]]], sep = "")
# add "p" to all predicition sites without observation levels
if(any(!pino)) PredSimDF[!pino,REnames[[ii]]] <- paste("p",
PredSimDF[!pino,REnames[[ii]]], sep = "")
PredSimDF[,REnames[[ii]]] <- as.factor(as.character(
PredSimDF[,REnames[[ii]]]))
# now get down to just levels with "o" & "p" added
blevels <- unique(c(levels(ObsSimDF[,REnames[[ii]]]),
levels(PredSimDF[,REnames[[ii]]])))
ObsSimDF[,REnames[[ii]]] <- factor(ObsSimDF[,REnames[[ii]]],
levels = blevels, ordered = FALSE)
PredSimDF[,REnames[[ii]]] <- factor(PredSimDF[,REnames[[ii]]],
levels = blevels, ordered = FALSE)
# now ordering of factors in Z matrices should be compatible
# with obs x obs Z matrices
REOs[[ii]] <- model.matrix(~ObsSimDF[distord,
REnames[[ii]]] - 1)
REPs[[ii]] <- model.matrix(~PredSimDF[,
REnames[[ii]]] - 1)
rownames(REOs[[ii]]) <- datao[distord,"pid"]
rownames(REPs[[ii]]) <- datap[,"pid"]
if(any(rownames(REOs[[ii]])!=rownames(a.mat)))
stop("rownames RE for obs do not match rownames of a.mat")
if(any(rownames(REPs[[ii]])!=colnames(a.mat)))
stop("rownames RE for preds do not match colnames of a.mat")
}
## corresponding block matrix
for(ii in 1:length(REnames)) REPs[[ii]] <-
REOs[[ii]] %*% t(REPs[[ii]])
}
Vpred <- makeCovMat(theta = theta, dist.hydro = dist.hydro,
a.mat = a.mat, b.mat = b.mat, w.matrix = w.matrix,
net.zero = net.zero, x.row = x.samp, y.row = y.samp,
x.col = x.pred, y.col = y.pred,
CorModels = CorModels,
useTailDownWeight = useTailDownWeight,
use.nugget = FALSE,
use.anisotropy = object$args$use.anisotropy, REs = REPs)
Vpred <- Vpred[ind, , drop = F]
}
# get a list of response and covariate names
response.col <- object$args$zcol
mod.names <- as.character(attr(terms(object$args$formula,
data = object$ssn.object@obspoints@SSNPoints[[1]]@point.data),"variables"))
# get the number of names ( + 1, as the first is always "list")
nc.tmp <- length(mod.names)
# if there are any covariates ...
ind.allcov <- rep(TRUE, times = length(datap[,1]))
if(nc.tmp > 2) {
# create a FALSE for a record with missing values of the covariates
for(i in 3:nc.tmp) ind.allcov <- ind.allcov & !is.na(datap[,mod.names[i]])
}
# prediction sample size without missing covariates
np.allcov <- sum(ind.allcov)
# add response variable as -1 in datap
datap[,response.col] <- NA
# add prediction standard errors as NAs in datap
datap[,paste(response.col,".predSE", sep = "")] <- NA
# remove records that had any missing values for any of the covariates
datap1 <- datap[ind.allcov,]
# add response variable as -1 in datap
datap1[,response.col] <- -1
# add prediction standard errors as NAs in datap
datap1[,paste(response.col,".predSE", sep = "")] <- NA
formula <- object$args$formula
# create design matrix for prediction data set
mf <- model.frame(formula, data = datap1)
mt <- attr(mf, "terms")
Xpred <- model.matrix(mt, mf)
Xpred <- Xpred[,object$sampinfo$cutX1toX2, drop = F]
# get the sum of partial sills
sumparsil <- sum(theta[attr(theta,"type") == "parsill"])
Vi <- object$estimates$Vi
covb <- object$estimates$covb
Xobs <- object$sampinfo$X
z <- object$sampinfo$z
n <- object$sampinfo$obs.sample.size
p <- object$sampinfo$rankX
parsilvec <- rep(sumparsil, times = length(Vpred[1,]))
if(object$args$family == "poisson") {
beta.hat <- object$estimates$betahat
eta.hatp <- Xpred[ind.allcov,] %*% beta.hat
eta.hato <- Xobs %*% beta.hat
#diagonal elements of Delta~^{-1} of my manuscript
Del.ip <- as.vector(1/exp(eta.hatp))
Del.io <- as.vector(1/exp(eta.hato))
#diagonal elements of A^(1/2) of my manuscript
A.5p <- as.vector(sqrt(exp(eta.hatp)))
A.5o <- as.vector(sqrt(exp(eta.hato)))
Vpred <- t((Del.ip*A.5p)*t(Del.io*A.5o*Vpred))
parsilvec <- sumparsil*(A.5p*Del.ip)^2
}
if(object$args$family == "binomial") {
beta.hat <- object$estimates$betahat
eta.hatp <- Xpred[ind.allcov,] %*% beta.hat
eta.hato <- Xobs %*% beta.hat
#diagonal elements of Delta~^{-1} of my manuscript
Del.ip <- as.vector((1 + exp(eta.hatp))^2/exp(eta.hatp))
Del.io <- as.vector((1 + exp(eta.hato))^2/exp(eta.hato))
#diagonal elements of A^(1/2) of my manuscript
A.5p <- as.vector(sqrt(exp(eta.hatp)/(1 +
exp(eta.hatp))^2))
A.5o <- as.vector(sqrt(exp(eta.hato)/(1 +
exp(eta.hato))^2/object$sampinfo$trialsvec))
Vpred <- t((Del.ip*A.5p)*t(Del.io*A.5o*Vpred))
parsilvec <- sumparsil*(A.5p*Del.ip)^2
}
M <- rbind(Vpred, t(Xpred), parsilvec)
XXSiXi <- Xobs %*% covb
XSi <- t(Xobs) %*% Vi
pred.out <- t(apply(M, 2, UK4Apply, covb = covb,
XXSiXi = XXSiXi, XSi = XSi, Vi = Vi, z = z, n = n, p = p))
datap1[,response.col] <- pred.out[,1]
datap1[,paste(response.col,".predSE", sep = "")] <- pred.out[,2]
datap[ind.allcov,] <- datap1
#put the predictions in the predicted data data.frame and return the SSN object
for(i in 1:length((object$ssn.object@predpoints@SSNPoints)))
if(object$ssn.object@predpoints@ID[i] == predpointsID){
object$ssn.object@predpoints@SSNPoints[[i]]@point.data <-
datap[order(distordp),]
}
object$args$predpointsID <- predpointsID
class(object) <- "glmssn.predict"
object
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.