Nothing
sensitivityJR <- function(z, s, y, beta0, beta1, phi, Pi, psi,
selection, groupings, ci=0.95,
ci.method=c("analytic", "bootstrap"), na.rm=FALSE,
N.boot=100, interval=c(-100,100),
oneSidedTest = FALSE, twoSidedTest = TRUE,
verbose=getOption("verbose"), isSlaveMode=FALSE)
{
withoutCdfs <- isSlaveMode && !missing(ci.method) && is.null(ci.method)
withoutCi <- ((!isSlaveMode && !missing(ci.method) && ci.method == "") ||
(isSlaveMode && !(!missing(ci.method) &&
!is.null(ci.method) &&
'analytic' %in% ci.method)))
calc.coefs <- function(y, beta, dF, RR, interval) {
coefs <- vector(length(beta), "list")
for(i in seq_along(beta)) {
alphahat <- .calc.alphahat(beta.y=beta[i]*y, dF=dF, C=RR,
interval=interval)
w <- .calc.w(alpha=alphahat, beta.y=beta[i]*y)
coefs[[i]] <- list(alphahat=alphahat, w=w)
}
return(coefs)
}
if(!isSlaveMode) {
## Not running a boot strap mode
## Run error checks on variables.
ErrMsg <- NULL
ErrMsg <- c(.CheckSelection(selection, s),
.CheckGroupings(groupings),
.CheckPhiPiPsi(phi=phi, Pi=Pi, psi=psi),
.CheckLength(z=z, s=s, y=y),
.CheckZ(z, groupings, na.rm=na.rm),
.CheckS(s, na.rm=na.rm),
.CheckY(y, s, selection, na.rm=na.rm))
if(length(ErrMsg) > 0L)
stop(paste(ErrMsg, collapse="\n "))
s <- s == selection
z <- z == groupings[2L]
if(na.rm == TRUE) {
naIndex <- !(is.na(s) | is.na(z) | (s & (is.na(y))))
z <- z[naIndex]
s <- s[naIndex]
y <- y[naIndex]
}
if(any(is.na(z) | is.na(s)))
stop("s, z cannot contain any NA values")
if(any(s & is.na(y)))
stop("selected y values cannot be NA")
ErrMsg <- .CheckPhiPiPsi(phi=phi, Pi=Pi, psi=psi, p0=sum(!z & s)/sum(!z), p1=sum(z & s)/sum(z))
if(length(ErrMsg) > 0L)
stop(paste(ErrMsg, collapse="\n "))
}
if(withoutCi)
ci.method <- NULL
else if(isSlaveMode)
ci.method <- "analytic"
else
ci.method <- sort(unique(match.arg(ci.method, several.ok=TRUE)))
z0.s1 <- !z & s
z1.s1 <- z & s
y0 <- y[z0.s1]
y1 <- y[z1.s1]
N <- length(z)
N0 <- sum(!z)
n0 <- sum(z0.s1)
p0 <- n0/N0
N1 <- sum(z)
n1 <- sum(z1.s1)
p1 <- n1/N1
RR <- p1/p0
VE <- 1 - RR
## Calc Pi because it will be needed later
tmp <- .calcPiPhiPsi(Pi=Pi, phi=phi, psi=psi, p0=p0, p1=p1)
Pi <- tmp$Pi
psi <- tmp$psi
phi <- tmp$phi
sens.var <- tmp$sens.var
Fn0 <- ecdf(y0)
y0.uniq <- knots(Fn0)
F0 <- Fn0(y0.uniq)
dF0 <- diff(c(0, F0))
Fn1 <- ecdf(y1)
y1.uniq <- knots(Fn1)
F1 <- Fn1(y1.uniq)
dF1 <- diff(c(0, F1))
ACE.dim <- c(length(beta0), length(beta1), length(Pi))
ACE.length <- prod(ACE.dim)
ACE.dimnames <- list(format(beta0, trim=TRUE),
format(beta1, trim=TRUE),
format(switch(sens.var,
Pi=Pi,
phi=phi,
psi=psi),
trim=TRUE, digits=4,
drop0trailing=TRUE))
names(ACE.dimnames) <- c("beta0", "beta1", sens.var)
ACE <- array(numeric(ACE.length), dim=ACE.dim, dimnames=ACE.dimnames)
FnAs0.dim <- ACE.dim[-2L]
FnAs0.length <- prod(FnAs0.dim)
FnAs0.dimnames <- ACE.dimnames[-2L]
mu0 <- alphahat0 <- array(numeric(FnAs0.length), dim=FnAs0.dim,
dimnames=FnAs0.dimnames)
FnAs1.dim <- ACE.dim[-1L]
FnAs1.length <- prod(FnAs1.dim)
FnAs1.dimnames <- ACE.dimnames[-1L]
mu1 <- alphahat1 <- array(numeric(FnAs1.length), dim=FnAs1.dim,
dimnames=FnAs1.dimnames)
if(!isSlaveMode) {
FnAs0 <- funArray(dim=FnAs0.dim,
dimnames=FnAs0.dimnames)
FnAs1 <- funArray(dim=FnAs1.dim,
dimnames=FnAs1.dimnames)
}
a0 <- Pi/p0
a1 <- Pi/p1
for(i in seq_along(Pi)) {
if(phi[i] == 1) {
ACE.info <- sensitivityGBH(z=z,s=s,y=y,beta=beta0,
groupings=FALSE,
ci.method=ci.method, isSlaveMode=TRUE)
if(!withoutCdfs) {
alphahat0[,i] <- ACE.info$alphahat
FnAs0[,i] <- ACE.info$FnAs0
alphahat1[,i] <- NA
FnAs1[,i] <- ACE.info$FnAs1
}
ACE[,,i] <- ACE.info$ACE
next
}
if(Pi[i] == 0) {
ACE[,,i] <- NA
alphahat0[,i] <- NA
alphahat1[,i] <- NA
next
}
q0c <- quantile(y0, probs=c(a0[i], 1-a0[i]))
q1c <- quantile(y1, probs=c(a1[i], 1-a1[i]))
for(j in seq_along(beta0)) {
if(is.finite(beta0[j])) {
alphahat0[j,i] <- .calc.alphahat(beta.y=beta0[j]*y0.uniq, dF=dF0,
C=a0[i],
interval=interval)
w0 <- .calc.w(alpha=alphahat0[j,i], beta.y=beta0[j]*y0.uniq)
dFas0 <- dF0*w0/a0[i]
if(!isSlaveMode)
Fas0 <- cumsum(dFas0)
mu0[j,i] <- sum(y0.uniq * dFas0)
} else {
if(beta0[j] == Inf) {
Fas0 <- ifelse(y0.uniq <= q0c[1L] & F0 < a0[i], F0/a0[1], 1)
} else if(beta0[j] == -Inf) {
Fas0 <- ifelse(y0.uniq >= q0c[2L], (F0 - (1 - a0[i]))/a0[i], 0)
} else {
stop("Invalid beta0 value ", beta0[j])
}
alphahat0[j,i] <- NA
mu0[j,i] <- sum(y0.uniq * diff(c(0, Fas0)))
}
if(!isSlaveMode) {
FnAs0[j,i] <- stepfun(y0.uniq, c(0, Fas0))
}
}
for(j in seq_along(beta1)) {
if(is.finite(beta1[j])) {
alphahat1[j,i] <- .calc.alphahat(beta.y=beta1[j]*y1.uniq, dF=dF1,
C=a1[i],
interval=interval)
w1 <- .calc.w(alpha=alphahat1[j,i], beta.y=beta1[j]*y1.uniq)
dFas1 <- dF1*w1/a1[i]
if(!isSlaveMode)
Fas1 <- cumsum(dFas1)
mu1[j,i] <- sum(y1.uniq * dFas1)
} else if(is.infinite(beta1[j])) {
if(beta1[j] == Inf) {
Fas1 <- ifelse(y1.uniq <= q1c[1L] & F1 < a1[i], F1/a1[i], 1L)
} else if(beta1[j] == -Inf) {
Fas1 <- ifelse(y1.uniq >= q1c[2L], (F1 - a1[i])/(a1[i]), 0L)
}
alphahat1[j,i] <- NA
mu1[j,i] <- sum(y1.uniq * diff(c(0L, Fas1)))
} else {
stop("Invalid beta1 value ", beta1[j])
}
if(!isSlaveMode)
FnAs1[j,i] <- stepfun(y1.uniq, c(0, Fas1))
}
ACE[,,i] <- outer(mu0[,i], mu1[,i], function(mu0, mu1) mu1 - mu0)
}
if(withoutCdfs)
return(list(ACE=ACE))
cdfs<-list(beta0=beta0, alphahat0=alphahat0, Fas0=FnAs0,
beta1=beta1, alphahat1=alphahat1, Fas1=FnAs1,
phi=phi, Pi=Pi, psi=psi)
if(withoutCi)
return(c(list(ACE=ACE), cdfs))
if(!isSlaveMode) {
if(twoSidedTest) {
if(ci < 0.5)
ci.probs <- c(ci, 1L) - ci/2L
else
ci.probs <- c(0L, ci) + (1-ci)/2
} else {
ci.probs <- NULL
}
if(oneSidedTest) {
ci.probs <- c(ci.probs, ci)
}
ci.probs <- sort(unique(ci.probs))
ci.probsLen <- length(ci.probs)
ACE.ci.dim <- c(ACE.dim, ci.probsLen, length(ci.method))
ACE.ci.length <- prod(ACE.ci.dim)
ACE.ci.dimnames <- c(ACE.dimnames,
list(ci.prob=sprintf("%s%%",
as.character(ci.probs*100)),
ci.method=ci.method))
ACE.ci <- array(numeric(ACE.ci.length), dim=ACE.ci.dim,
dimnames=ACE.ci.dimnames)
}
ACE.var.dim <- c(ACE.dim, length(ci.method))
ACE.var.length <- prod(ACE.var.dim)
ACE.var.dimnames <- c(ACE.dimnames, list(ci.method=ci.method))
ACE.var <- array(numeric(ACE.var.length), dim=ACE.var.dim,
dimnames=ACE.var.dimnames)
## run bootstrap method
if(any(ci.method == "analytic")) {
if(verbose)
cat("running Analytic")
Omega <- matrix(nrow=6,ncol=6)
for(k in seq_along(Pi)) {
if(phi[k] == 1) {
ACE.var[,, k, "analytic"] <- ACE.info$ACE.var
} else for(i in seq_along(beta0)) {
for(j in seq_along(beta1)) {
U <- rbind((1-z)*(p0 - s),
(z)*(p1 - s),
(1-z)*s*(1/(1+exp(-alphahat0[i,k] - beta0[i]*y)) - Pi[k]/p0),
(z)*s*(1/(1+exp(-alphahat1[j,k] - beta1[j]*y)) - Pi[k]/p1),
(1-z)*s*(mu0[i,k] - y*p0/Pi[k]/(1+exp(-alphahat0[i,k] - beta0[i]*y))),
(z)*s*(mu1[j,k] - y*p1/Pi[k]/(1+exp(-alphahat1[j,k] - beta1[j]*y))))
.sumCrossUpperTri(Omega) <- U
Omega <- Omega / N
Omega[lower.tri(Omega)] <- t(Omega)[lower.tri(Omega)]
Gamma <- matrix(colSums(cbind(1-z,
0,
s*(1-z)*Pi[k]/p0^2,
0,
-s*y*(1-z)/Pi[k]/(exp(-beta0[i]*y-alphahat0[i,k])+1),
0,
0,
z,
0,
s*z*Pi[k]/p1^2,
0,
-s*y*z/Pi[k]/(exp(-beta1[j]*y-alphahat1[j,k])+1),
0,
0,
s*exp(-beta0[i]*y-alphahat0[i,k])*(1-z)/(exp(-beta0[i]*y-alphahat0[i,k])+1)^2,
0,
-p0*s*y*exp(-beta0[i]*y-alphahat0[i,k])*(1-z)/Pi[k]/(exp(-beta0[i]*y-alphahat0[i,k])+1)^2,
0,
0,
0,
0,
s*exp(-beta1[j]*y-alphahat1[j,k])*z/(exp(-beta1[j]*y-alphahat1[j,k])+1)^2,
0,
-p1*s*y*exp(-beta1[j]*y-alphahat1[j,k])*z/Pi[k]/(exp(-beta1[j]*y-alphahat1[j,k])+1)^2,
0,
0,
0,
0,
s*(1-z),
0,
0,
0,
0,
0,
0,
s*z), na.rm=TRUE),
nrow=6,byrow=TRUE) / N
IGamma <- solve(Gamma)
## vartheta <- tcrossprod(IGamma %*% Omega, IGamma) / N
vartheta <- (t(IGamma) %*% Omega %*% IGamma) / N
## cat(all.equal(vartheta, alt), '\n',sep='')
ACE.var[i, j, k, "analytic"] <- vartheta[5, 5] + vartheta[6, 6] - 2 * vartheta[5, 6]
if(verbose) cat(".")
}
}
}
calculateCi <- function(i, norm, ACE, sqrt.ACE.var) {
ACE[i] + norm * sqrt.ACE.var[i]
}
if(!isSlaveMode) {
ACE.ci[,,,,"analytic"] <- outer(seq_along(ACE), qnorm(ci.probs),
FUN=calculateCi, ACE=ACE,
sqrt.ACE.var=sqrt(ACE.var[,,,'analytic']))
}
if(verbose) cat("\n")
}
if(isSlaveMode) {
return(c(list(ACE=ACE, ACE.var=ACE.var), cdfs))
}
if(any(ci.method == "bootstrap")) {
if(verbose) cat("running Bootstrap")
ACE.list <- array(dim=c(ACE.dim,N.boot))
for(i in seq_len(N.boot)) {
new.index <- sample(seq_len(N), N, replace=TRUE)
ACE.list[,,,i] <- Recall(z=z[new.index],s=s[new.index],y=y[new.index],
beta0=beta0, beta1=beta1, psi=psi,
groupings=FALSE,
ci.method=NULL, isSlaveMode=TRUE)$ACE
if(verbose) cat(".")
}
for(k in seq_along(Pi)) {
for(i in seq_along(beta0)) {
for(j in seq_along(beta1)) {
ACE.ci[i,j,k,,'bootstrap'] <- quantile(ACE.list[i,j,k,], probs=ci.probs)
ACE.var[i,j,k,'bootstrap'] <- var(ACE.list[i,j,k,])
}
}
}
if(verbose) cat("\n")
}
return(structure(c(list(ACE=ACE, ACE.ci=ACE.ci, ACE.var=ACE.var),
cdfs),
class=c("sensitivity.2.0d", "sensitivity.0d", "sensitivity"),
N.boot=N.boot,
parameters=list(z0=groupings[1], z1=groupings[2],
selected=selection)))
}
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.