Nothing
R/sand_var.R
In RCTrep: Validation of Estimates of Treatment Effects in Observational Data
Defines functions
sand_bin.modified
## Notation #########################
# single e stands for no augmentation
# ea stands for augmentation
# c represents conservative variance
#####################################
################## Sandwich Variance in binary group case####################################################
#' @importFrom numDeriv jacobian
#' @importFrom stats plogis
sand_bin.modified <- function(z, y, n, ftilt, thetaest, W = NULL, XY = NULL, eest, m0est, m1est, family = "gaussian", offset.e, type = "e", weight.external = NULL) {
p <- ncol(W)
q <- ncol(XY)
# matrix to calculate mu0,mu1 correlation
a <- matrix(0, 2, length(thetaest))
## score function
if (type == "e") {
phi <- function(theta) {
mu0 <- theta[1]
mu1 <- theta[2]
beta <- theta[3:(p + 2)]
e <- plogis(c(W %*% beta))
f1 <- (y - mu0) * (1 - z) * ftilt(e, weight.external) / (1 - e)
f2 <- (y - mu1) * z * ftilt(e, weight.external) / e
f3 <- W * (z - e)
f <- rbind(f1, f2, t(f3))
return(f)
}
a[1, 1] <- 1
a[2, 2] <- 1
} else if (type == "ec") {
phi <- function(theta) {
mu0 <- theta[1]
mu1 <- theta[2]
f1 <- (y - mu0) * (1 - z) * ftilt(eest, weight.external) / (1 - eest)
f2 <- (y - mu1) * z * ftilt(eest, weight.external) / eest
f <- rbind(f1, f2)
return(f)
}
a[1, 1] <- 1
a[2, 2] <- 1
} else if (type == "ea") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
beta <- theta[7:(p + 6)]
gamma0 <- theta[(p + 7):(p + 6 + q)]
gamma1 <- theta[(p + 7 + q):(p + 2 * q + 6)]
e <- plogis(c(W %*% beta))
if (family == "gaussian") {
m0 <- c(XY %*% gamma0)
m1 <- c(XY %*% gamma1)
} else if (family == "binomial") {
m0 <- plogis(c(XY %*% gamma0))
m1 <- plogis(c(XY %*% gamma1))
} else if (family == "poisson") {
m0 <- exp(c(XY %*% gamma0)) * offset.e
m1 <- exp(c(XY %*% gamma1)) * offset.e
}
f1 <- (1 - z) * (y - mu0) * ftilt(e, weight.external) / (1 - e)
f2 <- ftilt(e, weight.external) * (m0 - aug0h)
f3 <- (1 - z) * (m0 - aug0z) * ftilt(e, weight.external) / (1 - e)
f4 <- z * (y - mu1) * ftilt(e, weight.external) / e
f5 <- ftilt(e, weight.external) * (m1 - aug1h)
f6 <- z * (m1 - aug1z) * ftilt(e, weight.external) / e
f7 <- XY * (y - m0) * (1 - z)
f8 <- XY * (y - m1) * z
f9 <- W * (z - e)
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7), t(f8), t(f9))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "ecac") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
f1 <- (1 - z) * (y - mu0) * ftilt(eest, weight.external) / (1 - eest)
f2 <- ftilt(eest, weight.external) * (m0est - aug0h)
f3 <- (1 - z) * (m0est - aug0z) * ftilt(eest, weight.external) / (1 - eest)
f4 <- z * (y - mu1) * ftilt(eest, weight.external) / eest
f5 <- ftilt(eest, weight.external) * (m1est - aug1h)
f6 <- z * (m1est - aug1z) * ftilt(eest, weight.external) / eest
f <- rbind(f1, f2, f3, f4, f5, f6)
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "eca") {
# define the m estimator
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
gamma0 <- theta[7:(6 + q)]
gamma1 <- theta[(7 + q):(2 * q + 6)]
if (family == "gaussian") {
m1 <- c(XY %*% gamma1)
m0 <- c(XY %*% gamma0)
} else if (family == "binomial") {
m1 <- plogis(c(XY %*% gamma1))
m0 <- plogis(c(XY %*% gamma0))
} else if (family == "poisson") {
m1 <- exp(c(XY %*% gamma1)) * offset.e
m0 <- exp(c(XY %*% gamma0)) * offset.e
}
f1 <- (1 - z) * (y - mu0) * ftilt(eest, weight.external) / (1 - eest)
f2 <- ftilt(eest, weight.external) * (m0 - aug0h)
f3 <- (1 - z) * (m0 - aug0z) * ftilt(eest, weight.external) / (1 - eest)
f4 <- z * (y - mu1) * ftilt(eest, weight.external) / eest
f5 <- ftilt(eest, weight.external) * (m1 - aug1h)
f6 <- z * (m1 - aug1z) * ftilt(eest, weight.external) / eest
f7 <- XY * (y - m0) * (1 - z)
f8 <- XY * (y - m1) * z
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7), t(f8))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "eac") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
beta <- theta[7:(p + 6)]
e <- plogis(c(W %*% beta))
f1 <- (1 - z) * (y - mu0) * ftilt(e) / (1 - e)
f2 <- ftilt(e) * (m0est - aug0h)
f3 <- (1 - z) * (m0est - aug0z) * ftilt(e) / (1 - e)
f4 <- z * (y - mu1) * ftilt(e) / e
f5 <- ftilt(e) * (m1est - aug1h)
f6 <- z * (m1est - aug1z) * ftilt(e) / e
f7 <- W * (z - e)
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
}
mphi <- function(theta) {
rowMeans(phi(theta))
}
# define the meat B, covariance operator
Omega <- function(theta) {
phis <- phi(theta)
return(tcrossprod(phis) / n)
}
Atheta <- jacobian(mphi, thetaest)
invAtheta <- solve(Atheta)
# calculate the sandwich variance
Vtmp <- invAtheta %*% Omega(thetaest) %*% t(invAtheta) / n
covmu <- a %*% Vtmp %*% t(a)
return(covmu)
}
################## Sandwich Variance in multiple group case####################################################
# sand_mul <- function(z, y, n, ncate, ftilt, thetaest, W = NULL, XY = NULL, eest, mest, family = "gaussian", offset.e, type = "e") {
# p <- ncol(W)
# q <- ncol(XY)
#
# # matrix to calucalte correlation
# a <- matrix(0, ncate, length(thetaest))
#
# ## score function
# if (type == "e") {
# p <- ncol(W)
# phi <- function(theta) {
# mu <- theta[1:ncate]
# beta <- theta[(ncate + 1):(ncate + (ncate - 1) * p)]
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# wtt <- (1 / e) * ftilt(e)
# fmu <- c()
# fbeta <- c()
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, t(fbeta))
# return(f)
# }
# diag(a) <- 1
# } else if (type == "ec") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# wtt <- (1 / eest) * ftilt(eest, weight.external)
# fmu <- c()
# fbeta <- c()
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# }
#
# f <- rbind(fmu)
# return(f)
# }
# diag(a) <- 1
# } else if (type == "ea") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# beta <- theta[(3 * ncate + 1):(3 * ncate + (ncate - 1) * p)]
# gamma <- theta[(3 * ncate + (ncate - 1) * p + 1):(3 * ncate + (ncate - 1) * p + ncate * q)]
#
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# wtt <- (1 / e) * ftilt(e)
#
# # outcome aug
# m <- c()
# if (family == "gaussian") {
# for (i in 1:ncate) {
# mtmp <- c(XY %*% gamma[((i - 1) * q + 1):(i * q)])
# m <- cbind(m, mtmp)
# }
# } else if (family == "binomial") {
# for (i in 1:ncate) {
# mtmp <- stat::plogis(c(XY %*% gamma[((i - 1) * q + 1):(i * q)]))
# m <- cbind(m, mtmp)
# }
# } else if (family == "poisson") {
# for (i in 1:ncate) {
# mtmp <- exp(c(XY %*% gamma[((i - 1) * q + 1):(i * q)])) * offset.e
# m <- cbind(m, mtmp)
# }
# }
#
# fmu <- c()
# fbeta <- c()
# fgamma <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# fgammatmp <- XY * (y - m[, i]) * (z == i)
# fgamma <- cbind(fgamma, fgammatmp)
# faugztmp <- (m[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- ftilt(e) * (m[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fbeta), t(fgamma))
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "ecac") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# htilt <- ftilt(eest, weight.external)
# wtt <- (1 / eest) * htilt
#
# fmu <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# faugztmp <- (mest[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (mest[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
#
# f <- rbind(fmu, faugz, faugh)
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "eca") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# gamma <- theta[(3 * ncate + 1):(3 * ncate + ncate * q)]
#
# htilt <- ftilt(eest, weight.external)
# wtt <- (1 / eest) * htilt
#
# m <- c()
#
# if (family == "gaussian") {
# for (i in 1:ncate) {
# mtmp <- c(XY %*% gamma[((i - 1) * q + 1):(i * q)])
# m <- cbind(m, mtmp)
# }
# } else if (family == "binomial") {
# for (i in 1:ncate) {
# mtmp <- stat::plogis(c(XY %*% gamma[((i - 1) * q + 1):(i * q)]))
# m <- cbind(m, mtmp)
# }
# } else if (family == "poisson") {
# for (i in 1:ncate) {
# mtmp <- exp(c(XY %*% gamma[((i - 1) * q + 1):(i * q)])) * offset.e
# m <- cbind(m, mtmp)
# }
# }
#
# fmu <- c()
# fgamma <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# fgammatmp <- XY * (y - m[, i]) * (z == i)
# fgamma <- cbind(fgamma, fgammatmp)
# faugztmp <- (m[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (m[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fgamma))
# return(f)
# }
#
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "eac") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# beta <- theta[(3 * ncate + 1):(3 * ncate + (ncate - 1) * p)]
#
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# htilt <- ftilt(e)
# wtt <- (1 / e) * htilt
#
# fmu <- c()
# fbeta <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# faugztmp <- (mest[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (mest[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fbeta))
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# }
#
# mphi <- function(theta) {
# rowMeans(phi(theta))
# }
#
# # define the meat B, covariance operator
# Omega <- function(theta) {
# phis <- phi(theta)
# return(tcrossprod(phis) / n)
# }
#
# Atheta <- jacobian(mphi, thetaest)
# invAtheta <- solve(Atheta)
#
# # calculate the sandwich variance
# Vtmp <- invAtheta %*% Omega(thetaest) %*% t(invAtheta) / n
#
# covmu <- a %*% Vtmp %*% t(a)
#
# return(covmu)
# }
Try the RCTrep package in your browser
Any scripts or data that you put into this service are public.
RCTrep documentation built on Nov. 2, 2023, 6:10 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 sand_bin.modified
## Notation #########################
# single e stands for no augmentation
# ea stands for augmentation
# c represents conservative variance
#####################################
################## Sandwich Variance in binary group case####################################################
#' @importFrom numDeriv jacobian
#' @importFrom stats plogis
sand_bin.modified <- function(z, y, n, ftilt, thetaest, W = NULL, XY = NULL, eest, m0est, m1est, family = "gaussian", offset.e, type = "e", weight.external = NULL) {
p <- ncol(W)
q <- ncol(XY)
# matrix to calculate mu0,mu1 correlation
a <- matrix(0, 2, length(thetaest))
## score function
if (type == "e") {
phi <- function(theta) {
mu0 <- theta[1]
mu1 <- theta[2]
beta <- theta[3:(p + 2)]
e <- plogis(c(W %*% beta))
f1 <- (y - mu0) * (1 - z) * ftilt(e, weight.external) / (1 - e)
f2 <- (y - mu1) * z * ftilt(e, weight.external) / e
f3 <- W * (z - e)
f <- rbind(f1, f2, t(f3))
return(f)
}
a[1, 1] <- 1
a[2, 2] <- 1
} else if (type == "ec") {
phi <- function(theta) {
mu0 <- theta[1]
mu1 <- theta[2]
f1 <- (y - mu0) * (1 - z) * ftilt(eest, weight.external) / (1 - eest)
f2 <- (y - mu1) * z * ftilt(eest, weight.external) / eest
f <- rbind(f1, f2)
return(f)
}
a[1, 1] <- 1
a[2, 2] <- 1
} else if (type == "ea") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
beta <- theta[7:(p + 6)]
gamma0 <- theta[(p + 7):(p + 6 + q)]
gamma1 <- theta[(p + 7 + q):(p + 2 * q + 6)]
e <- plogis(c(W %*% beta))
if (family == "gaussian") {
m0 <- c(XY %*% gamma0)
m1 <- c(XY %*% gamma1)
} else if (family == "binomial") {
m0 <- plogis(c(XY %*% gamma0))
m1 <- plogis(c(XY %*% gamma1))
} else if (family == "poisson") {
m0 <- exp(c(XY %*% gamma0)) * offset.e
m1 <- exp(c(XY %*% gamma1)) * offset.e
}
f1 <- (1 - z) * (y - mu0) * ftilt(e, weight.external) / (1 - e)
f2 <- ftilt(e, weight.external) * (m0 - aug0h)
f3 <- (1 - z) * (m0 - aug0z) * ftilt(e, weight.external) / (1 - e)
f4 <- z * (y - mu1) * ftilt(e, weight.external) / e
f5 <- ftilt(e, weight.external) * (m1 - aug1h)
f6 <- z * (m1 - aug1z) * ftilt(e, weight.external) / e
f7 <- XY * (y - m0) * (1 - z)
f8 <- XY * (y - m1) * z
f9 <- W * (z - e)
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7), t(f8), t(f9))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "ecac") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
f1 <- (1 - z) * (y - mu0) * ftilt(eest, weight.external) / (1 - eest)
f2 <- ftilt(eest, weight.external) * (m0est - aug0h)
f3 <- (1 - z) * (m0est - aug0z) * ftilt(eest, weight.external) / (1 - eest)
f4 <- z * (y - mu1) * ftilt(eest, weight.external) / eest
f5 <- ftilt(eest, weight.external) * (m1est - aug1h)
f6 <- z * (m1est - aug1z) * ftilt(eest, weight.external) / eest
f <- rbind(f1, f2, f3, f4, f5, f6)
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "eca") {
# define the m estimator
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
gamma0 <- theta[7:(6 + q)]
gamma1 <- theta[(7 + q):(2 * q + 6)]
if (family == "gaussian") {
m1 <- c(XY %*% gamma1)
m0 <- c(XY %*% gamma0)
} else if (family == "binomial") {
m1 <- plogis(c(XY %*% gamma1))
m0 <- plogis(c(XY %*% gamma0))
} else if (family == "poisson") {
m1 <- exp(c(XY %*% gamma1)) * offset.e
m0 <- exp(c(XY %*% gamma0)) * offset.e
}
f1 <- (1 - z) * (y - mu0) * ftilt(eest, weight.external) / (1 - eest)
f2 <- ftilt(eest, weight.external) * (m0 - aug0h)
f3 <- (1 - z) * (m0 - aug0z) * ftilt(eest, weight.external) / (1 - eest)
f4 <- z * (y - mu1) * ftilt(eest, weight.external) / eest
f5 <- ftilt(eest, weight.external) * (m1 - aug1h)
f6 <- z * (m1 - aug1z) * ftilt(eest, weight.external) / eest
f7 <- XY * (y - m0) * (1 - z)
f8 <- XY * (y - m1) * z
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7), t(f8))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
} else if (type == "eac") {
phi <- function(theta) {
mu0 <- theta[1]
aug0h <- theta[2]
aug0z <- theta[3]
mu1 <- theta[4]
aug1h <- theta[5]
aug1z <- theta[6]
beta <- theta[7:(p + 6)]
e <- plogis(c(W %*% beta))
f1 <- (1 - z) * (y - mu0) * ftilt(e) / (1 - e)
f2 <- ftilt(e) * (m0est - aug0h)
f3 <- (1 - z) * (m0est - aug0z) * ftilt(e) / (1 - e)
f4 <- z * (y - mu1) * ftilt(e) / e
f5 <- ftilt(e) * (m1est - aug1h)
f6 <- z * (m1est - aug1z) * ftilt(e) / e
f7 <- W * (z - e)
f <- rbind(f1, f2, f3, f4, f5, f6, t(f7))
return(f)
}
a[1, 1:3] <- c(1, 1, -1)
a[2, 4:6] <- c(1, 1, -1)
}
mphi <- function(theta) {
rowMeans(phi(theta))
}
# define the meat B, covariance operator
Omega <- function(theta) {
phis <- phi(theta)
return(tcrossprod(phis) / n)
}
Atheta <- jacobian(mphi, thetaest)
invAtheta <- solve(Atheta)
# calculate the sandwich variance
Vtmp <- invAtheta %*% Omega(thetaest) %*% t(invAtheta) / n
covmu <- a %*% Vtmp %*% t(a)
return(covmu)
}
################## Sandwich Variance in multiple group case####################################################
# sand_mul <- function(z, y, n, ncate, ftilt, thetaest, W = NULL, XY = NULL, eest, mest, family = "gaussian", offset.e, type = "e") {
# p <- ncol(W)
# q <- ncol(XY)
#
# # matrix to calucalte correlation
# a <- matrix(0, ncate, length(thetaest))
#
# ## score function
# if (type == "e") {
# p <- ncol(W)
# phi <- function(theta) {
# mu <- theta[1:ncate]
# beta <- theta[(ncate + 1):(ncate + (ncate - 1) * p)]
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# wtt <- (1 / e) * ftilt(e)
# fmu <- c()
# fbeta <- c()
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, t(fbeta))
# return(f)
# }
# diag(a) <- 1
# } else if (type == "ec") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# wtt <- (1 / eest) * ftilt(eest, weight.external)
# fmu <- c()
# fbeta <- c()
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# }
#
# f <- rbind(fmu)
# return(f)
# }
# diag(a) <- 1
# } else if (type == "ea") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# beta <- theta[(3 * ncate + 1):(3 * ncate + (ncate - 1) * p)]
# gamma <- theta[(3 * ncate + (ncate - 1) * p + 1):(3 * ncate + (ncate - 1) * p + ncate * q)]
#
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# wtt <- (1 / e) * ftilt(e)
#
# # outcome aug
# m <- c()
# if (family == "gaussian") {
# for (i in 1:ncate) {
# mtmp <- c(XY %*% gamma[((i - 1) * q + 1):(i * q)])
# m <- cbind(m, mtmp)
# }
# } else if (family == "binomial") {
# for (i in 1:ncate) {
# mtmp <- stat::plogis(c(XY %*% gamma[((i - 1) * q + 1):(i * q)]))
# m <- cbind(m, mtmp)
# }
# } else if (family == "poisson") {
# for (i in 1:ncate) {
# mtmp <- exp(c(XY %*% gamma[((i - 1) * q + 1):(i * q)])) * offset.e
# m <- cbind(m, mtmp)
# }
# }
#
# fmu <- c()
# fbeta <- c()
# fgamma <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# fgammatmp <- XY * (y - m[, i]) * (z == i)
# fgamma <- cbind(fgamma, fgammatmp)
# faugztmp <- (m[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- ftilt(e) * (m[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fbeta), t(fgamma))
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "ecac") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# htilt <- ftilt(eest, weight.external)
# wtt <- (1 / eest) * htilt
#
# fmu <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# faugztmp <- (mest[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (mest[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
#
# f <- rbind(fmu, faugz, faugh)
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "eca") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# gamma <- theta[(3 * ncate + 1):(3 * ncate + ncate * q)]
#
# htilt <- ftilt(eest, weight.external)
# wtt <- (1 / eest) * htilt
#
# m <- c()
#
# if (family == "gaussian") {
# for (i in 1:ncate) {
# mtmp <- c(XY %*% gamma[((i - 1) * q + 1):(i * q)])
# m <- cbind(m, mtmp)
# }
# } else if (family == "binomial") {
# for (i in 1:ncate) {
# mtmp <- stat::plogis(c(XY %*% gamma[((i - 1) * q + 1):(i * q)]))
# m <- cbind(m, mtmp)
# }
# } else if (family == "poisson") {
# for (i in 1:ncate) {
# mtmp <- exp(c(XY %*% gamma[((i - 1) * q + 1):(i * q)])) * offset.e
# m <- cbind(m, mtmp)
# }
# }
#
# fmu <- c()
# fgamma <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# fgammatmp <- XY * (y - m[, i]) * (z == i)
# fgamma <- cbind(fgamma, fgammatmp)
# faugztmp <- (m[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (m[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fgamma))
# return(f)
# }
#
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# } else if (type == "eac") {
# phi <- function(theta) {
# mu <- theta[1:ncate]
# augz <- theta[(ncate + 1):(2 * ncate)]
# augh <- theta[(2 * ncate + 1):(3 * ncate)]
# beta <- theta[(3 * ncate + 1):(3 * ncate + (ncate - 1) * p)]
#
# e <- rep(1, n)
# for (i in 1:(ncate - 1)) {
# etmp <- exp(W %*% beta[((i - 1) * p + 1):(i * p)])
# e <- cbind(e, etmp)
# }
# e <- e / (apply(e, 1, sum))
# htilt <- ftilt(e)
# wtt <- (1 / e) * htilt
#
# fmu <- c()
# fbeta <- c()
# faugz <- c()
# faugh <- c()
#
# for (i in 1:ncate) {
# fmutmp <- (y - mu[i]) * wtt[, i] * (z == i)
# fmu <- rbind(fmu, fmutmp)
# faugztmp <- (mest[, i] - augz[i]) * wtt[, i] * (z == i)
# faughtmp <- htilt * (mest[, i] - augh[i])
# faugz <- rbind(faugz, faugztmp)
# faugh <- rbind(faugh, faughtmp)
# }
# for (i in 2:ncate) {
# fbetatmp <- ((z == i) - e[, i]) * W
# fbeta <- cbind(fbeta, fbetatmp)
# }
#
# f <- rbind(fmu, faugz, faugh, t(fbeta))
# return(f)
# }
# for (i in 1:ncate) {
# a[i, c(i, 2 * ncate + i)] <- 1
# a[i, ncate + i] <- -1
# }
# }
#
# mphi <- function(theta) {
# rowMeans(phi(theta))
# }
#
# # define the meat B, covariance operator
# Omega <- function(theta) {
# phis <- phi(theta)
# return(tcrossprod(phis) / n)
# }
#
# Atheta <- jacobian(mphi, thetaest)
# invAtheta <- solve(Atheta)
#
# # calculate the sandwich variance
# Vtmp <- invAtheta %*% Omega(thetaest) %*% t(invAtheta) / n
#
# covmu <- a %*% Vtmp %*% t(a)
#
# return(covmu)
# }
Try the RCTrep 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.