Nothing
R/sandwich_ex_log.R
In RegCalReliab: Regression Calibration Using Reliability Studies
Defines functions
sandwich_estimator_ex_log
#' Sandwich Variance Estimator for External Logistic Regression Calibration
#'
#' \code{sandwich_estimator_ex_log()} computes robust (sandwich) standard errors
#' for regression calibration in logistic regression using external reliability data.
#' It propagates uncertainty from estimating the calibration (variance components)
#' into the final coefficient variances, yielding valid confidence intervals and
#' p-values even when measurement error is present.
#'
#' @param xhat Numeric matrix of calibrated exposures (\eqn{n_m \times t}),
#' typically from \code{\link{reg_calibration_ex_log}}.
#' @param z.main.std Numeric matrix of standardized main-study exposures
#' (\eqn{n_m \times t}).
#' @param z.rep.std Named list of standardized replicate measurements from the
#' reliability study. Each list element is a matrix of dimension
#' \eqn{n_r \times r_i}.
#' @param r Integer vector of replicate counts for reliability-study subjects
#' (length \eqn{n_m + n_r}), with 1 for main-study subjects.
#' @param Y Binary (0/1) outcome vector of length \eqn{n_m}.
#' @param indicator Binary vector of length \eqn{n_m + n_r} indicating main-study
#' (1) vs. reliability-study (0) subjects.
#' @param v12star Calibration matrix (\eqn{\Sigma_X \Sigma_Z^{-1}\Sigma_{XZ}})
#' used to transform observed measurements to their calibrated values.
#' @param beta.fit2 Numeric vector of logistic regression coefficients from the
#' calibrated model (\code{fit2}).
#' @param W.main.std Optional numeric matrix of standardized error-free
#' covariates (\eqn{n_m \times q}).
#' @param sigma Within-person variance-covariance matrix of the exposures.
#' @param sigmaz Total variance-covariance matrix of the observed exposures
#' (and covariates, if applicable).
#' @param sigmawithin Estimated within-person covariance matrix (same dimension
#' as \code{sigmaz}).
#' @param sdz Numeric vector of standard deviations of the unstandardized exposures.
#' @param sdw Optional numeric vector of standard deviations of the unstandardized covariates.
#' @param muz Numeric vector of means of the unstandardized exposures.
#' @param muw Optional numeric vector of means of the unstandardized covariates.
#' @param fit2 The \code{glm} object returned by \code{reg_calibration_ex_log()}.
#' @param zbar Matrix of subject-level means of replicate measurements
#' in the reliability study.
#'
#' @return A list with one component:
#' \describe{
#' \item{\code{Sandwich Corrected estimates}}{Matrix of regression calibration
#' estimates with sandwich-corrected standard errors, z-values,
#' p-values, odds ratios, and 95\% confidence intervals on the
#' original scale.}
#' }
#'
#' @details
#' This function implements the three-block estimating equation described in
#' Carroll et al. (2006) and White (1980):
#' \enumerate{
#' \item \strong{Variance components}: estimates total and within-person covariance.
#' \item \strong{Calibration}: computes \eqn{X^\text{hat}} using the
#' estimated variance components.
#' \item \strong{Outcome model}: fits a logistic regression and applies a
#' sandwich (robust) variance formula to combine the uncertainty from
#' all stages.
#' }
#'
#' @examples
#' set.seed(123)
#' # Simulate main-study data
#' z.main <- matrix(rnorm(60), ncol = 1)
#' colnames(z.main) <- "sbp"
#' Y <- rbinom(60, 1, plogis(0.3 * z.main))
#'
#' # Simulate external reliability data (30 subjects, 2 replicates)
#' z.rep <- list(sbp = matrix(rnorm(30 * 2), nrow = 30))
#' zbar <- sapply(z.rep, rowMeans)
#'
#' # Prepare standardization inputs
#' sdz <- apply(z.main, 2, sd)
#' z.main.std <- scale(z.main)
#' z.rep.std <- list(sbp = scale(z.rep$sbp))
#' r <- c(rep(1, 60), rep(2, 30))
#' indicator <- c(rep(1, 60), rep(0, 30))
#'
#' # (In practice, run reg_calibration_ex_log() first to get xhat, fit2, etc.)
#' # Here we show a simplified call using mock objects:
#' # sandwich_estimator_ex_log(xhat, z.main.std, z.rep.std, r, Y, indicator,
#' # v12star, beta.fit2, sigma, sigmaz, sigmawithin, sdz, NULL, muz, NULL, fit2, zbar)
#'
#' @references
#' Carroll RJ, Ruppert D, Stefanski LA, Crainiceanu CM. \emph{Measurement Error
#' in Nonlinear Models: A Modern Perspective}. Chapman & Hall/CRC, 2006.
#' White H. A heteroskedasticity-consistent covariance matrix estimator and a
#' direct test for heteroskedasticity. \emph{Econometrica}. 1980;48(4):817–838.
#'
#' @noRd
sandwich_estimator_ex_log = function(xhat,z.main.std,z.rep.std,r,Y,indicator, v12star,beta.fit2,W.main.std = NULL,sigma,
sigmaz,sigmawithin,sdz,sdw,muz,muw,fit2,zbar){
# -----------------------------------------------
# 0) Basic dimensions
# -----------------------------------------------
nm = nrow(z.main.std) # number of subjects in the main study
nr = sum(indicator == 0)
n = nm+nr # total number of subjects (main + reliability)
t = ncol(z.main.std) # number of variables in z.main.std
r0 = r[indicator==0]
muz_std <- as.numeric(colMeans(z.main.std)) # length t
if (!is.null(W.main.std)) muw_std <- as.numeric(colMeans(W.main.std)) # length q
if(is.null(W.main.std)){
# Compute fitted probabilities p from the corrected (calibrated) logistic model
p = as.vector(exp(beta.fit2 %*% t(cbind(1,xhat)))/(1+exp(beta.fit2 %*% t(cbind(1,xhat)))))
# Compute the inverse of the estimated covariance matrix of Z, sigmaz
sigmaz.inv = solve(sigmaz)
# Create matrix Z as the transpose of the main study data, from nm x t to t x nm
Z = t(cbind(z.main.std))
# Initialize two t x t zero matrices, they will be used later
# to construct off–diagonal derivative matrices.
m = matrix(0,nrow = t,ncol = t)
c = matrix(0,nrow = t,ncol = t)
##sigmaX
###diag
# ddv.x and ddm.x are essentially indicator matrices with a 1 in the (x,x) entry. S and T
ddv.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
ddm.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
# the matrix of partial derivatives (with respect to the (x,x) diagonal element of sigmaz,) of the calibration
# Diagonal Partials wrt sigmaz
db.x = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%ddm.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
# Build off–diagonal derivative matrices
# For each x from 1 to t-1, for each y from x+1 to t, update matrix c
if(t>1){
odv.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y]<-c[x,y]+1
c[y,x]<-c[y,x]+1
c
},simplify = FALSE),simplify = FALSE)
odm.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = FALSE),simplify = F)
# Compute the off-diagonal derivative contributions
ob.x<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
}
# Off‐Diagonal Partials wrt sigmaz
##sigma (the within–group variance)
#Partials wrt Another Variance W
###diag
db.0 = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
# correspond to partial derivatives wrt the diagonal of the other covariance matrix sigmaZW
###off-diag
if(t>1){ob.0 = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)}
# correspond to partial derivatives wrt the off-diagonal of the other covariance matrix sigmaZW
m = (t)*(t+1)/2+(t*(t+1))/2 #number of the covariance estimates
s = t+1 #number of beta estimates (regression parameters)
# If t > 1, include both diagonal and off-diagonal pieces otherwise, only the diagonal pieces
# When t=1 we have only one “error‐prone” covariate—so there are no off‐diagonal elements
# (1×1 covariance matrix has no off‐diagonal entries)
b = if(t>1){
list(db.x,ob.x,db.0,ob.0)
}else list(db.x,db.0)
# Reshape the list b into a list of m matrices
# For each index from 1 to m, extract a matrix with nm rows and t columns from the unlisted b
b = sapply(1:m,function(x) matrix(unlist(b),nrow=nm)[,(t*x-t+1):(t*x)],simplify = FALSE)
d = as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)])) / nm
bstar = sapply(b,function(x) -rowSums(x*d)) # - sum_j b_{ij} * d_{ij}
a = matrix(NA,ncol=m,nrow=s)
a[1,] = colSums(bstar)
a[2:(t+1),] = t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar))) #correspond to A2
w = diag(p * (1 - p))
Astar = t(cbind(1,xhat)) %*% w %*% cbind(1,xhat) /nm # correspond to A3
A = rbind(cbind(diag(rep(-1,m)),
matrix(0,nrow=m,ncol=s)),cbind(a,Astar)) # building up A2 and A3 together
A = rbind(cbind(matrix(diag(rep(-1,t),nrow=t),t),matrix(0,nrow=t,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t),A)) # combining A1, A2, A3 together building the A matrix
###sigmas
###diag
dmgi <- sapply(1:t, function(x) (Z[x,] - muz_std[x]) / nm)
dgi <- t(((Z - muz_std)^2 - diag(sigmaz)) / nm)
###off-diag
if (t > 1) {
ogi <- sapply(1:(t-1), function(x) sapply((x+1):t, function(y)
((Z[x,] - muz_std[x]) * (Z[y,] - muz_std[y]) - sigmaz[x,y]) / nm))
}
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r0-1))/sum(r0-1))
}else sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r0-1))/sum(r0-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.rep.std[[x]]-zbar[,x])*(z.rep.std[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r0-1))/sum(r0-1)))
}
###betas
gi.beta = cbind(1,xhat)*(Y-p) / nm # Construct the residualfrom the logistic model.
# Combine the residuals for the main‐study portion into one matrix
gi.1 = if(t==1){
cbind(dmgi,dgi,gi.beta)
}else cbind(dmgi,dgi,matrix(unlist(ogi),nrow=nm),gi.beta) # adding the off-diagonal
B1 = t(gi.1)%*%gi.1 #the sigmaV_iV_i^T pattern but only for the main study
gi.2 = if(t==1){
dgi.0
}else cbind(dgi.0,matrix(unlist(ogi.0),nrow=nr))# Combine the residuals for the reliability study into one matrix gi.2
B2 = t(gi.2)%*%gi.2 # the sigmaV_iV_i^T pattern but only for the reliable study
m1 = (t)*(t+1)/2 # number of unique covariance parameters (lower‐triangular elements) in main study
m2 = t*(t+1)/2 # number of unique covariance parameters (lower‐triangular elements) in reliable study
# Block [1 : (t+m1), 1 : (t+m1)], from B1 (main study)
# Block (2,2): B2 = reliability part.
# Block [ (t+m1+1):(t+m1+s), (t+m1+1):(t+m1+s) ] from B1 = logistic part.
B = rbind(cbind(B1[1:(t+m1),1:(t+m1)],matrix(0,t+m1,m2),B1[1:(t+m1),(t+m1+1):(t+m1+s)]),
cbind(matrix(0,m2,t+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+m1+1):(t+m1+s),1:(t+m1)],matrix(0,s,m2),B1[(t+m1+1):(t+m1+s),(t+m1+1):(t+m1+s)]))
cov2 = solve(A)%*%B%*%t(solve(A))
var3 = cov2[(t+m+1):(t+m+s),(t+m+1):(t+m+s)] # Extract the submatrix corresponding to the logistic coefficients
tab3 = summary(fit2)$coefficients
tab3[,2] = sqrt(diag(cov2)[(t+m+1):(t+m+s)])
tab3[,1:2] = tab3[,1:2]/c(1,sdz)
CI.low = tab3[,1]-1.96*tab3[,2]
CI.high = tab3[,1]+1.96*tab3[,2]
tab3[,3] = tab3[,1]/tab3[,2] # Z‐value
tab3[,4] = 2*pnorm(tab3[,3],lower.tail=FALSE) # two‐sided p‐value from the normal distribution
tab3 = cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
return(list(
`Sandwich Corrected estimates` = tab3
))
} else{
colnames(xhat) = colnames(z.main.std)
colnames(W.main.std) = colnames(W.main.std)
nm = sum(indicator)
nr = n-nm
q = ncol(W.main.std)
model_matrix = cbind(1, xhat, W.main.std)
colnames(model_matrix) = names(beta.fit2)
p = as.vector(exp(beta.fit2 %*% t(model_matrix))/(1+exp(beta.fit2 %*% t(model_matrix))))
# invert the total covariance matrix of (zmain, Wmain)
# the inv. in the passage, use to correct for the dependencies between Z and W when
# computing sensitivity measures. It appear in the sandwich esitmator to ensure proper variance
sigmaz.inv = solve(sigmaz)
# the transpose of the design matrix for the original measured
# lines up with computing sum(Zi-mu_z) or similar terms in the score equations
Z = t(cbind(z.main.std,W.main.std))
# initialize zero matrices used to mark partial-derivative positions for off-diagonal blocks.
m = matrix(0,nrow = t+q,ncol = t+q)
c = matrix(0,nrow = t,ncol = t+q)
##sigmaX
###diag
# S matrix for diagonal elements of the top-left block sigmaz
ddv.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
cbind(diag(a),matrix(0,nrow = t,ncol = q))
},simplify = FALSE)
# T matrix for diagonal elements
ddm.x = sapply(1:t,function(x){
a = rep(0,t+q)
a[x]<-a[x]+1
diag(a)
},simplify = FALSE)
# compute the final partial derivative derivaritve x_hat w.r.t to sigmaz using the formula
# bi = (S*Inv - v*Inv*T*Inv)*Ci
db.x = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%ddm.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
# Off-Diagonal version:
###off-diag
if(t>1){odv.x = sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y] = c[x,y]+1
c[y,x] = c[y,x]+1
c
},simplify = FALSE),simplify = FALSE)
odm.x = sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y] = m[x,y]+1
m[y,x] = m[y,x]+1
m
},simplify = FALSE),simplify = FALSE)
ob.x = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
}
##sigmaXW (ZW)
##all off-diag
# S for the cross block
odv.xw = sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
c[x,y]<-c[x,y]+1
c
},simplify = FALSE),simplify = FALSE)
# T
odm.xw = sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = FALSE),simplify = FALSE)
# Cross-Covariance, sigma ZW
ob.xw = sapply(1:t,function(x) sapply(1:q,function(y)
t((odv.xw[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.xw[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
##sigmaW
###diag
ddm.w = sapply((t+1):(t+q),function(x){
a = rep(0,t+q)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
# sigmaW
db.w = sapply(1:q,function(x)
t((-v12star%*%sigmaz.inv%*%ddm.w[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
if(q>1){odm.w<-sapply((t+1):(t+q-1),function(x) sapply(min((x+1),t+q):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.w<-sapply(1:(q-1),function(x) sapply(1:(q-x),function(y)
t((-v12star%*%sigmaz.inv%*%odm.w[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = F),simplify = F)
}
##sigma within-person sigma
###diag
db.0<-sapply(1:t,function(x) t((-ddv.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
if(t>1){ob.0 = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((-odv.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)}
m = (t+q)*(t+q+1)/2+(t*(t+1))/2 #number of the covariance estimates
s = t+q+1 #number of beta estimates
#put partial derivatives into a list b, base on t and q
b = if(t>1 & q>1){
list(db.x,db.w,ob.x,ob.xw,ob.w,db.0,ob.0)
}else if(q>1){
list(db.x,db.w,ob.xw,ob.w,db.0)
}else if(t>1){
list(db.x,db.w,ob.x,ob.xw,db.0,ob.0)
}else list(db.x,db.w,ob.xw,db.0)
#Reformat all derivative blocks into m submatrices, each corresponding to a specific covariance parameter.
b = sapply(1:m,function(x) matrix(unlist(b),nrow=nm)[,(t*x-t+1):(t*x)],simplify = FALSE)
# logistic derivative w.r.t. X: p_i(1-p_i) * betaZ
# matches di = p_i(1-p_i)*betaZ in the passage
d = as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)]))/nm
bstar = sapply(b,function(x) -rowSums(x*d))
#Build the cross-derivative block in the big Jacobian A
a = matrix(NA,ncol=m,nrow=s)
a[1,] = colSums(bstar)
a[2:(t+1),] = t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar)))
a[(t+2):(t+q+1),] = t(apply(as.matrix(W.main.std), 2, function(x) colSums(x * bstar))) # A2
w = diag(p * (1 - p))
#∑pi(pi)[1,xhat,Wi][1,xhat,Wi]^T
Astar = t(cbind(1,xhat, W.main.std)) %*% w %*% cbind(1,xhat, W.main.std) /nm #A3
# assemble the entire big matrix A
A = rbind(cbind(diag(rep(-1,m)),matrix(0,nrow=m,ncol=s)),#A1
cbind(a,Astar))
A = rbind(cbind(diag(rep(-1,t+q),nrow=t+q),matrix(0,nrow=t+q,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t+q),A))
##mus
mu = c(muz_std, muw_std)
dmgi = sapply(1:(t+q),function(x) (Z[x,]-mu[x])/nm) # first line of g(u)
###sigmas
###diag
dgi <- t(((Z - mu)^2 - diag(sigmaz)) / nm)
ogi = sapply(1:(t+q-1),function(x) sapply((x+1):(t+q), function(y) # thrid line of g(u)
((Z[x,]-mu[x])*(Z[y,]-mu[y])-sigmaz[x,y])/nm))
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r0-1))/sum(r0-1))
}else sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r0-1))/sum(r0-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.rep.std[[x]]-zbar[,x])*(z.rep.std[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r0-1))/sum(r0-1)))
}
###betas
# last three line of g(u)
gi.beta = cbind(1,xhat,W.main.std)*(Y-p)/nm
gi.1 = cbind(dmgi,dgi,matrix(unlist(ogi),nrow=nm),gi.beta)
B1 = t(gi.1)%*%gi.1
gi.2 = if(t==1){
dgi.0
} else cbind(dgi.0,matrix(unlist(ogi.0),nrow=nr))
B2 = t(gi.2)%*%gi.2
m1 = (t+q)*(t+q+1)/2
m2 = t*(t+1)/2
#Building full B matrix
B = rbind(cbind(B1[1:(t+q+m1),1:(t+q+m1)],matrix(0,t+q+m1,m2),B1[1:(t+q+m1),(t+q+m1+1):(t+q+m1+s)]),
cbind(matrix(0,m2,t+q+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+q+m1+1):(t+q+m1+s),1:(t+q+m1)],matrix(0,s,m2),B1[(t+q+m1+1):(t+q+m1+s),(t+q+m1+1):(t+q+m1+s)]))
# Sandwich formula
cov2 = solve(A)%*%B%*%t(solve(A))
var3 = cov2[(t+q+m+1):(t+q+m+s),(t+q+m+1):(t+q+m+s)]
tab3 = summary(fit2)$coefficients
tab3[,2] = sqrt(diag(cov2)[(t+q+m+1):(t+q+m+s)])
tab3[,1:2] = tab3[,1:2]/c(1,sdz,sdw)
CI.low = tab3[,1]-1.96*tab3[,2]
CI.high = tab3[,1]+1.96*tab3[,2]
tab3[,3] = tab3[,1]/tab3[,2]
tab3[,4] = 2*pnorm(tab3[,3],lower.tail=FALSE)
tab3 = cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
return(list(
`Sandwich Corrected estimates` = tab3
))
}
}
Try the RegCalReliab package in your browser
Any scripts or data that you put into this service are public.
RegCalReliab documentation built on Nov. 6, 2025, 1:18 a.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 sandwich_estimator_ex_log
#' Sandwich Variance Estimator for External Logistic Regression Calibration
#'
#' \code{sandwich_estimator_ex_log()} computes robust (sandwich) standard errors
#' for regression calibration in logistic regression using external reliability data.
#' It propagates uncertainty from estimating the calibration (variance components)
#' into the final coefficient variances, yielding valid confidence intervals and
#' p-values even when measurement error is present.
#'
#' @param xhat Numeric matrix of calibrated exposures (\eqn{n_m \times t}),
#' typically from \code{\link{reg_calibration_ex_log}}.
#' @param z.main.std Numeric matrix of standardized main-study exposures
#' (\eqn{n_m \times t}).
#' @param z.rep.std Named list of standardized replicate measurements from the
#' reliability study. Each list element is a matrix of dimension
#' \eqn{n_r \times r_i}.
#' @param r Integer vector of replicate counts for reliability-study subjects
#' (length \eqn{n_m + n_r}), with 1 for main-study subjects.
#' @param Y Binary (0/1) outcome vector of length \eqn{n_m}.
#' @param indicator Binary vector of length \eqn{n_m + n_r} indicating main-study
#' (1) vs. reliability-study (0) subjects.
#' @param v12star Calibration matrix (\eqn{\Sigma_X \Sigma_Z^{-1}\Sigma_{XZ}})
#' used to transform observed measurements to their calibrated values.
#' @param beta.fit2 Numeric vector of logistic regression coefficients from the
#' calibrated model (\code{fit2}).
#' @param W.main.std Optional numeric matrix of standardized error-free
#' covariates (\eqn{n_m \times q}).
#' @param sigma Within-person variance-covariance matrix of the exposures.
#' @param sigmaz Total variance-covariance matrix of the observed exposures
#' (and covariates, if applicable).
#' @param sigmawithin Estimated within-person covariance matrix (same dimension
#' as \code{sigmaz}).
#' @param sdz Numeric vector of standard deviations of the unstandardized exposures.
#' @param sdw Optional numeric vector of standard deviations of the unstandardized covariates.
#' @param muz Numeric vector of means of the unstandardized exposures.
#' @param muw Optional numeric vector of means of the unstandardized covariates.
#' @param fit2 The \code{glm} object returned by \code{reg_calibration_ex_log()}.
#' @param zbar Matrix of subject-level means of replicate measurements
#' in the reliability study.
#'
#' @return A list with one component:
#' \describe{
#' \item{\code{Sandwich Corrected estimates}}{Matrix of regression calibration
#' estimates with sandwich-corrected standard errors, z-values,
#' p-values, odds ratios, and 95\% confidence intervals on the
#' original scale.}
#' }
#'
#' @details
#' This function implements the three-block estimating equation described in
#' Carroll et al. (2006) and White (1980):
#' \enumerate{
#' \item \strong{Variance components}: estimates total and within-person covariance.
#' \item \strong{Calibration}: computes \eqn{X^\text{hat}} using the
#' estimated variance components.
#' \item \strong{Outcome model}: fits a logistic regression and applies a
#' sandwich (robust) variance formula to combine the uncertainty from
#' all stages.
#' }
#'
#' @examples
#' set.seed(123)
#' # Simulate main-study data
#' z.main <- matrix(rnorm(60), ncol = 1)
#' colnames(z.main) <- "sbp"
#' Y <- rbinom(60, 1, plogis(0.3 * z.main))
#'
#' # Simulate external reliability data (30 subjects, 2 replicates)
#' z.rep <- list(sbp = matrix(rnorm(30 * 2), nrow = 30))
#' zbar <- sapply(z.rep, rowMeans)
#'
#' # Prepare standardization inputs
#' sdz <- apply(z.main, 2, sd)
#' z.main.std <- scale(z.main)
#' z.rep.std <- list(sbp = scale(z.rep$sbp))
#' r <- c(rep(1, 60), rep(2, 30))
#' indicator <- c(rep(1, 60), rep(0, 30))
#'
#' # (In practice, run reg_calibration_ex_log() first to get xhat, fit2, etc.)
#' # Here we show a simplified call using mock objects:
#' # sandwich_estimator_ex_log(xhat, z.main.std, z.rep.std, r, Y, indicator,
#' # v12star, beta.fit2, sigma, sigmaz, sigmawithin, sdz, NULL, muz, NULL, fit2, zbar)
#'
#' @references
#' Carroll RJ, Ruppert D, Stefanski LA, Crainiceanu CM. \emph{Measurement Error
#' in Nonlinear Models: A Modern Perspective}. Chapman & Hall/CRC, 2006.
#' White H. A heteroskedasticity-consistent covariance matrix estimator and a
#' direct test for heteroskedasticity. \emph{Econometrica}. 1980;48(4):817–838.
#'
#' @noRd
sandwich_estimator_ex_log = function(xhat,z.main.std,z.rep.std,r,Y,indicator, v12star,beta.fit2,W.main.std = NULL,sigma,
sigmaz,sigmawithin,sdz,sdw,muz,muw,fit2,zbar){
# -----------------------------------------------
# 0) Basic dimensions
# -----------------------------------------------
nm = nrow(z.main.std) # number of subjects in the main study
nr = sum(indicator == 0)
n = nm+nr # total number of subjects (main + reliability)
t = ncol(z.main.std) # number of variables in z.main.std
r0 = r[indicator==0]
muz_std <- as.numeric(colMeans(z.main.std)) # length t
if (!is.null(W.main.std)) muw_std <- as.numeric(colMeans(W.main.std)) # length q
if(is.null(W.main.std)){
# Compute fitted probabilities p from the corrected (calibrated) logistic model
p = as.vector(exp(beta.fit2 %*% t(cbind(1,xhat)))/(1+exp(beta.fit2 %*% t(cbind(1,xhat)))))
# Compute the inverse of the estimated covariance matrix of Z, sigmaz
sigmaz.inv = solve(sigmaz)
# Create matrix Z as the transpose of the main study data, from nm x t to t x nm
Z = t(cbind(z.main.std))
# Initialize two t x t zero matrices, they will be used later
# to construct off–diagonal derivative matrices.
m = matrix(0,nrow = t,ncol = t)
c = matrix(0,nrow = t,ncol = t)
##sigmaX
###diag
# ddv.x and ddm.x are essentially indicator matrices with a 1 in the (x,x) entry. S and T
ddv.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
ddm.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
# the matrix of partial derivatives (with respect to the (x,x) diagonal element of sigmaz,) of the calibration
# Diagonal Partials wrt sigmaz
db.x = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%ddm.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
# Build off–diagonal derivative matrices
# For each x from 1 to t-1, for each y from x+1 to t, update matrix c
if(t>1){
odv.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y]<-c[x,y]+1
c[y,x]<-c[y,x]+1
c
},simplify = FALSE),simplify = FALSE)
odm.x<-sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = FALSE),simplify = F)
# Compute the off-diagonal derivative contributions
ob.x<-sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
}
# Off‐Diagonal Partials wrt sigmaz
##sigma (the within–group variance)
#Partials wrt Another Variance W
###diag
db.0 = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
# correspond to partial derivatives wrt the diagonal of the other covariance matrix sigmaZW
###off-diag
if(t>1){ob.0 = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)}
# correspond to partial derivatives wrt the off-diagonal of the other covariance matrix sigmaZW
m = (t)*(t+1)/2+(t*(t+1))/2 #number of the covariance estimates
s = t+1 #number of beta estimates (regression parameters)
# If t > 1, include both diagonal and off-diagonal pieces otherwise, only the diagonal pieces
# When t=1 we have only one “error‐prone” covariate—so there are no off‐diagonal elements
# (1×1 covariance matrix has no off‐diagonal entries)
b = if(t>1){
list(db.x,ob.x,db.0,ob.0)
}else list(db.x,db.0)
# Reshape the list b into a list of m matrices
# For each index from 1 to m, extract a matrix with nm rows and t columns from the unlisted b
b = sapply(1:m,function(x) matrix(unlist(b),nrow=nm)[,(t*x-t+1):(t*x)],simplify = FALSE)
d = as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)])) / nm
bstar = sapply(b,function(x) -rowSums(x*d)) # - sum_j b_{ij} * d_{ij}
a = matrix(NA,ncol=m,nrow=s)
a[1,] = colSums(bstar)
a[2:(t+1),] = t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar))) #correspond to A2
w = diag(p * (1 - p))
Astar = t(cbind(1,xhat)) %*% w %*% cbind(1,xhat) /nm # correspond to A3
A = rbind(cbind(diag(rep(-1,m)),
matrix(0,nrow=m,ncol=s)),cbind(a,Astar)) # building up A2 and A3 together
A = rbind(cbind(matrix(diag(rep(-1,t),nrow=t),t),matrix(0,nrow=t,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t),A)) # combining A1, A2, A3 together building the A matrix
###sigmas
###diag
dmgi <- sapply(1:t, function(x) (Z[x,] - muz_std[x]) / nm)
dgi <- t(((Z - muz_std)^2 - diag(sigmaz)) / nm)
###off-diag
if (t > 1) {
ogi <- sapply(1:(t-1), function(x) sapply((x+1):t, function(y)
((Z[x,] - muz_std[x]) * (Z[y,] - muz_std[y]) - sigmaz[x,y]) / nm))
}
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r0-1))/sum(r0-1))
}else sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r0-1))/sum(r0-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.rep.std[[x]]-zbar[,x])*(z.rep.std[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r0-1))/sum(r0-1)))
}
###betas
gi.beta = cbind(1,xhat)*(Y-p) / nm # Construct the residualfrom the logistic model.
# Combine the residuals for the main‐study portion into one matrix
gi.1 = if(t==1){
cbind(dmgi,dgi,gi.beta)
}else cbind(dmgi,dgi,matrix(unlist(ogi),nrow=nm),gi.beta) # adding the off-diagonal
B1 = t(gi.1)%*%gi.1 #the sigmaV_iV_i^T pattern but only for the main study
gi.2 = if(t==1){
dgi.0
}else cbind(dgi.0,matrix(unlist(ogi.0),nrow=nr))# Combine the residuals for the reliability study into one matrix gi.2
B2 = t(gi.2)%*%gi.2 # the sigmaV_iV_i^T pattern but only for the reliable study
m1 = (t)*(t+1)/2 # number of unique covariance parameters (lower‐triangular elements) in main study
m2 = t*(t+1)/2 # number of unique covariance parameters (lower‐triangular elements) in reliable study
# Block [1 : (t+m1), 1 : (t+m1)], from B1 (main study)
# Block (2,2): B2 = reliability part.
# Block [ (t+m1+1):(t+m1+s), (t+m1+1):(t+m1+s) ] from B1 = logistic part.
B = rbind(cbind(B1[1:(t+m1),1:(t+m1)],matrix(0,t+m1,m2),B1[1:(t+m1),(t+m1+1):(t+m1+s)]),
cbind(matrix(0,m2,t+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+m1+1):(t+m1+s),1:(t+m1)],matrix(0,s,m2),B1[(t+m1+1):(t+m1+s),(t+m1+1):(t+m1+s)]))
cov2 = solve(A)%*%B%*%t(solve(A))
var3 = cov2[(t+m+1):(t+m+s),(t+m+1):(t+m+s)] # Extract the submatrix corresponding to the logistic coefficients
tab3 = summary(fit2)$coefficients
tab3[,2] = sqrt(diag(cov2)[(t+m+1):(t+m+s)])
tab3[,1:2] = tab3[,1:2]/c(1,sdz)
CI.low = tab3[,1]-1.96*tab3[,2]
CI.high = tab3[,1]+1.96*tab3[,2]
tab3[,3] = tab3[,1]/tab3[,2] # Z‐value
tab3[,4] = 2*pnorm(tab3[,3],lower.tail=FALSE) # two‐sided p‐value from the normal distribution
tab3 = cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
return(list(
`Sandwich Corrected estimates` = tab3
))
} else{
colnames(xhat) = colnames(z.main.std)
colnames(W.main.std) = colnames(W.main.std)
nm = sum(indicator)
nr = n-nm
q = ncol(W.main.std)
model_matrix = cbind(1, xhat, W.main.std)
colnames(model_matrix) = names(beta.fit2)
p = as.vector(exp(beta.fit2 %*% t(model_matrix))/(1+exp(beta.fit2 %*% t(model_matrix))))
# invert the total covariance matrix of (zmain, Wmain)
# the inv. in the passage, use to correct for the dependencies between Z and W when
# computing sensitivity measures. It appear in the sandwich esitmator to ensure proper variance
sigmaz.inv = solve(sigmaz)
# the transpose of the design matrix for the original measured
# lines up with computing sum(Zi-mu_z) or similar terms in the score equations
Z = t(cbind(z.main.std,W.main.std))
# initialize zero matrices used to mark partial-derivative positions for off-diagonal blocks.
m = matrix(0,nrow = t+q,ncol = t+q)
c = matrix(0,nrow = t,ncol = t+q)
##sigmaX
###diag
# S matrix for diagonal elements of the top-left block sigmaz
ddv.x = sapply(1:t,function(x){
a = rep(0,t)
a[x] = a[x]+1
cbind(diag(a),matrix(0,nrow = t,ncol = q))
},simplify = FALSE)
# T matrix for diagonal elements
ddm.x = sapply(1:t,function(x){
a = rep(0,t+q)
a[x]<-a[x]+1
diag(a)
},simplify = FALSE)
# compute the final partial derivative derivaritve x_hat w.r.t to sigmaz using the formula
# bi = (S*Inv - v*Inv*T*Inv)*Ci
db.x = sapply(1:t,function(x) t((ddv.x[[x]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%ddm.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
# Off-Diagonal version:
###off-diag
if(t>1){odv.x = sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
c[x,y] = c[x,y]+1
c[y,x] = c[y,x]+1
c
},simplify = FALSE),simplify = FALSE)
odm.x = sapply(1:(t-1),function(x) sapply(min((x+1),t):t,function(y){
m[x,y] = m[x,y]+1
m[y,x] = m[y,x]+1
m
},simplify = FALSE),simplify = FALSE)
ob.x = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((odv.x[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
}
##sigmaXW (ZW)
##all off-diag
# S for the cross block
odv.xw = sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
c[x,y]<-c[x,y]+1
c
},simplify = FALSE),simplify = FALSE)
# T
odm.xw = sapply(1:t,function(x) sapply(max((x+1),t+1):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = FALSE),simplify = FALSE)
# Cross-Covariance, sigma ZW
ob.xw = sapply(1:t,function(x) sapply(1:q,function(y)
t((odv.xw[[x]][[y]]%*%sigmaz.inv-
v12star%*%sigmaz.inv%*%odm.xw[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)
##sigmaW
###diag
ddm.w = sapply((t+1):(t+q),function(x){
a = rep(0,t+q)
a[x] = a[x]+1
diag(a)
},simplify = FALSE)
# sigmaW
db.w = sapply(1:q,function(x)
t((-v12star%*%sigmaz.inv%*%ddm.w[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
if(q>1){odm.w<-sapply((t+1):(t+q-1),function(x) sapply(min((x+1),t+q):(t+q),function(y){
m[x,y]<-m[x,y]+1
m[y,x]<-m[y,x]+1
m
},simplify = F),simplify = F)
ob.w<-sapply(1:(q-1),function(x) sapply(1:(q-x),function(y)
t((-v12star%*%sigmaz.inv%*%odm.w[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = F),simplify = F)
}
##sigma within-person sigma
###diag
db.0<-sapply(1:t,function(x) t((-ddv.x[[x]]%*%sigmaz.inv)%*%Z),simplify = FALSE)
###off-diag
if(t>1){ob.0 = sapply(1:(t-1),function(x) sapply(1:(t-x),function(y)
t((-odv.x[[x]][[y]]%*%sigmaz.inv)%*%Z),simplify = FALSE),simplify = FALSE)}
m = (t+q)*(t+q+1)/2+(t*(t+1))/2 #number of the covariance estimates
s = t+q+1 #number of beta estimates
#put partial derivatives into a list b, base on t and q
b = if(t>1 & q>1){
list(db.x,db.w,ob.x,ob.xw,ob.w,db.0,ob.0)
}else if(q>1){
list(db.x,db.w,ob.xw,ob.w,db.0)
}else if(t>1){
list(db.x,db.w,ob.x,ob.xw,db.0,ob.0)
}else list(db.x,db.w,ob.xw,db.0)
#Reformat all derivative blocks into m submatrices, each corresponding to a specific covariance parameter.
b = sapply(1:m,function(x) matrix(unlist(b),nrow=nm)[,(t*x-t+1):(t*x)],simplify = FALSE)
# logistic derivative w.r.t. X: p_i(1-p_i) * betaZ
# matches di = p_i(1-p_i)*betaZ in the passage
d = as.matrix(p*(1-p)%*%t(beta.fit2[2:(t+1)]))/nm
bstar = sapply(b,function(x) -rowSums(x*d))
#Build the cross-derivative block in the big Jacobian A
a = matrix(NA,ncol=m,nrow=s)
a[1,] = colSums(bstar)
a[2:(t+1),] = t(apply(as.matrix(xhat), 2, function(x) colSums(x * bstar)))
a[(t+2):(t+q+1),] = t(apply(as.matrix(W.main.std), 2, function(x) colSums(x * bstar))) # A2
w = diag(p * (1 - p))
#∑pi(pi)[1,xhat,Wi][1,xhat,Wi]^T
Astar = t(cbind(1,xhat, W.main.std)) %*% w %*% cbind(1,xhat, W.main.std) /nm #A3
# assemble the entire big matrix A
A = rbind(cbind(diag(rep(-1,m)),matrix(0,nrow=m,ncol=s)),#A1
cbind(a,Astar))
A = rbind(cbind(diag(rep(-1,t+q),nrow=t+q),matrix(0,nrow=t+q,ncol=m+s)),
cbind(matrix(0,nrow=m+s,ncol=t+q),A))
##mus
mu = c(muz_std, muw_std)
dmgi = sapply(1:(t+q),function(x) (Z[x,]-mu[x])/nm) # first line of g(u)
###sigmas
###diag
dgi <- t(((Z - mu)^2 - diag(sigmaz)) / nm)
ogi = sapply(1:(t+q-1),function(x) sapply((x+1):(t+q), function(y) # thrid line of g(u)
((Z[x,]-mu[x])*(Z[y,]-mu[y])-sigmaz[x,y])/nm))
###within variance
###diag
dgi.0<-if(t==1){
sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-sigma*(r0-1))/sum(r0-1))
}else sapply(1:t,function(y) (rowSums((z.rep.std[[y]]-zbar[,y])^2,na.rm = T)-diag(sigma)[y]*(r0-1))/sum(r0-1))
###off-diag
if(t>1){ogi.0<-sapply(1:(t-1),function(x) sapply((x+1):t, function(y)
(rowSums((z.rep.std[[x]]-zbar[,x])*(z.rep.std[[y]]-zbar[,y]),na.rm = T)-sigma[x,y]*(r0-1))/sum(r0-1)))
}
###betas
# last three line of g(u)
gi.beta = cbind(1,xhat,W.main.std)*(Y-p)/nm
gi.1 = cbind(dmgi,dgi,matrix(unlist(ogi),nrow=nm),gi.beta)
B1 = t(gi.1)%*%gi.1
gi.2 = if(t==1){
dgi.0
} else cbind(dgi.0,matrix(unlist(ogi.0),nrow=nr))
B2 = t(gi.2)%*%gi.2
m1 = (t+q)*(t+q+1)/2
m2 = t*(t+1)/2
#Building full B matrix
B = rbind(cbind(B1[1:(t+q+m1),1:(t+q+m1)],matrix(0,t+q+m1,m2),B1[1:(t+q+m1),(t+q+m1+1):(t+q+m1+s)]),
cbind(matrix(0,m2,t+q+m1),B2,matrix(0,m2,s)),
cbind(B1[(t+q+m1+1):(t+q+m1+s),1:(t+q+m1)],matrix(0,s,m2),B1[(t+q+m1+1):(t+q+m1+s),(t+q+m1+1):(t+q+m1+s)]))
# Sandwich formula
cov2 = solve(A)%*%B%*%t(solve(A))
var3 = cov2[(t+q+m+1):(t+q+m+s),(t+q+m+1):(t+q+m+s)]
tab3 = summary(fit2)$coefficients
tab3[,2] = sqrt(diag(cov2)[(t+q+m+1):(t+q+m+s)])
tab3[,1:2] = tab3[,1:2]/c(1,sdz,sdw)
CI.low = tab3[,1]-1.96*tab3[,2]
CI.high = tab3[,1]+1.96*tab3[,2]
tab3[,3] = tab3[,1]/tab3[,2]
tab3[,4] = 2*pnorm(tab3[,3],lower.tail=FALSE)
tab3 = cbind(tab3,exp(cbind(OR = tab3[, 1],CI.low,CI.high)))
return(list(
`Sandwich Corrected estimates` = tab3
))
}
}
Try the RegCalReliab 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.