#Calculate attributable fraction
### attrdl function
attrdl <- function(x,basis,cases,model=NULL,coef=NULL,vcov=NULL,type="af",
dir="back",tot=TRUE,cen,range=NULL,sim=FALSE,nsim=5000) {
#
# CHECK VERSION OF THE DLNM PACKAGE
if(packageVersion("dlnm")<"2.2.0")
stop("update dlnm package to version >= 2.2.0")
## EXTRACT NAME AND CHECK type AND dir
name <- deparse(substitute(basis))
type <- match.arg(type,c("an","af"))
dir <- match.arg(dir,c("back","forw"))
#
# DEFINE CENTERING
if(missing(cen) && is.null(cen <- attr(basis,"argvar")$cen))
stop("'cen' must be provided")
if(!is.numeric(cen) && length(cen)>1L) stop("'cen' must be a numeric scalar")
attributes(basis)$argvar$cen <- NULL
## SELECT RANGE (FORCE TO CENTERING VALUE OTHERWISE, MEANING NULL RISK)
if(!is.null(range)) x[x<range[1]|x>range[2]] <- cen
## COMPUTE THE MATRIX OF
# - LAGGED EXPOSURES IF dir="back"
# - CONSTANT EXPOSURES ALONG LAGS IF dir="forw"
lag <- attr(basis,"lag")
if(NCOL(x)==1L) {
at <- if(dir=="back") tsModel:::Lag(x,seq(lag[1],lag[2])) else
matrix(rep(x,diff(lag)+1),length(x))
} else {
if(dir=="forw") stop("'x' must be a vector when dir='forw'")
if(ncol(at <- x)!=diff(lag)+1)
stop("dimension of 'x' not compatible with 'basis'")
}
#
# NUMBER USED FOR THE CONTRIBUTION AT EACH TIME IN FORWARD TYPE
# - IF cases PROVIDED AS A MATRIX, TAKE THE ROW AVERAGE
# - IF PROVIDED AS A TIME SERIES, COMPUTE THE FORWARD MOVING AVERAGE
# - THIS EXCLUDES MISSING ACCORDINGLY
# ALSO COMPUTE THE DENOMINATOR TO BE USED BELOW
if(NROW(cases)!=NROW(at)) stop("'x' and 'cases' not consistent")
if(NCOL(cases)>1L) {
if(dir=="back") stop("'cases' must be a vector if dir='back'")
if(ncol(cases)!=diff(lag)+1) stop("dimension of 'cases' not compatible")
den <- sum(rowMeans(cases,na.rm=TRUE),na.rm=TRUE)
cases <- rowMeans(cases)
} else {
den <- sum(cases,na.rm=TRUE)
if(dir=="forw")
cases <- rowMeans(as.matrix(tsModel:::Lag(cases,-seq(lag[1],lag[2]))))
}
#
# EXTRACT COEF AND VCOV IF MODEL IS PROVIDED
if(!is.null(model)) {
cond <- paste0(name,"[[:print:]]*v[0-9]{1,2}\\.l[0-9]{1,2}")
if(ncol(basis)==1L) cond <- name
model.class <- class(model)
coef <- dlnm:::getcoef(model,model.class)
ind <- grep(cond,names(coef))
coef <- coef[ind]
vcov <- dlnm:::getvcov(model,model.class)[ind,ind,drop=FALSE]
model.link <- dlnm:::getlink(model,model.class)
if(model.link!="log") stop("'model' must have a log link function")
}
#
# IF REDUCED ESTIMATES ARE PROVIDED
typebasis <- ifelse(length(coef)!=ncol(basis),"one","cb")
#
# PREPARE THE ARGUMENTS FOR TH BASIS TRANSFORMATION
predvar <- if(typebasis=="one") x else seq(NROW(at))
predlag <- if(typebasis=="one") 0 else dlnm:::seqlag(lag)
#
# CREATE THE MATRIX OF TRANSFORMED CENTRED VARIABLES (DEPENDENT ON typebasis)
if(typebasis=="cb") {
Xpred <- dlnm:::mkXpred(typebasis,basis,at,predvar,predlag,cen)
Xpredall <- 0
for (i in seq(length(predlag))) {
ind <- seq(length(predvar))+length(predvar)*(i-1)
Xpredall <- Xpredall + Xpred[ind,,drop=FALSE]
}
} else {
basis <- do.call(onebasis,c(list(x=x),attr(basis,"argvar")))
Xpredall <- dlnm:::mkXpred(typebasis,basis,x,predvar,predlag,cen)
}
#
# CHECK DIMENSIONS
if(length(coef)!=ncol(Xpredall))
stop("arguments 'basis' do not match 'model' or 'coef'-'vcov'")
if(any(dim(vcov)!=c(length(coef),length(coef))))
stop("arguments 'coef' and 'vcov' do no match")
if(typebasis=="one" && dir=="back")
stop("only dir='forw' allowed for reduced estimates")
#
# COMPUTE AF AND AN
af <- 1-exp(-drop(as.matrix(Xpredall%*%coef)))
an <- af*cases
#
# TOTAL
# - SELECT NON-MISSING OBS CONTRIBUTING TO COMPUTATION
# - DERIVE TOTAL AF
# - COMPUTE TOTAL AN WITH ADJUSTED DENOMINATOR (OBSERVED TOTAL NUMBER)
if(tot) {
isna <- is.na(an)
af <- sum(an[!isna])/sum(cases[!isna])
an <- af*den
}
#
# EMPIRICAL CONFIDENCE INTERVALS
if(!tot && sim) {
sim <- FALSE
warning("simulation samples only returned for tot=T")
}
if(sim) {
# SAMPLE COEF
k <- length(coef)
eigen <- eigen(vcov)
X <- matrix(rnorm(length(coef)*nsim),nsim)
coefsim <- coef + eigen$vectors %*% diag(sqrt(eigen$values),k) %*% t(X)
# RUN THE LOOP
# pre_afsim <- (1 - exp(- Xpredall %*% coefsim)) * cases # a matrix
# afsim <- colSums(pre_afsim,na.rm=TRUE) / sum(cases[!isna],na.rm=TRUE)
afsim <- apply(coefsim,2, function(coefi) {
ani <- (1-exp(-drop(Xpredall%*%coefi)))*cases
sum(ani[!is.na(ani)])/sum(cases[!is.na(ani)])
})
ansim <- afsim*den
}
#
res <- if(sim) {
if(type=="an") ansim else afsim
} else {
if(type=="an") an else af
}
#
return(res)
}
##Build up the matrix and calculate the AF.
matsim <- matrix(NA,1,3,dimnames=list("shanghai",c("glob","cold","heat")))
nsim <- 1000
matsim[1,"glob"] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28)
matsim[1,"cold"] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28, range=c(-100,28))
matsim[1,"heat"] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28, range=c(28,100))
arraysim <- array(NA, dim = c(1, 3, nsim), dimnames = list("shanghai", c("glob", "cold", "heat")))
arraysim[1,"glob",] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28,sim=T,nsim=nsim)
arraysim[1,"cold",] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28,
range=c(-100,28),sim=T,nsim=nsim)
arraysim[1,"heat",] <- attrdl(data$temp,cb,data$total,coef=predcum$coefficients, vcov=predcum$vcov,type="an",dir="forw",cen=28,
range=c(28,100),sim=T,nsim=nsim)
## Calculating total deaths
tottotal <- sum(data$total,na.rm=T)
ancity <- matsim
ancitylow <- apply(arraysim,c(1,2),quantile,0.025)
ancityhigh <- apply(arraysim,c(1,2),quantile,0.975)
rownames(ancity) <- rownames(ancitylow) <- rownames(ancityhigh) <- "shanghai"
# Calculating AF
afcity <- ancity/tottotal*100
afcitylow <- ancitylow/tottotal*100
afcityhigh <- ancityhigh/tottotal*100
print(afcity)
