Nothing
###
### R routines for the R package dlnm (c)
#
crosspred <-
function(basis, model=NULL, coef=NULL, vcov=NULL, model.link=NULL, at=NULL,
from=NULL, to=NULL, by=NULL, lag, bylag=1, cen=NULL, ci.level=0.95,
cumul=FALSE) {
#
################################################################################
# DETERMINE THE TYPE OF MODEL AND CHECKS
#
# TYPE OF PREDICTION: CROSSBASIS, ONEBASIS, OR PENALIZED GAM
type <- if(any(class(basis)%in%"crossbasis")) "cb" else
if(any(class(basis)%in%"onebasis")) "one" else "gam"
#
# CHECKS ON TYPE, AND SET name, basis AND RESET type
errormes <- "arguments 'basis' and 'model' not consistent. See help(crosspred)"
if(type=="gam") {
if(!is.character(basis) || length(basis)>1L) stop(errormes)
if(is.null(model) || !any(class(model)%in%"gam")) stop(errormes)
name <- basis
sterms <- sapply(model$smooth, function(x) x$term[1])
if(name%in%sterms) basis <- model$smooth[[which(sterms==name)[1]]] else
stop(errormes)
if(length(which(sterms==name)) > 1)
warning(paste(name,"included in multiple smoothers, only the first one taken"))
if(!"cb.smooth"%in%class(basis) && basis$dim>1L)
stop("predictions not provided for multi-dimensional smoothers other than 'cb'")
} else name <- deparse(substitute(basis))
#
# EXTRACT ORIGINAL lag (DEPENDENT ON TYPE)
origlag <- switch(type,
cb = attr(basis,"lag"),
one = c(0,0),
gam = if(is.null(basis$lag)) c(0,0) else basis$lag
)
lag <- if(missing(lag)) origlag else mklag(lag)
#
# CHECKS ON lag AND bylag
if(!all(lag==origlag) && cumul)
stop("cumulative prediction not allowed for lag sub-period")
lagfun <- switch(type, cb=attr(basis,"arglag")$fun, one=NULL,
gam=if(basis$dim==1L) NULL else basis$margin[[2]]$fun)
if(bylag!=1L && !is.null(lagfun) && lagfun=="integer")
stop("prediction for non-integer lags not allowed for type 'integer'")
#
# OTHER COHERENCE CHECKS
if(is.null(model) && (is.null(coef) || is.null(vcov)))
stop("At least 'model' or 'coef'-'vcov' must be provided")
if(!is.numeric(ci.level) || ci.level>=1 || ci.level<=0)
stop("'ci.level' must be numeric and between 0 and 1")
#
################################################################################
# SET COEF, VCOV CLASS AND LINK FOR EVERY TYPE OF MODELS
#
# WRITE CONDITIONS (DEPENDENT ON TYPE AND IF MATRIX/VECTOR)
cond <- if(type=="gam") with(basis,first.para:last.para) else
if(ncol(basis)==1L) name else
if(type=="one") paste(name,"[[:print:]]*b[0-9]{1,2}",sep="") else
paste(name,"[[:print:]]*v[0-9]{1,2}\\.l[0-9]{1,2}",sep="")
#
# IF MODEL PROVIDED, EXTRACT FROM HERE, OTHERWISE DIRECTLY FROM COEF AND VCOV
if(!is.null(model)) {
model.class <- class(model)
coef <- getcoef(model, model.class)
vcov <- getvcov(model, model.class)
indcoef <- if(type=="gam") cond else grep(cond,names(coef))
indvcov <- if(type=="gam") cond else grep(cond,rownames(vcov))
coef <- coef[indcoef]
vcov <- vcov[indvcov,indvcov,drop=FALSE]
model.link <- getlink(model, model.class, model.link)
} else model.class <- NA
#
# CHECK
npar <- if(type=="gam") length(indcoef) else ncol(basis)
if(length(coef)!=npar || length(coef)!=dim(vcov)[1] || any(is.na(coef)) ||
any(is.na(vcov)))
stop("coef/vcov not consistent with basis matrix. See help(crosspred)")
#
##########################################################################
# AT, PREDVAR, PREDLAG AND CENTERING
#
# RANGE
range <- if(type=="gam") range(model$model[[basis$term[1]]]) else
attr(basis,"range")
#
# SET at, predvar AND predlag
at <- mkat(at,from,to,by,range,lag,bylag)
predvar <- if(is.matrix(at)) rownames(at) else at
predlag <- seqlag(lag,bylag)
#
# DEFINE CENTERING VALUE (NULL IF UNCENTERED), AND REMOVE INFO FROM BASIS
cen <- mkcen(cen, type, basis, range)
if(type=="one") attributes(basis)$cen <- NULL
if(type=="cb") attributes(basis)$argvar$cen <- NULL
#
################################################################################
# PREDICTION OF LAG-SPECIFIC EFFECTS
#
# CREATE THE MATRIX OF TRANSFORMED CENTRED VARIABLES (DEPENDENT ON TYPE)
Xpred <- mkXpred(type,basis,at,predvar,predlag,cen)
#
# CREATE LAG-SPECIFIC EFFECTS AND SE
matfit <- matrix(Xpred%*%coef, length(predvar), length(predlag))
matse <- matrix(sqrt(pmax(0,rowSums((Xpred%*%vcov)*Xpred))), length(predvar),
length(predlag))
#
# NAMES
rownames(matfit) <- rownames(matse) <- predvar
colnames(matfit) <- colnames(matse) <- outer("lag",predlag,paste,sep="")
#
################################################################################
# PREDICTION OF OVERALL+CUMULATIVE EFFECTS
#
# RE-CREATE LAGGED VALUES (NB: ONLY LAG INTEGERS)
predlag <- seqlag(lag)
#
# CREATE THE MATRIX OF TRANSFORMED VARIABLES (DEPENDENT ON TYPE)
Xpred <- mkXpred(type,basis,at,predvar,predlag,cen)
#
# CREATE OVERALL AND (OPTIONAL) CUMULATIVE EFFECTS AND SE
Xpredall <- 0
if(cumul) {
cumfit <- cumse <- matrix(0,length(predvar),length(predlag))
}
for (i in seq(length(predlag))) {
ind <- seq(length(predvar)) + length(predvar)*(i-1)
Xpredall <- Xpredall + Xpred[ind,,drop=FALSE]
if(cumul) {
cumfit[, i] <- Xpredall %*% coef
cumse[, i] <- sqrt(pmax(0,rowSums((Xpredall%*%vcov)*Xpredall)))
}
}
allfit <- as.vector(Xpredall %*% coef)
allse <- sqrt(pmax(0,rowSums((Xpredall%*%vcov)*Xpredall)))
#
# NAMES
names(allfit) <- names(allse) <- predvar
if(cumul) {
rownames(cumfit) <- rownames(cumse) <- predvar
colnames(cumfit) <- colnames(cumse) <- outer("lag",seqlag(lag),paste,sep="")
}
#
################################################################################
# CREATE THE OBJECT
#
# INITIAL LIST, THEN ADD COMPONENTS
list <- list(predvar=predvar)
if(!is.null(cen)) list$cen <- cen
list <- c(list, list(lag=lag, bylag=bylag, coefficients=coef, vcov=vcov,
matfit=matfit, matse=matse, allfit=allfit, allse=allse))
if(cumul) list <- c(list, list(cumfit=cumfit, cumse=cumse))
#
# MATRICES AND VECTORS WITH EXPONENTIATED EFFECTS AND CONFIDENCE INTERVALS
z <- qnorm(1-(1-ci.level)/2)
if(!is.null(model.link) && model.link %in% c("log","logit")) {
list$matRRfit <- exp(matfit)
list$matRRlow <- exp(matfit-z*matse)
list$matRRhigh <- exp(matfit+z*matse)
list$allRRfit <- exp(allfit)
list$allRRlow <- exp(allfit-z*allse)
names(list$allRRlow) <- names(allfit)
list$allRRhigh <- exp(allfit+z*allse)
names(list$allRRhigh) <- names(allfit)
if(cumul) {
list$cumRRfit <- exp(cumfit)
list$cumRRlow <- exp(cumfit-z*cumse)
list$cumRRhigh <- exp(cumfit+z*cumse)
}
} else {
list$matlow <- matfit-z*matse
list$mathigh <- matfit+z*matse
list$alllow <- allfit-z*allse
names(list$alllow) <- names(allfit)
list$allhigh <- allfit+z*allse
names(list$allhigh) <- names(allfit)
if(cumul) {
list$cumlow <- cumfit-z*cumse
list$cumhigh <- cumfit+z*cumse
}
}
#
list$ci.level <- ci.level
list$model.class <- model.class
list$model.link <- model.link
#
class(list) <- "crosspred"
#
return(list)
}
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.