Nothing
BAPC <- function(APCList, predict = list(npredict = 0, retro = TRUE),
model=list(age=list(model="rw2", prior="loggamma", param=c(1, 0.00005), initial=4, scale.model=FALSE),
period=list(include=TRUE, model="rw2", prior="loggamma", param=c(1, 0.00005), initial=4, scale.model=FALSE),
cohort=list(include=TRUE, model="rw2", prior="loggamma", param=c(1, 0.00005), initial=4, scale.model=FALSE),
overdis=list(include=TRUE, model="iid", prior="loggamma", param=c(1, 0.005), initial=4)),
secondDiff=FALSE, stdweight=NULL, verbose=FALSE){
if(!is.null(stdweight)){
if(length(stdweight) != nage(APCList))
stop("\n\n\tThe weight vector for age-standardization is not of correct length!\n")
}
if(!is.null(stdweight) & (sum(stdweight)!= 1)){
stdweight=stdweight/sum(stdweight)
cat("\n\n\tCaution: Your age-specific weights were normalised so that they sum to one.\n")
}
if(predict$npredict < 0){
stop("\n\n\tThe argument npredict must be an integer larger or equal to zero\n")
}
if(!is.logical(secondDiff)){
stop("\n\n\tThe argument secondDiff must be either TRUE or FALSE\n")
}
if(!is.logical(verbose)){
stop("\n\n\tThe argument verbose must be either TRUE or FALSE\n")
}
model$age <- .checkmodel(model$age)
model$period <- .checkmodel(model$period)
model$cohort <- .checkmodel(model$cohort)
model$overdis <- .checkoverdis(model$overdis)
APCList@npred <- predict$npredict
# number of a data.inla <- data.frame(y=y, n=n, i=i, j=j, k=k, z=z)
I <- nage(APCList)
# number of periods
J <- nperiod(APCList)
# number of cohorts
K <- ncohort(APCList)
y <- APCList@epi
# retrospective prediction
if(predict$npredict > 0){
if(predict$retro){
y <- epi(APCList)
y[(J-predict$npredict + 1):J, ] <- NA
y <- c(y)
}
}
n <- APCList@pyrs
# age index
i <- ageindex(APCList)
# period index
j <- periodindex(APCList)
# cohort index
k <- cohortindex(APCList)
# overdispersion parameter
z <- overdisindex(APCList)
# data frame for inla
data.inla <- data.frame(y=y, n=n, i=i, j=j, k=k, z=z)
# formula object
formula <- paste("y ~ f(i, model=\"", model$age$model , "\",
hyper=list(prec=list(prior=\"", model$age$prior, "\", param=c(", model$age$param[1], ",", model$age$param[2], "),
initial=", model$age$initial, ")), constr=T, scale.model=", model$age$scale.model, ")", sep="")
if(model$period$include){
if(model$period$model != "drift"){
formula <- paste(formula, paste("f(j, model=\"", model$period$model , "\",
hyper=list(prec=list(prior=\"", model$period$prior, "\", param=c(", model$period$param[1], ",", model$period$param[2], "),
initial=", model$period$initial, ")), scale.model=", model$period$scale.model, ")", sep=""), collapse="+", sep="+")
} else {
formula <- paste(formula, "j", collapse="+", sep="+")
}
}
if(model$cohort$include){
if(model$cohort$model != "drift"){
formula <- paste(formula, paste("f(k, model=\"", model$cohort$model , "\",
hyper=list(prec=list(prior=\"", model$cohort$prior, "\", param=c(", model$cohort$param[1], ",", model$cohort$param[2], "),
initial=", model$cohort$initial, ")), scale.model=", model$cohort$scale.model, ")", sep=""), collapse="+", sep="+")
} else {
formula <- paste(formula, "k", collapse="+", sep="+")
}
}
if(model$overdis$include){
formula <- paste(formula, paste("f(z, model=\"", model$overdis$model , "\",
hyper=list(prec=list(prior=\"", model$overdis$prior, "\", param=c(", model$overdis$param[1], ",", model$overdis$param[2], "),
initial=", model$overdis$initial, ")))", sep=""), collapse="+", sep="+")
}
if(!is.null(stdweight)){
# make linear combinations which are the nPred linear predictors
len <- length(y)
p <- matrix(NA, len, len)
diag(p) <- 1
lincomb <- INLA::inla.make.lincombs(Predictor=p)
names(lincomb) <- unlist(strsplit(sprintf("lc%.6f", 1:len/10000), ".", fixed=T))[seq(2, 2*len, 2)]
config <- TRUE
} else {
lincomb <- NULL
config <- TRUE
}
if(secondDiff){
# define linear combinations for the second differences
# AGE
for(l in 3:I)
{
idx <- rep(NA, I)
idx[c(l, l-1, l-2)] <- c(1,-2,1)
lc <- INLA::inla.make.lincomb(i=idx)
names(lc) <- paste("diff_i.", .num(l, width=4), sep="")
lincomb <- c(lincomb, lc)
}
# PERIOD
for(l in 3:J)
{
idx <- rep(NA, J)
idx[c(l, l-1, l-2)] <- c(1,-2, 1)
lc <- INLA::inla.make.lincomb(j=idx)
names(lc) <- paste("diff_j.", .num(l, width=4), sep="")
lincomb <- c(lincomb, lc)
}
# COHORT
for(l in 3:K)
{
idx <- rep(NA, K)
idx[c(l, l-1, l-2)] <- c(1,-2, 1)
lc <- INLA::inla.make.lincomb(k=idx)
names(lc) <- paste("diff_k.", .num(l, width=4), sep="")
lincomb <- c(lincomb, lc)
}
}
res <- INLA::inla(as.formula(formula), family="Poisson", data=data.inla, E=n,
control.predictor=list(compute=TRUE, link=1),
lincomb=lincomb,
#quantiles=c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975),
control.compute=list(config=config, dic=TRUE, cpo=TRUE,
return.marginals.predictor=TRUE),
control.inla=list(#lincomb.derived.only=TRUE,
lincomb.derived.correlation.matrix=TRUE),
verbose=verbose)
inlares(APCList) <- res
# lambda=sapply(1:(I*J), function(i){INLA::inla.emarginal(function(x){exp(x)}, res$marginals.linear.predictor[[i]])})
# lambda2=sapply(1:(I*J), function(i){INLA::inla.emarginal(function(x){exp(x)^2}, res$marginals.linear.predictor[[i]])})
#
# # variance of lambda
# slambda = lambda2 - lambda^2
lambda = res$summary.fitted.values$mean
# variance of lambda
slambda = (res$summary.fitted.values$sd)^2
# mean of the predictive distribution
mu = n*lambda
# variance of the predictive distribution
sig = mu + n^2*slambda
mu_m = matrix(mu, ncol=I, nrow=J, byrow=F)
## lci_m = matrix(mu-1.96*sqrt(sig), ncol=I, nrow=J, byrow=F)
# lci_m = matrix(qnorm(0.025, mean=mu, sd=sqrt(sig)), ncol=I, nrow=J, byrow=F)
# do not allow for negative projection rates and set them to zero instead
#lci_m[lci_m < 0] = 0
##uci_m = matrix(mu+1.96*sqrt(sig), ncol=I, nrow=J, byrow=F)
#uci_m = matrix(qnorm(0.975, mean=mu, sd=sqrt(sig)), ncol=I, nrow=J, byrow=F)
sd_m = matrix(sqrt(sig), ncol=I, nrow=J, byrow=F)
#inlaq <- t(sapply(1:(I*J), function(i){INLA::inla.qmarginal(c(0.025, 0.975), res$marginals.linear.predictor[[i]])}))
#pre <- exp(inlaq)
#pre <- cbind("0.025Q"=pre[,1], "mean"=lambda, "0.975Q"=pre[,2], "sd"=sqrt(slambda))
pre <- cbind("mean"=lambda, "sd"=sqrt(slambda))
plab=periodlabels(APCList)
agespec.rate(APCList) <- lapply(1:I, function(m){tmp=pre[((m-1)*J+1):(m*J),]
rownames(tmp)=plab
return(tmp)})
names(agespec.rate(APCList)) <- agelabels(APCList)
#agespec.proj(APCList) <- lapply(1:I, function(m){tmp=cbind("0.025Q"=lci_m[,m], "mean"=mu_m[,m], "0.975Q"=uci_m[,m], "sd"=sd_m[,m]);
# rownames(tmp)=plab
# return(tmp)})
agespec.proj(APCList) <- lapply(1:I, function(m){tmp=cbind("mean"=mu_m[,m], "sd"=sd_m[,m]);
rownames(tmp)=plab
return(tmp)})
names(agespec.proj(APCList)) <- agelabels(APCList)
if(!is.null(stdweight)){
stdweight(APCList) <- stdweight
my.wm <- matrix(rep(stdweight, each=J), byrow=F, nrow=J)
stdobs(APCList)=rowSums(my.wm * epi(APCList)/pyrs(APCList), na.rm=T)*rowSums(pyrs(APCList))
idx.sort <- sort(res$summary.lincomb.derived$ID[1:(I*J)], index.return=T)$ix
lc <- res$summary.lincomb.derived[idx.sort,]
mcor <- res$misc$lincomb.derived.correlation.matrix[1:(I*J), 1:(I*J)]
if(sum(colnames(mcor)!=rownames(mcor)) > 0){
message("WARNING")}
mcor.idx <- sort(as.numeric(colnames(mcor)), index.return=T)$ix
mcor <- mcor[mcor.idx,]
mcor <- mcor[, mcor.idx]
if(sum(colnames(mcor) != rownames(lc)) > 0){
message("WARNING")}
# get the covariance matrix for the linear combinations
mcov <- .cor2cov(mcor, lc$sd)
# use the multivariate delta rule to get the covariance matrix
# for exp(linear combinations)
D_matrix <- diag(exp(lc$mean))
mcov_exp <- D_matrix %*% mcov %*% t(D_matrix)
w_matrix <- matrix(0, ncol=(I*J), nrow=J)
for(m in 1:J){
w_matrix[m, j==m] <- stdweight
}
## get the variances of the age standardize quantities
my_sd <- sqrt(diag(w_matrix %*% mcov_exp %*% t(w_matrix)))
mi <- res$marginals.lincomb.derived[idx.sort]
#new_mean <- matrix(sapply(1:(I*J), function(u){inla.emarginal(function(y){exp(y)}, mi[[u]])}), nrow=J, byrow=F)
new_mean <- exp(lc$mean)
my.wm <- matrix(rep(stdweight, each=J), byrow=F, nrow=J)
new_mean <- rowSums(my.wm*new_mean)
##new_l <- (new_mean - 1.96*my_sd)
##new_u <- (new_mean + 1.96*my_sd)
#new_l <- qnorm(0.025, mean=new_mean, sd=my_sd)
#new_u <- qnorm(0.975, mean=new_mean, sd=my_sd)
#tmp = cbind("0.025Q"=new_l, "mean"=new_mean, "0.975Q"=new_u, "sd"=my_sd)
tmp = cbind("mean"=new_mean, "sd"=my_sd)
rownames(tmp) = plab
agestd.rate(APCList) <- tmp
agg.n <- rowSums(pyrs(APCList))
agg.mean <- agg.n * new_mean
agg.std <- sqrt(agg.mean + agg.n^2*my_sd^2)
#tmp = cbind("0.025Q"=agg.mean-1.96*agg.std, "mean"=agg.mean,
# "0.975Q"=agg.mean+1.96*agg.std, "sd"=agg.std)
tmp = cbind("mean"=agg.mean, "sd"=agg.std)
rownames(tmp) = plab
agestd.proj(APCList) <- tmp
}
return(APCList)
}
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.