Nothing
R/useful_commands.R
In spTDyn: Spatially Varying and Spatio-Temporal Dynamic Linear Models
Defines functions
prob.below.threshold
reverse.truncated.fnc
truncated.fnc
ObsGridData
gridTodata
ObsGridLoc
tp
sp
fnc.time
spT.check.sites.inside
spT.check.locations
PMCC
spT.keep.morethan.dist
spT.hit.false
spT.segment.plot
spT.Summary.Stat
tp.dimname.fnc
sp.dimname.fnc
plot.spTD
summary.spTD
as.mcmc.spTD
residuals.spTD
coef.spTD
confint.spTD
fitted.spTD
print.spTD
Documented in
as.mcmc.spTD
coef.spTD
confint.spTD
fitted.spTD
fnc.time
gridTodata
ObsGridData
ObsGridLoc
plot.spTD
PMCC
print.spTD
prob.below.threshold
residuals.spTD
reverse.truncated.fnc
sp
sp.dimname.fnc
spT.check.locations
spT.check.sites.inside
spT.hit.false
spT.keep.morethan.dist
spT.segment.plot
spT.Summary.Stat
summary.spTD
tp
tp.dimname.fnc
truncated.fnc
##
## Print function
##
print.spTD<-function(x, ...) {
cat("-----------------------------------------------------"); cat('\n');
cat("Model: "); cat(x$model); cat('\n');
cat("Call: "); print(x$call); #cat('\n')
cat("Iterations: "); cat(x$iterations); cat("\n")
cat("nBurn: "); cat(x$nBurn); cat("\n")
cat("Acceptance rate for phi (%): "); cat(x$accept); cat("\n")
cat("-----------------------------------------------------"); cat('\n');
#cat("PMCC: "); cat("\n");
print(x$PMCC);
cat("-----------------------------------------------------"); cat('\n');
#cat("Parameters:\n")
#print(x$parameter); #cat("\n");
#cat("-----------------------------------------------------"); cat('\n');
cat("Computation time: "); cat(x$computation.time); cat("\n")
}
##
## fitted function
##
fitted.spTD<-function(object, ...){
x<-data.frame(object$fitted)
x
}
##
## use of confint()
##
confint.spTD<-function(object, parm, level=0.95, ...){
#
x<-as.mcmc(object)
up<-level+(1-level)/2
low<-(1-level)/2
FUN <- function(x){quantile(x,probs=c(low,up))}
out<-apply(x,2,FUN=FUN)
out<-t(out)
if(missing(parm)){
out
}
else{
if(length(parm)>1){
out<-as.matrix(out[dimnames(out)[[1]] %in% parm,])
out
}
else{
out<-t(as.matrix(out[dimnames(out)[[1]] %in% parm,]))
dimnames(out)[[1]]<-parm
out
}
}
}
##
## use of coefficients
##
coef.spTD<-function(object, digits=4, ...){
round(t(object$parameter)[1,],digits=digits)
}
##
## use of residuals
##
residuals.spTD<-function(object, ...){
#if(object$combined.fit.pred==TRUE){
# stop("\n# Error: not useful for output with combined fit and predict")
#}
#else{
if(object$scale.transform=="NONE"){
tmp<-object$Y-object$fitted[,1]
tmp
}
else if(object$scale.transform=="SQRT"){
tmp<-sqrt(object$Y)-object$fitted[,1]
tmp
}
else if(object$scale.transform=="LOG"){
tmp<-log(object$Y)-object$fitted[,1]
tmp
}
else{
}
#}
}
##
## To use in coda package
##
as.mcmc.spTD<-function(x, ...){
if (x$combined.fit.pred == TRUE) {
stop("\n# Error: not available for combined fit-prediction option. Please read the MCMC samples from the *.txt file manually.")
}
model <- x$model
if (is.null(model) == TRUE) {
stop("\n# Error: need to define the model")
}
else if (model == "GP" | model == "truncatedGP") {
r <- x$r
p <- x$p
#
if(x$cov.fnc=="matern"){
if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi", "nu")
}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi", "nu")
}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi", "nu")
}
else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi", "nu")
}
else {
stop("Error")
}
}
else {
if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi")
}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi")
}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi")
}
else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi")
}
else {
stop("Error")
}
}
#
para<-t(para)
para<-mcmc(para)
para
}
else {
}
}
##
## use of summary
##
summary.spTD<-function(object, digits=4, package="spTDyn", coefficient=NULL, ...){
#
if(package=="coda"){
if(object$combined.fit.pred==TRUE){
stop("\n# Error: coda package is not useful for output with combined fit and predict \n")
}
else{
if(is.null(coefficient)){
cat("\n#### MCMC summary statistics using coda package ####\n")
if(!is.null(object$sp.covariate.names)) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
tmp<-as.mcmc(object)
summary(tmp, ...)
}
else if(coefficient=="spatial"){
if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betasp))
if(object$model=="GPP"){n<-object$knots}
else{n<-object$n}
tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q)
summary(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the temporally varying dynamic coefficients ####\n")
tmp<-as.mcmc(t(object$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T)
summary(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
tmp<-as.mcmc(t(object$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="")
summary(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
}
}
}
else{
if(package=="spTDyn"){
coefficient <- coefficient
}
#
else if("Xsp"%in%names(object) | "Xtp"%in%names(object)){
coefficient <- coefficient
}
else{
coefficient <- NULL
}
#
if(is.null(coefficient)){
if((!is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
else if((!is.null(object$sp.covariate.names)) & (is.null(object$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
else if((is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
else { }
#ans<-NULL
#ans$Model=object$model;ans$Call=object$call;ans$Iterations=object$iterations;ans$nBurn=object$nBurn;ans$PMCC=object$pmcc;ans$Parameters=round(object$parameter,digits=digits)
#ans
print(object)
cat("-----------------------------------------------------"); cat('\n');
cat("Parameters:\n")
print(round(object$parameter,digits=digits)); #cat("\n");
cat("-----------------------------------------------------"); cat('\n');
}
else if(coefficient=="spatial"){
if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betasp))
if(object$model=="GPP"){n<-object$knots}
else{n<-object$n}
tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q)
summary(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T)
summary(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
tmp<-as.mcmc(t(object$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="")
summary(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes character 'spatial' and 'temporal'.")
}
}
}
##
## use of plot
##
plot.spTD<-function(x, residuals=FALSE, coefficient=NULL, ...){
if(as.logical(residuals)==FALSE){
if(x$combined.fit.pred==TRUE){
if(!is.null(x$sp.covariate.names)) {cat("## \n# Spatially varying parameters are not included.\n##\n ")}
cat("\n## Only residual plots are available for output with combined fit and predict option \n\n")
plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
par(ask=TRUE)
qqnorm(residuals(x));qqline(residuals(x),lty=2)
}
else{
if(is.null(coefficient)){
if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
else { }
if(x$model=="GP"){
tmp<-as.mcmc(x)
plot(tmp, ...)
}
else{
x$model="GP"
tmp<-as.mcmc(x)
plot(tmp, ...)
}
}
else if(coefficient=="spatial"){
if(is.null(x$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
tmp<-as.mcmc(t(x$betasp))
if(x$model=="GPP"){n<-x$knots}
else{n<-x$n}
tmp<-sp.dimname.fnc(x=tmp,names=x$sp.covariate.names,n=n,q=x$q)
plot(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(x$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
tmp<-as.mcmc(t(x$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=x$tp.covariate.names,u=x$u,T=x$T)
plot(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(x$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
tmp<-as.mcmc(t(x$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:x$u,sep="")
plot(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
}
}
}
else{
if(!is.null(x$sp.covariate.names)) {cat("## \n# Models fitted with spatially varying parameters.\n## \n ")}
plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
par(ask=TRUE)
qqnorm(residuals(x));qqline(residuals(x),lty=2)
} # end of 2nd if
}
##
## to include dimnames: function used in function summary statistics for spatially varying coef
##
sp.dimname.fnc<-function(x,names,n,q){
dimnames(x)[[2]][1:(n*q)]<-1:(n*q)
for(i in 1:q){
dimnames(x)[[2]][(1+(i-1)*n):(n*i)]<-rep(paste(names[i],"site",1:n,sep=""))
}
x
}
##
## to include dimnames: function used in function summary statistics for temporally varying coef
##
tp.dimname.fnc<-function(x,names,u,T){
dimnames(x)[[2]][1:(u*T)]<-1:(u*T)
for(i in 1:u){
dimnames(x)[[2]][(1+(i-1)*T):(T*i)]<-rep(paste(names[i],"time",1:T,sep=""))
}
x
}
##
## Summary Statistics
##
spT.Summary.Stat <- function(y)
{
## define Upper and lower limit function
up.low.limit<-function(y,limit)
{
#
y<-sort(y)
y<-y[limit]
y
#
}
#
if(is.vector(y)==TRUE){
y <- matrix(y[!is.na(y)])
}
else{
y <- t(y)
}
N <- length(y[1, ])
nItr <- length(y[, 1])
z <- matrix(nrow = N, ncol = 5)
dimnames(z) <- list(dimnames(y)[[2]], c("Mean","Median","SD","Low2.5p","Up97.5p"))
#
if (nItr < 40) {
stop("\n##\n# Error: number of samples must be >= 40\n##")
}
z[, 1] <- apply(y,2,mean)
z[, 2] <- apply(y,2,median)
z[, 3] <- apply(y,2,sd)
#z[, 4] <- apply(y,2,quantile,0.025)
#z[, 5] <- apply(y,2,quantile,0.975)
nl <- as.integer(nItr * 0.025)
nu <- as.integer(nItr * 0.975)
z[, 4] <- apply(y,2,up.low.limit,limit=nl)
z[, 5] <- apply(y,2,up.low.limit,limit=nu)
as.data.frame(z)
}
##
## segment plot for upper and lower limit
##
spT.segment.plot<-function(obs, est, up, low, limit=NULL){
#
tmp<-cbind(obs,est,up,low)
tmp<-na.omit(tmp)
if(is.null(limit)==TRUE){
plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions", pch="*",
xlim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)),
ylim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)))
}
#
else{
plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions",
xlim=c(limit[1],limit[2]),ylim=c(limit[1],limit[2]),pch="*")
}
#
segments(tmp[,1],tmp[,2],tmp[,1],tmp[,3])
segments(tmp[,1],tmp[,2],tmp[,1],tmp[,4])
#
}
##
## hit and false alarm function for forecast
##
spT.hit.false<-function(obs,fore,tol){
#
tmp<-cbind(obs,fore)
tmp<-na.omit(tmp)
#
c11<-tmp[tmp[,1]<=tol & tmp[,2]<=tol,]
c11<-length(c11)/2
c12<-tmp[tmp[,1]<=tol & tmp[,2]>tol,]
c12<-length(c12)/2
c21<-tmp[tmp[,1]>tol & tmp[,2]<=tol,]
c21<-length(c21)/2
c22<-tmp[tmp[,1]>tol & tmp[,2]>tol,]
c22<-length(c22)/2
mat<-matrix(c(c11,c21,c12,c22),2,2)
dimnames(mat)[[1]]<-c(paste("[Obs:<=",tol,"]"),paste("[Obs:> ",tol,"]"))
dimnames(mat)[[2]]<-c(paste("[Forecast:<=",tol,"]"),paste("[Forecast:>",tol,"]"))
POD<-round(mat[1,1]/sum(diag(mat)),4)
FAR<-round(mat[1,2]/sum(mat[1,]),4)
HAR<-round(sum(diag(mat))/sum(mat),4)
top<-2*(mat[1,1]*mat[2,2]-mat[1,2]*mat[2,1])
bot<-mat[1,2]^2+mat[2,1]^2+2*mat[1,1]*mat[2,2]+(mat[1,2]+mat[2,1])*sum(diag(mat))
S<-round(top/bot,4)
x<-list(False.Alarm=FAR,Hit.Rate=HAR,Probability.of.Detection=POD,
Heidke.Skill=S,cross.table=mat,tolerance.limit=tol)
x
}
##
## To keep the long lat positions based on tol.dist
##
spT.keep.morethan.dist <- function(coords, tol.dist=100)
{
#
# coords must have two columns named long and lat
#
a <- as.data.frame(coords)
names(a) <- c("long","lat")
n <- nrow(coords)
c1 <- rep(1:n, each=n)
c2 <- rep(1:n, n)
b <- matrix(NA, nrow=n, ncol=n)
w <- as.vector(upper.tri(b))
bigmat <- matrix(0, nrow=n*n, ncol=7)
bigmat[, 1] <- c1
bigmat[, 2] <- c2
bigmat[, 3] <- a$long[c1]
bigmat[, 4] <- a$lat[c1]
bigmat[, 5] <- a$long[c2]
bigmat[, 6] <- a$lat[c2]
ubmat <- bigmat[w, ]
ubmat[,7] <- as.vector(apply(ubmat[,3:6], 1, spT.geo_dist))
v <- ubmat[,7]
w <- ubmat[v<tol.dist, ]
z <- unique(w[,1])
a <- coords[-z, ]
a
}
##
## Predictive Model Choice Criteria
##
PMCC<-function(z=NULL, z.mean=NULL, z.sd=NULL, z.samples=NULL)
{
#
# Predictive Model Choice Criteria
#
if(is.null(z)){
stop("Error: need to provide z values.")
}
#
if(is.null(z.mean) | is.null(z.sd)){
if(!is.null(z.mean)){
stop("Error: need to provide z.sd value.")
}
if(!is.null(z.sd)){
stop("Error: need to provide z.mean value.")
}
if(is.null(z.samples)){
stop("Error: need to provide z.samples value.")
}
}
#
if(!is.null(z.samples)){
if ( !is.matrix(z.samples)) {
stop("Error: z.samples must be a (N x nItr) matrix")
}
if (dim(z.samples)[1] != length(z)) {
stop("Error: observations in z.samples in each iteration must be equal to length of z")
}
if ( dim(z.samples)[2] < 40) {
stop("Error: samples are too small to obtain summary statistics")
}
sum.stat = matrix(NA,length(c(z)),6)
sum.stat[,1:5] = as.matrix(spT.Summary.Stat(z.samples))
sum.stat[,6] = c(z)
sum.stat = sum.stat[!is.na(sum.stat[,6]),]
goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,6])^2),2)
penalty = round(sum(sum.stat[,3]^2),2)
pmcc = round(goodness.of.fit + penalty,2)
out = NULL
out$pmcc = pmcc;
out$goodness.of.fit = goodness.of.fit
out$penalty = penalty
#class(out) <- "PMCC"
out
}
else{
if(is.null(z.mean) | is.null(z.sd)){
stop("Error: need to provide z.mean and/or z.sd values.")
}
if(length(c(z)) != length(c(z.mean))){
stop("Error: z and z.mean should be in same length.")
}
if(length(c(z)) != length(c(z.sd))){
stop("Error: z and z.sd should be in same length.")
}
sum.stat = matrix(NA,length(c(z)),3)
sum.stat[,1] = c(z)
sum.stat[,2] = c(z.mean)
sum.stat[,3] = c(z.sd)
sum.stat = sum.stat[!is.na(sum.stat[,1]),]
goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,2])^2),2)
penalty = round(sum(sum.stat[,3]^2),2)
pmcc = round(goodness.of.fit + penalty,2)
out = NULL
out$pmcc = pmcc;
out$goodness.of.fit = goodness.of.fit
out$penalty = penalty
#class(out) <- "PMCC"
out
}
}
##
## To identify the Integer Number
##
#is.wholenumber<-function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
##
## To check sites that has <tol km of distances
##
spT.check.locations<-function(fit.locations, pred.locations,
method="geodetic:km", tol=5){
#
#
if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
coords.all <- rbind(fit.locations,pred.locations)
tn.fitsites <- length(fit.locations[, 1])
nfit.sites <- 1:tn.fitsites
tn.predsites <- length(coords.all[, 1]) - tn.fitsites
npred.sites <- (tn.fitsites + 1):(length(coords.all[, 1]))
#
if(method=="geodetic:km"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
}
else if(method=="geodetic:mile"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
}
else{
coords.D <- as.matrix(dist(coords.all, method, diag = T, upper = T))
}
coords.D[is.na(coords.D)]<-0
#
diag(coords.D)<-NA
#
fdmis<-coords.D[nfit.sites, npred.sites]
if(is.matrix(fdmis)==TRUE){
fdmis<-cbind(c(t(fdmis)),1:dim(fdmis)[[2]],sort(rep(1:dim(fdmis)[[1]],dim(fdmis)[[2]]))) # add pred sites and fitted sites
fdmis<-fdmis[fdmis[,1] < tol,]
#
if(!is.na(fdmis[1])==TRUE){
cat("#\n# Tolerance Limit:", paste(tol))
cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
fdmis<-matrix(fdmis,(length(fdmis)/3),3)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",fdmis[i,2]," Fitted location:", fdmis[i,3],""))
}
cat("#\n# Romove the locations and run again. \n#\n")
dimnames(fdmis)[[2]]<-c('distance','pred_location','fit_location')
fdmis
}
else{
cat("#\n# Tolerance Limit (unit):", paste(tol))
#cat("\n# Fitted and Predicted location distances are alright \n#\n")
cat("\n# Location distances are alright \n#\n")
}
}
else{
fdmis<-cbind(c(fdmis),1:length(fdmis)) #
fdmis<-fdmis[fdmis[,1] < tol,]
if(!is.na(fdmis[1])==TRUE){
cat("#\n# Tolerance Limit:", paste(tol))
cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
fdmis<-matrix(fdmis)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",1," Fitted location:", fdmis[i,2],""))
}
cat("#\n# Romove the locations and run again. \n#\n")
dimnames(fdmis)[[2]]<-c('distance','fit_location')
fdmis
}
else{
cat("#\n# Tolerance Limit (unit):", paste(tol))
#cat("\n# Fitted and Predicted location distances are alright \n#\n")
cat("\n# Location distances are alright \n#\n")
}
}
}
##
## To check sites inside codes
##
spT.check.sites.inside<-function(coords, method, tol=0.1){
#
#
if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
#
if(method=="geodetic:km"){
fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=TRUE))
}
else if(method=="geodetic:mile"){
fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=FALSE))
}
else{
fdm<- as.matrix(dist(coords, method, diag = TRUE, upper = TRUE))
}
#
diag(fdm)<-NA
#
fdm<-cbind(c(fdm),1:dim(fdm)[[2]],sort(rep(1:dim(fdm)[[1]],dim(fdm)[[2]])))
#
fdm<-fdm[!is.na(fdm[,1]),]
#
#tol <- 0.01
fdmis<-fdm[fdm[,1] < tol,]
#
if(!is.na(fdmis[1])==TRUE){
cat("#\n# There are some locations very close:")
print(paste("( < ",tol," unit) to each other."))
fdmis<-matrix(fdmis,(length(fdmis)/3),3)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance (unit):", round(fdmis[i,1],2)," site:",fdmis[i,2]," site:", fdmis[i,3],""))
}
dimnames(fdmis)[[2]]<-c('dist_km','pred_site','fit_site')
fdmis
cat("#\n# Romove the sites and run again. \n#\n")
stop("Error: Termination.")
}
#
#tol <- 0.5
#fdmis<-fdm[fdm[,1] < tol,]
#
#if(!is.na(fdmis[1])==TRUE){
# cat("#\n# Warnings: There are some locations very close ( < 0.5 unit) to each other.\n#\n")
#}
#
}
##
##
#check.sites.inside.pred<-spT.check.locations
##
## convert seconds into min. hour. and day
##
fnc.time<-function(t)
{
#
if(t < 60){
t <- round(t,2)
tt <- paste(t," - Sec.")
cat(paste("##\n# Elapsed time:",t,"Sec.\n##\n"))
}
#
if(t < (60*60) && t >= 60){
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t1,"Min.",t,"Sec.\n##\n"))
}
#
if(t < (60*60*24) && t >= (60*60)){
t2 <- as.integer(t/(60*60))
t <- t-t2*60*60
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t2," - Hour/s.",t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t2,"Hour/s.",t1,"Min.",t,"Sec.\n##\n"))
}
#
if(t >= (60*60*24)){
t3 <- as.integer(t/(60*60*24))
t <- t-t3*60*60*24
t2 <- as.integer(t/(60*60))
t <- t-t2*60*60
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t3," - Day/s.",t2," - Hour/s.",t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t3,"Day/s.",t2,"Hour/s.",t1,"Mins.",t,"Sec.\n##\n"))
}
#
tt
}
##
## For spatially varying beta
##
sp<-function(x){
class(x)<-"spBeta"
x
}
##
## For temporally varying beta
##
tp<-function(x){
class(x)<-"tpBeta"
x
}
##
## find close locations
##
ObsGridLoc<-function(obsLoc, gridLoc, distance.method="geodetic:km",
plot=FALSE)
{
#
#
if(!distance.method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
coords.all <- rbind(obsLoc, gridLoc)
n1.sites <- 1:dim(obsLoc)[[1]]
n2.sites <- (dim(obsLoc)[[1]] + 1):(dim(coords.all)[[1]])
#
if(distance.method=="geodetic:km"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
}
else if(distance.method=="geodetic:mile"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
}
else{
coords.D <- as.matrix(dist(coords.all, distance.method, diag = T, upper = T))
}
#
coords.D[is.na(coords.D)]<-0
diag(coords.D)<-NA
#
fdmis<-coords.D[n1.sites, n2.sites]
min.fnc<-function(x){
x<-cbind(x,1:length(x))
x<-x[order(x[,1]),]
x<-x[1,]
x
}
fdmis<-t(apply(fdmis,1,min.fnc))
fdmis<-cbind(obsLoc,gridLoc[fdmis[,2],],fdmis[,2],fdmis[,1])
dimnames(fdmis)[[2]]<-c("obsLon","obsLat","gridLon","gridLat","gridNum","dist")
if(plot==TRUE){
plot(fdmis[,1:2],pch="*")
points(fdmis[,3:4],pch=19,col=2)
legend("bottomleft",pch=c(8,19),col=c(1,2),legend=c("Observation locations", "Grid locations"),cex=0.9)
}
fdmis
}
##
## combine grid data and coords for spTimer input, gridData should be in array with dim = lon, lat, day, year
##
gridTodata<-function(gridData, gridLoc=NULL, gridLon=NULL, gridLat=NULL)
{
# gridData should be in array with dim = lon, lat, day, year
# for morethan one array use list argument
# gridLoc should be obtain using expand.grid() argument
options(warn=-1)
if(!is.null(gridLoc)){
lonlat<-gridLoc
dimnames(lonlat)[[2]]<-c("lon","lat")
}
else if((!is.null(gridLon)) & (!is.null(gridLat))){
lonlat<-expand.grid(gridLon, gridLat)
dimnames(lonlat)[[2]]<-c("lon","lat")
}
else{
stop("Error: check grid location input in the function")
}
if(is.list(gridData)){
n<-length(gridData)
if(n>1){
dat<-NULL
for(i in 1:n){
grid<-c(gridData[[i]])
dat<-cbind(dat,grid)
}
for(i in 1:n){ dimnames(dat)[[2]][i]<-paste("grid",i,sep="") }
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),dat)
}
else{
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),grid=c(gridData))
}
}
else{
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData))/dim(lonlat)[[1]]),grid=c(gridData))
}
dat<-dat[order(dat[,1],dat[,2]),]
dat
}
##
## combine observation and grid data for model fitting
##
ObsGridData<-function(obsData, gridData, obsLoc, gridLoc, distance.method="geodetic:km")
{
# obsData should be in data form
# gridData should be in array form with dim = lon, lat, day, year
obsLoc<-as.matrix(obsLoc)
gridLoc<-as.matrix(gridLoc)
dimnames(obsLoc)<-NULL
dimnames(gridLoc)<-NULL
if((is.list(gridData)) & (is.list(gridLoc))){
n1<-length(gridData)
n2<-length(gridLoc)
if(n1 != n2){ stop("\n Error: list of gridData and gridLoc should be same") }
datt<-NULL
for(i in 1:n1){
gLoc<-gridLoc[[i]]
gLoc<-gLoc[order(gLoc[,1],gLoc[,2]),]
comb<-ObsGridLoc(obsLoc=obsLoc, gridLoc=gLoc, distance.method=distance.method)
grid<-gridTodata(gridData=gridData[[i]], gridLoc=gLoc)
grid<-cbind(gridNum=sort(rep(1:dim(gLoc)[[1]],dim(grid)[[1]]/dim(gLoc)[[1]])),grid)
dat<-NULL
for(j in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==j,]) }
dimnames(dat)[[2]]<-c(paste("gridNum",i,sep=""),paste("grid.lon",i,sep=""),paste("grid.lat",i,sep=""),paste("grid",i,sep=""))
if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n obsData and gridData time points mismatch.\n")}
datt<-cbind(datt,dat)
}
dat<-cbind(obsData,datt); rm(datt);
dat
}
else{
comb<-ObsGridLoc(obsLoc=obsLoc, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),], distance.method=distance.method)
grid<-gridTodata(gridData=gridData, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),])
grid<-cbind(gridNum=sort(rep(1:dim(gridLoc)[[1]],dim(grid)[[1]]/dim(gridLoc)[[1]])),grid)
dat<-NULL
for(i in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==i,]) }
dimnames(dat)[[2]][2:3]<-c("grid.lon","grid.lat")
if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n obsData and gridData time points mismatch.\n")}
dat<-cbind(obsData,dat)
dat
}
}
##
## Truncated Gaussian code
##
truncated.fnc<-function(Y, at=0, lambda=NULL, both=FALSE){
#
#
# at is left-tailed
#
if(is.null(lambda)){
stop("Error: define truncation parameter lambda properly using list ")
}
if(is.null(at)){
stop("Error: define truncation point properly using list ")
}
if(at < 0){
stop("Error: currently truncation point only can take value >= zero ")
}
zm <- cbind(Y,Y-at,1)
zm[zm[, 2] <= 0, 3] <- 0
zm[, 2] <- zm[, 2]^(1/lambda)
zm[zm[, 3] == 0, 2] <- -(rgamma(nrow(zm[zm[,3]==0,]), shape=1, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#mat <- matrix(NA,nrow(zm[zm[,3]==0,]),10)
#mat[,] <- (rgamma(nrow(zm[zm[,3]==0,])*10, shape=1, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#mat[,] <- -(rexp(nrow(zm[zm[,3]==0,])*10, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#zm[zm[, 3] == 0, 2] <- -apply(mat,1,max)
if (both == TRUE) {
zm
}
else {
c(zm[,2])
}
#
}
##
## for truncated model
##
reverse.truncated.fnc<-function(Y, at=0, lambda=NULL){
#
# at is left-tailed
#
if(is.null(lambda)){
stop("Error: define truncation parameter lambda properly using list ")
}
zm <- Y
zm[zm <= 0]<- 0
zm <- zm^(lambda)
zm <- zm + at
zm
#
}
##
## probability below threshold (for truncated)
##
prob.below.threshold <- function(out, at){
# out is the N x nItr MCMC samples
fnc<-function(x, at){
length(x[x <= at])/length(x)
}
apply(out,1,fnc,at=at)
}
##
##
##
Try the spTDyn package in your browser
Any scripts or data that you put into this service are public.
spTDyn documentation built on Nov. 22, 2022, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions prob.below.threshold reverse.truncated.fnc truncated.fnc ObsGridData gridTodata ObsGridLoc tp sp fnc.time spT.check.sites.inside spT.check.locations PMCC spT.keep.morethan.dist spT.hit.false spT.segment.plot spT.Summary.Stat tp.dimname.fnc sp.dimname.fnc plot.spTD summary.spTD as.mcmc.spTD residuals.spTD coef.spTD confint.spTD fitted.spTD print.spTD
Documented in as.mcmc.spTD coef.spTD confint.spTD fitted.spTD fnc.time gridTodata ObsGridData ObsGridLoc plot.spTD PMCC print.spTD prob.below.threshold residuals.spTD reverse.truncated.fnc sp sp.dimname.fnc spT.check.locations spT.check.sites.inside spT.hit.false spT.keep.morethan.dist spT.segment.plot spT.Summary.Stat summary.spTD tp tp.dimname.fnc truncated.fnc
##
## Print function
##
print.spTD<-function(x, ...) {
cat("-----------------------------------------------------"); cat('\n');
cat("Model: "); cat(x$model); cat('\n');
cat("Call: "); print(x$call); #cat('\n')
cat("Iterations: "); cat(x$iterations); cat("\n")
cat("nBurn: "); cat(x$nBurn); cat("\n")
cat("Acceptance rate for phi (%): "); cat(x$accept); cat("\n")
cat("-----------------------------------------------------"); cat('\n');
#cat("PMCC: "); cat("\n");
print(x$PMCC);
cat("-----------------------------------------------------"); cat('\n');
#cat("Parameters:\n")
#print(x$parameter); #cat("\n");
#cat("-----------------------------------------------------"); cat('\n');
cat("Computation time: "); cat(x$computation.time); cat("\n")
}
##
## fitted function
##
fitted.spTD<-function(object, ...){
x<-data.frame(object$fitted)
x
}
##
## use of confint()
##
confint.spTD<-function(object, parm, level=0.95, ...){
#
x<-as.mcmc(object)
up<-level+(1-level)/2
low<-(1-level)/2
FUN <- function(x){quantile(x,probs=c(low,up))}
out<-apply(x,2,FUN=FUN)
out<-t(out)
if(missing(parm)){
out
}
else{
if(length(parm)>1){
out<-as.matrix(out[dimnames(out)[[1]] %in% parm,])
out
}
else{
out<-t(as.matrix(out[dimnames(out)[[1]] %in% parm,]))
dimnames(out)[[1]]<-parm
out
}
}
}
##
## use of coefficients
##
coef.spTD<-function(object, digits=4, ...){
round(t(object$parameter)[1,],digits=digits)
}
##
## use of residuals
##
residuals.spTD<-function(object, ...){
#if(object$combined.fit.pred==TRUE){
# stop("\n# Error: not useful for output with combined fit and predict")
#}
#else{
if(object$scale.transform=="NONE"){
tmp<-object$Y-object$fitted[,1]
tmp
}
else if(object$scale.transform=="SQRT"){
tmp<-sqrt(object$Y)-object$fitted[,1]
tmp
}
else if(object$scale.transform=="LOG"){
tmp<-log(object$Y)-object$fitted[,1]
tmp
}
else{
}
#}
}
##
## To use in coda package
##
as.mcmc.spTD<-function(x, ...){
if (x$combined.fit.pred == TRUE) {
stop("\n# Error: not available for combined fit-prediction option. Please read the MCMC samples from the *.txt file manually.")
}
model <- x$model
if (is.null(model) == TRUE) {
stop("\n# Error: need to define the model")
}
else if (model == "GP" | model == "truncatedGP") {
r <- x$r
p <- x$p
#
if(x$cov.fnc=="matern"){
if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi", "nu")
}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi", "nu")
}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi", "nu")
}
else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip), t(x$nup))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi", "nu")
}
else {
stop("Error")
}
}
else {
if((is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "phi")
}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "phi")
}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2deltap", "sig2op", "phi")
}
else if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){
para <- rbind((x$betap), t(x$sig2ep), t(x$sig2etap), t(x$sig2betap), t(x$sig2deltap), t(x$sig2op), t(x$phip))
dimnames(para)[[1]] <- c(dimnames(x$X)[[2]], "sig2eps", "sig2eta", "sig2beta", "sig2deltap", "sig2op", "phi")
}
else {
stop("Error")
}
}
#
para<-t(para)
para<-mcmc(para)
para
}
else {
}
}
##
## use of summary
##
summary.spTD<-function(object, digits=4, package="spTDyn", coefficient=NULL, ...){
#
if(package=="coda"){
if(object$combined.fit.pred==TRUE){
stop("\n# Error: coda package is not useful for output with combined fit and predict \n")
}
else{
if(is.null(coefficient)){
cat("\n#### MCMC summary statistics using coda package ####\n")
if(!is.null(object$sp.covariate.names)) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
tmp<-as.mcmc(object)
summary(tmp, ...)
}
else if(coefficient=="spatial"){
if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betasp))
if(object$model=="GPP"){n<-object$knots}
else{n<-object$n}
tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q)
summary(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the temporally varying dynamic coefficients ####\n")
tmp<-as.mcmc(t(object$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T)
summary(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
tmp<-as.mcmc(t(object$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="")
summary(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
}
}
}
else{
if(package=="spTDyn"){
coefficient <- coefficient
}
#
else if("Xsp"%in%names(object) | "Xtp"%in%names(object)){
coefficient <- coefficient
}
else{
coefficient <- NULL
}
#
if(is.null(coefficient)){
if((!is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
else if((!is.null(object$sp.covariate.names)) & (is.null(object$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
else if((is.null(object$sp.covariate.names)) & (!is.null(object$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
else { }
#ans<-NULL
#ans$Model=object$model;ans$Call=object$call;ans$Iterations=object$iterations;ans$nBurn=object$nBurn;ans$PMCC=object$pmcc;ans$Parameters=round(object$parameter,digits=digits)
#ans
print(object)
cat("-----------------------------------------------------"); cat('\n');
cat("Parameters:\n")
print(round(object$parameter,digits=digits)); #cat("\n");
cat("-----------------------------------------------------"); cat('\n');
}
else if(coefficient=="spatial"){
if(is.null(object$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betasp))
if(object$model=="GPP"){n<-object$knots}
else{n<-object$n}
tmp<-sp.dimname.fnc(x=tmp,names=object$sp.covariate.names,n=n,q=object$q)
summary(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(object$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
cat("\n#### MCMC summary statistics for the spatially varying coefficients ####\n")
tmp<-as.mcmc(t(object$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=object$tp.covariate.names,u=object$u,T=object$T)
summary(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(object$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
cat("\n#### MCMC summary statistics for the rho coefficients ####\n")
tmp<-as.mcmc(t(object$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:object$u,sep="")
summary(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes character 'spatial' and 'temporal'.")
}
}
}
##
## use of plot
##
plot.spTD<-function(x, residuals=FALSE, coefficient=NULL, ...){
if(as.logical(residuals)==FALSE){
if(x$combined.fit.pred==TRUE){
if(!is.null(x$sp.covariate.names)) {cat("## \n# Spatially varying parameters are not included.\n##\n ")}
cat("\n## Only residual plots are available for output with combined fit and predict option \n\n")
plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
par(ask=TRUE)
qqnorm(residuals(x));qqline(residuals(x),lty=2)
}
else{
if(is.null(coefficient)){
if((!is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))){cat("\n## \n# Spatially and temporally varying parameters are not included.\n##\n ")}
else if((!is.null(x$sp.covariate.names)) & (is.null(x$tp.covariate.names))) {cat("\n## \n# Spatially varying parameters are not included.\n##\n ")}
else if((is.null(x$sp.covariate.names)) & (!is.null(x$tp.covariate.names))) {cat("\n## \n# Temporally varying dynamic parameters are not included.\n##\n ")}
else { }
if(x$model=="GP"){
tmp<-as.mcmc(x)
plot(tmp, ...)
}
else{
x$model="GP"
tmp<-as.mcmc(x)
plot(tmp, ...)
}
}
else if(coefficient=="spatial"){
if(is.null(x$sp.covariate.names)) {stop("\n## \n# Spatially varying parameters are not used in the model.\n##\n ")}
tmp<-as.mcmc(t(x$betasp))
if(x$model=="GPP"){n<-x$knots}
else{n<-x$n}
tmp<-sp.dimname.fnc(x=tmp,names=x$sp.covariate.names,n=n,q=x$q)
plot(tmp, ...)
}
else if(coefficient=="temporal"){
if(is.null(x$tp.covariate.names)) {stop("\n## \n# Temporally varying dynamic parameters are not used in the model.\n##\n ")}
tmp<-as.mcmc(t(x$betatp))
tmp<-tp.dimname.fnc(x=tmp,names=x$tp.covariate.names,u=x$u,T=x$T)
plot(tmp, ...)
}
else if(coefficient=="rho"){
if(is.null(x$rhotp)) {stop("\n## \n# rho parameter has not been sampled in the temporally varying dynamic model.\n##\n ")}
tmp<-as.mcmc(t(x$rhotp))
dimnames(tmp)[[2]]<-paste("rho",1:x$u,sep="")
plot(tmp, ...)
}
else{
stop("Error: the argument coefficient only takes charecter 'spatial' and 'temporal'.")
}
}
}
else{
if(!is.null(x$sp.covariate.names)) {cat("## \n# Models fitted with spatially varying parameters.\n## \n ")}
plot(x$fitted[,1],residuals(x),ylab="Residuals",xlab="Fitted values");abline(h=0,lty=2);title("Residuals vs Fitted")
par(ask=TRUE)
qqnorm(residuals(x));qqline(residuals(x),lty=2)
} # end of 2nd if
}
##
## to include dimnames: function used in function summary statistics for spatially varying coef
##
sp.dimname.fnc<-function(x,names,n,q){
dimnames(x)[[2]][1:(n*q)]<-1:(n*q)
for(i in 1:q){
dimnames(x)[[2]][(1+(i-1)*n):(n*i)]<-rep(paste(names[i],"site",1:n,sep=""))
}
x
}
##
## to include dimnames: function used in function summary statistics for temporally varying coef
##
tp.dimname.fnc<-function(x,names,u,T){
dimnames(x)[[2]][1:(u*T)]<-1:(u*T)
for(i in 1:u){
dimnames(x)[[2]][(1+(i-1)*T):(T*i)]<-rep(paste(names[i],"time",1:T,sep=""))
}
x
}
##
## Summary Statistics
##
spT.Summary.Stat <- function(y)
{
## define Upper and lower limit function
up.low.limit<-function(y,limit)
{
#
y<-sort(y)
y<-y[limit]
y
#
}
#
if(is.vector(y)==TRUE){
y <- matrix(y[!is.na(y)])
}
else{
y <- t(y)
}
N <- length(y[1, ])
nItr <- length(y[, 1])
z <- matrix(nrow = N, ncol = 5)
dimnames(z) <- list(dimnames(y)[[2]], c("Mean","Median","SD","Low2.5p","Up97.5p"))
#
if (nItr < 40) {
stop("\n##\n# Error: number of samples must be >= 40\n##")
}
z[, 1] <- apply(y,2,mean)
z[, 2] <- apply(y,2,median)
z[, 3] <- apply(y,2,sd)
#z[, 4] <- apply(y,2,quantile,0.025)
#z[, 5] <- apply(y,2,quantile,0.975)
nl <- as.integer(nItr * 0.025)
nu <- as.integer(nItr * 0.975)
z[, 4] <- apply(y,2,up.low.limit,limit=nl)
z[, 5] <- apply(y,2,up.low.limit,limit=nu)
as.data.frame(z)
}
##
## segment plot for upper and lower limit
##
spT.segment.plot<-function(obs, est, up, low, limit=NULL){
#
tmp<-cbind(obs,est,up,low)
tmp<-na.omit(tmp)
if(is.null(limit)==TRUE){
plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions", pch="*",
xlim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)),
ylim=c(min(c(tmp),na.rm=TRUE),max(c(tmp),na.rm=TRUE)))
}
#
else{
plot(tmp[,1],tmp[,2],xlab="Observations",ylab="Predictions",
xlim=c(limit[1],limit[2]),ylim=c(limit[1],limit[2]),pch="*")
}
#
segments(tmp[,1],tmp[,2],tmp[,1],tmp[,3])
segments(tmp[,1],tmp[,2],tmp[,1],tmp[,4])
#
}
##
## hit and false alarm function for forecast
##
spT.hit.false<-function(obs,fore,tol){
#
tmp<-cbind(obs,fore)
tmp<-na.omit(tmp)
#
c11<-tmp[tmp[,1]<=tol & tmp[,2]<=tol,]
c11<-length(c11)/2
c12<-tmp[tmp[,1]<=tol & tmp[,2]>tol,]
c12<-length(c12)/2
c21<-tmp[tmp[,1]>tol & tmp[,2]<=tol,]
c21<-length(c21)/2
c22<-tmp[tmp[,1]>tol & tmp[,2]>tol,]
c22<-length(c22)/2
mat<-matrix(c(c11,c21,c12,c22),2,2)
dimnames(mat)[[1]]<-c(paste("[Obs:<=",tol,"]"),paste("[Obs:> ",tol,"]"))
dimnames(mat)[[2]]<-c(paste("[Forecast:<=",tol,"]"),paste("[Forecast:>",tol,"]"))
POD<-round(mat[1,1]/sum(diag(mat)),4)
FAR<-round(mat[1,2]/sum(mat[1,]),4)
HAR<-round(sum(diag(mat))/sum(mat),4)
top<-2*(mat[1,1]*mat[2,2]-mat[1,2]*mat[2,1])
bot<-mat[1,2]^2+mat[2,1]^2+2*mat[1,1]*mat[2,2]+(mat[1,2]+mat[2,1])*sum(diag(mat))
S<-round(top/bot,4)
x<-list(False.Alarm=FAR,Hit.Rate=HAR,Probability.of.Detection=POD,
Heidke.Skill=S,cross.table=mat,tolerance.limit=tol)
x
}
##
## To keep the long lat positions based on tol.dist
##
spT.keep.morethan.dist <- function(coords, tol.dist=100)
{
#
# coords must have two columns named long and lat
#
a <- as.data.frame(coords)
names(a) <- c("long","lat")
n <- nrow(coords)
c1 <- rep(1:n, each=n)
c2 <- rep(1:n, n)
b <- matrix(NA, nrow=n, ncol=n)
w <- as.vector(upper.tri(b))
bigmat <- matrix(0, nrow=n*n, ncol=7)
bigmat[, 1] <- c1
bigmat[, 2] <- c2
bigmat[, 3] <- a$long[c1]
bigmat[, 4] <- a$lat[c1]
bigmat[, 5] <- a$long[c2]
bigmat[, 6] <- a$lat[c2]
ubmat <- bigmat[w, ]
ubmat[,7] <- as.vector(apply(ubmat[,3:6], 1, spT.geo_dist))
v <- ubmat[,7]
w <- ubmat[v<tol.dist, ]
z <- unique(w[,1])
a <- coords[-z, ]
a
}
##
## Predictive Model Choice Criteria
##
PMCC<-function(z=NULL, z.mean=NULL, z.sd=NULL, z.samples=NULL)
{
#
# Predictive Model Choice Criteria
#
if(is.null(z)){
stop("Error: need to provide z values.")
}
#
if(is.null(z.mean) | is.null(z.sd)){
if(!is.null(z.mean)){
stop("Error: need to provide z.sd value.")
}
if(!is.null(z.sd)){
stop("Error: need to provide z.mean value.")
}
if(is.null(z.samples)){
stop("Error: need to provide z.samples value.")
}
}
#
if(!is.null(z.samples)){
if ( !is.matrix(z.samples)) {
stop("Error: z.samples must be a (N x nItr) matrix")
}
if (dim(z.samples)[1] != length(z)) {
stop("Error: observations in z.samples in each iteration must be equal to length of z")
}
if ( dim(z.samples)[2] < 40) {
stop("Error: samples are too small to obtain summary statistics")
}
sum.stat = matrix(NA,length(c(z)),6)
sum.stat[,1:5] = as.matrix(spT.Summary.Stat(z.samples))
sum.stat[,6] = c(z)
sum.stat = sum.stat[!is.na(sum.stat[,6]),]
goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,6])^2),2)
penalty = round(sum(sum.stat[,3]^2),2)
pmcc = round(goodness.of.fit + penalty,2)
out = NULL
out$pmcc = pmcc;
out$goodness.of.fit = goodness.of.fit
out$penalty = penalty
#class(out) <- "PMCC"
out
}
else{
if(is.null(z.mean) | is.null(z.sd)){
stop("Error: need to provide z.mean and/or z.sd values.")
}
if(length(c(z)) != length(c(z.mean))){
stop("Error: z and z.mean should be in same length.")
}
if(length(c(z)) != length(c(z.sd))){
stop("Error: z and z.sd should be in same length.")
}
sum.stat = matrix(NA,length(c(z)),3)
sum.stat[,1] = c(z)
sum.stat[,2] = c(z.mean)
sum.stat[,3] = c(z.sd)
sum.stat = sum.stat[!is.na(sum.stat[,1]),]
goodness.of.fit = round(sum((sum.stat[,1]-sum.stat[,2])^2),2)
penalty = round(sum(sum.stat[,3]^2),2)
pmcc = round(goodness.of.fit + penalty,2)
out = NULL
out$pmcc = pmcc;
out$goodness.of.fit = goodness.of.fit
out$penalty = penalty
#class(out) <- "PMCC"
out
}
}
##
## To identify the Integer Number
##
#is.wholenumber<-function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
##
## To check sites that has <tol km of distances
##
spT.check.locations<-function(fit.locations, pred.locations,
method="geodetic:km", tol=5){
#
#
if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
coords.all <- rbind(fit.locations,pred.locations)
tn.fitsites <- length(fit.locations[, 1])
nfit.sites <- 1:tn.fitsites
tn.predsites <- length(coords.all[, 1]) - tn.fitsites
npred.sites <- (tn.fitsites + 1):(length(coords.all[, 1]))
#
if(method=="geodetic:km"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
}
else if(method=="geodetic:mile"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
}
else{
coords.D <- as.matrix(dist(coords.all, method, diag = T, upper = T))
}
coords.D[is.na(coords.D)]<-0
#
diag(coords.D)<-NA
#
fdmis<-coords.D[nfit.sites, npred.sites]
if(is.matrix(fdmis)==TRUE){
fdmis<-cbind(c(t(fdmis)),1:dim(fdmis)[[2]],sort(rep(1:dim(fdmis)[[1]],dim(fdmis)[[2]]))) # add pred sites and fitted sites
fdmis<-fdmis[fdmis[,1] < tol,]
#
if(!is.na(fdmis[1])==TRUE){
cat("#\n# Tolerance Limit:", paste(tol))
cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
fdmis<-matrix(fdmis,(length(fdmis)/3),3)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",fdmis[i,2]," Fitted location:", fdmis[i,3],""))
}
cat("#\n# Romove the locations and run again. \n#\n")
dimnames(fdmis)[[2]]<-c('distance','pred_location','fit_location')
fdmis
}
else{
cat("#\n# Tolerance Limit (unit):", paste(tol))
#cat("\n# Fitted and Predicted location distances are alright \n#\n")
cat("\n# Location distances are alright \n#\n")
}
}
else{
fdmis<-cbind(c(fdmis),1:length(fdmis)) #
fdmis<-fdmis[fdmis[,1] < tol,]
if(!is.na(fdmis[1])==TRUE){
cat("#\n# Tolerance Limit:", paste(tol))
cat("\n# There are some Prediction locations very close to the Fitted locations.\n#\n")
fdmis<-matrix(fdmis)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance:", round(fdmis[i,1],2)," Predicted location:",1," Fitted location:", fdmis[i,2],""))
}
cat("#\n# Romove the locations and run again. \n#\n")
dimnames(fdmis)[[2]]<-c('distance','fit_location')
fdmis
}
else{
cat("#\n# Tolerance Limit (unit):", paste(tol))
#cat("\n# Fitted and Predicted location distances are alright \n#\n")
cat("\n# Location distances are alright \n#\n")
}
}
}
##
## To check sites inside codes
##
spT.check.sites.inside<-function(coords, method, tol=0.1){
#
#
if(!method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
#
if(method=="geodetic:km"){
fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=TRUE))
}
else if(method=="geodetic:mile"){
fdm<- as.matrix(spT.geodist(Lon=coords[,1],Lat=coords[,2], KM=FALSE))
}
else{
fdm<- as.matrix(dist(coords, method, diag = TRUE, upper = TRUE))
}
#
diag(fdm)<-NA
#
fdm<-cbind(c(fdm),1:dim(fdm)[[2]],sort(rep(1:dim(fdm)[[1]],dim(fdm)[[2]])))
#
fdm<-fdm[!is.na(fdm[,1]),]
#
#tol <- 0.01
fdmis<-fdm[fdm[,1] < tol,]
#
if(!is.na(fdmis[1])==TRUE){
cat("#\n# There are some locations very close:")
print(paste("( < ",tol," unit) to each other."))
fdmis<-matrix(fdmis,(length(fdmis)/3),3)
for(i in 1:dim(fdmis)[[1]]){
print(paste("Distance (unit):", round(fdmis[i,1],2)," site:",fdmis[i,2]," site:", fdmis[i,3],""))
}
dimnames(fdmis)[[2]]<-c('dist_km','pred_site','fit_site')
fdmis
cat("#\n# Romove the sites and run again. \n#\n")
stop("Error: Termination.")
}
#
#tol <- 0.5
#fdmis<-fdm[fdm[,1] < tol,]
#
#if(!is.na(fdmis[1])==TRUE){
# cat("#\n# Warnings: There are some locations very close ( < 0.5 unit) to each other.\n#\n")
#}
#
}
##
##
#check.sites.inside.pred<-spT.check.locations
##
## convert seconds into min. hour. and day
##
fnc.time<-function(t)
{
#
if(t < 60){
t <- round(t,2)
tt <- paste(t," - Sec.")
cat(paste("##\n# Elapsed time:",t,"Sec.\n##\n"))
}
#
if(t < (60*60) && t >= 60){
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t1,"Min.",t,"Sec.\n##\n"))
}
#
if(t < (60*60*24) && t >= (60*60)){
t2 <- as.integer(t/(60*60))
t <- t-t2*60*60
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t2," - Hour/s.",t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t2,"Hour/s.",t1,"Min.",t,"Sec.\n##\n"))
}
#
if(t >= (60*60*24)){
t3 <- as.integer(t/(60*60*24))
t <- t-t3*60*60*24
t2 <- as.integer(t/(60*60))
t <- t-t2*60*60
t1 <- as.integer(t/60)
t <- round(t-t1*60,2)
tt <- paste(t3," - Day/s.",t2," - Hour/s.",t1," - Mins.",t," - Sec.")
cat(paste("##\n# Elapsed time:",t3,"Day/s.",t2,"Hour/s.",t1,"Mins.",t,"Sec.\n##\n"))
}
#
tt
}
##
## For spatially varying beta
##
sp<-function(x){
class(x)<-"spBeta"
x
}
##
## For temporally varying beta
##
tp<-function(x){
class(x)<-"tpBeta"
x
}
##
## find close locations
##
ObsGridLoc<-function(obsLoc, gridLoc, distance.method="geodetic:km",
plot=FALSE)
{
#
#
if(!distance.method %in% c("geodetic:km", "geodetic:mile", "euclidean",
"maximum", "manhattan", "canberra")){
stop("\n# Error: correctly define distance.method \n")
}
#
coords.all <- rbind(obsLoc, gridLoc)
n1.sites <- 1:dim(obsLoc)[[1]]
n2.sites <- (dim(obsLoc)[[1]] + 1):(dim(coords.all)[[1]])
#
if(distance.method=="geodetic:km"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=TRUE))
}
else if(distance.method=="geodetic:mile"){
coords.D <- as.matrix(spT.geodist(Lon=coords.all[,1],Lat=coords.all[,2], KM=FALSE))
}
else{
coords.D <- as.matrix(dist(coords.all, distance.method, diag = T, upper = T))
}
#
coords.D[is.na(coords.D)]<-0
diag(coords.D)<-NA
#
fdmis<-coords.D[n1.sites, n2.sites]
min.fnc<-function(x){
x<-cbind(x,1:length(x))
x<-x[order(x[,1]),]
x<-x[1,]
x
}
fdmis<-t(apply(fdmis,1,min.fnc))
fdmis<-cbind(obsLoc,gridLoc[fdmis[,2],],fdmis[,2],fdmis[,1])
dimnames(fdmis)[[2]]<-c("obsLon","obsLat","gridLon","gridLat","gridNum","dist")
if(plot==TRUE){
plot(fdmis[,1:2],pch="*")
points(fdmis[,3:4],pch=19,col=2)
legend("bottomleft",pch=c(8,19),col=c(1,2),legend=c("Observation locations", "Grid locations"),cex=0.9)
}
fdmis
}
##
## combine grid data and coords for spTimer input, gridData should be in array with dim = lon, lat, day, year
##
gridTodata<-function(gridData, gridLoc=NULL, gridLon=NULL, gridLat=NULL)
{
# gridData should be in array with dim = lon, lat, day, year
# for morethan one array use list argument
# gridLoc should be obtain using expand.grid() argument
options(warn=-1)
if(!is.null(gridLoc)){
lonlat<-gridLoc
dimnames(lonlat)[[2]]<-c("lon","lat")
}
else if((!is.null(gridLon)) & (!is.null(gridLat))){
lonlat<-expand.grid(gridLon, gridLat)
dimnames(lonlat)[[2]]<-c("lon","lat")
}
else{
stop("Error: check grid location input in the function")
}
if(is.list(gridData)){
n<-length(gridData)
if(n>1){
dat<-NULL
for(i in 1:n){
grid<-c(gridData[[i]])
dat<-cbind(dat,grid)
}
for(i in 1:n){ dimnames(dat)[[2]][i]<-paste("grid",i,sep="") }
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),dat)
}
else{
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData[[1]]))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData[[1]]))/dim(lonlat)[[1]]),grid=c(gridData))
}
}
else{
dat<-cbind(lon=rep(lonlat[,1],length(c(gridData))/dim(lonlat)[[1]]),
lat=rep(lonlat[,2],length(c(gridData))/dim(lonlat)[[1]]),grid=c(gridData))
}
dat<-dat[order(dat[,1],dat[,2]),]
dat
}
##
## combine observation and grid data for model fitting
##
ObsGridData<-function(obsData, gridData, obsLoc, gridLoc, distance.method="geodetic:km")
{
# obsData should be in data form
# gridData should be in array form with dim = lon, lat, day, year
obsLoc<-as.matrix(obsLoc)
gridLoc<-as.matrix(gridLoc)
dimnames(obsLoc)<-NULL
dimnames(gridLoc)<-NULL
if((is.list(gridData)) & (is.list(gridLoc))){
n1<-length(gridData)
n2<-length(gridLoc)
if(n1 != n2){ stop("\n Error: list of gridData and gridLoc should be same") }
datt<-NULL
for(i in 1:n1){
gLoc<-gridLoc[[i]]
gLoc<-gLoc[order(gLoc[,1],gLoc[,2]),]
comb<-ObsGridLoc(obsLoc=obsLoc, gridLoc=gLoc, distance.method=distance.method)
grid<-gridTodata(gridData=gridData[[i]], gridLoc=gLoc)
grid<-cbind(gridNum=sort(rep(1:dim(gLoc)[[1]],dim(grid)[[1]]/dim(gLoc)[[1]])),grid)
dat<-NULL
for(j in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==j,]) }
dimnames(dat)[[2]]<-c(paste("gridNum",i,sep=""),paste("grid.lon",i,sep=""),paste("grid.lat",i,sep=""),paste("grid",i,sep=""))
if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n obsData and gridData time points mismatch.\n")}
datt<-cbind(datt,dat)
}
dat<-cbind(obsData,datt); rm(datt);
dat
}
else{
comb<-ObsGridLoc(obsLoc=obsLoc, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),], distance.method=distance.method)
grid<-gridTodata(gridData=gridData, gridLoc=gridLoc[order(gridLoc[,1],gridLoc[,2]),])
grid<-cbind(gridNum=sort(rep(1:dim(gridLoc)[[1]],dim(grid)[[1]]/dim(gridLoc)[[1]])),grid)
dat<-NULL
for(i in comb[,5]){ dat<-rbind(dat,grid[grid[,1]==i,]) }
dimnames(dat)[[2]][2:3]<-c("grid.lon","grid.lat")
if(dim(obsData)[[1]] != dim(dat)[[1]]){ stop("Error: check the day and year for gridData\n obsData and gridData time points mismatch.\n")}
dat<-cbind(obsData,dat)
dat
}
}
##
## Truncated Gaussian code
##
truncated.fnc<-function(Y, at=0, lambda=NULL, both=FALSE){
#
#
# at is left-tailed
#
if(is.null(lambda)){
stop("Error: define truncation parameter lambda properly using list ")
}
if(is.null(at)){
stop("Error: define truncation point properly using list ")
}
if(at < 0){
stop("Error: currently truncation point only can take value >= zero ")
}
zm <- cbind(Y,Y-at,1)
zm[zm[, 2] <= 0, 3] <- 0
zm[, 2] <- zm[, 2]^(1/lambda)
zm[zm[, 3] == 0, 2] <- -(rgamma(nrow(zm[zm[,3]==0,]), shape=1, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#mat <- matrix(NA,nrow(zm[zm[,3]==0,]),10)
#mat[,] <- (rgamma(nrow(zm[zm[,3]==0,])*10, shape=1, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#mat[,] <- -(rexp(nrow(zm[zm[,3]==0,])*10, rate=1/(range(zm[zm[,3]==1,2],na.rm=TRUE)[[2]]/4.1)))
#zm[zm[, 3] == 0, 2] <- -apply(mat,1,max)
if (both == TRUE) {
zm
}
else {
c(zm[,2])
}
#
}
##
## for truncated model
##
reverse.truncated.fnc<-function(Y, at=0, lambda=NULL){
#
# at is left-tailed
#
if(is.null(lambda)){
stop("Error: define truncation parameter lambda properly using list ")
}
zm <- Y
zm[zm <= 0]<- 0
zm <- zm^(lambda)
zm <- zm + at
zm
#
}
##
## probability below threshold (for truncated)
##
prob.below.threshold <- function(out, at){
# out is the N x nItr MCMC samples
fnc<-function(x, at){
length(x[x <= at])/length(x)
}
apply(out,1,fnc,at=at)
}
##
##
##
Try the spTDyn package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.