Type: Package
Title: Estimation for Binary Emax Models with Missing Responses and Bias Reduction
Version: 0.1.0
Description: Provides estimation utilities for binary Emax dose-response models. Includes Expectation-Maximization based maximum likelihood estimation when the binary response is missing, as well as bias-reduced estimators including Jeffreys-penalized likelihood, Firth-score, and Cox-Snell corrections.The methodology is described in Zhang, Pradhan, and Zhao (2025) <doi:10.1177/09622802251403356> and Zhang, Pradhan, and Zhao (2026) <doi:10.1080/10543406.2026.2627387>.
License: MIT + file LICENSE
Encoding: UTF-8
Depends: R (≥ 4.0.0), clinDR (≥ 2.5.2)
Imports: BB, brglm, boot, formula.tools, MASS, maxLik, numDeriv, stats
Suggests: knitr, rmarkdown, testthat (≥ 3.0.0)
RoxygenNote: 7.3.3
NeedsCompilation: no
Packaged: 2026-03-12 21:13:36 UTC; jiangshanzhang
Author: Jiangshan Zhang [aut, cre], Vivek Pradhan [aut], Yuxi Zhao [aut]
Maintainer: Jiangshan Zhang <jiszhang@ucdavis.edu>
Repository: CRAN
Date/Publication: 2026-03-17 19:10:20 UTC

Augment missing response with observed information

Description

Augment missing response with 0 and 1, and remain all other variables the same.

Usage

Augment_Missing(data)

Arguments

data

Dataset contain missing response indicated as'NA', including response variable as 'y' and dose variable as 'dose'.

Details

DETAILS

Value

A complete dataset with augmentation.

Compute analytical form of Hessian matrix of Binary Emax model

Description

Compute Hessian matrix of Binary Emax model

Usage

Comp_Hess(data, theta, weight)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

weight

A vector of working weight from the EM iterations.

Details

Compute Hessian matrix of Binary Emax model which will be used in finding derivative of Jeffery's prior.

Value

A matrix for Heissian.

Compute analytical form of expected information matrix of Binary Emax model

Description

Compute expected information matrix of Binary Emax model

Usage

Comp_I(data, weight, theta)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

weight

A vector of working weight from the EM iterations.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

Details

Compute expected information matrix of Binary Emax model which will be used in finding derivative of Jeffery's prior.

Value

A matrix of expected information

Maximization function estimation of EM algorithm with defined weight.

Description

Estimate Maximization function with given parameters, weight, and data for the EM Emax model.

Usage

comp_Q(data, theta, alpha, weight, mis_form)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

alpha

Parameters of logisitic missing indicator model.

weight

A vector of working weight from the EM iterations.

mis_form

an object of class "formula": a symbolic description of the model to be fitted.

Details

Calculating the maximization function of EM algorithm for each iteration.

Value

A value of function estimation

Maximization function estimation of bias reduced EM algorithm with defined weight.

Description

Estimate Maximization function with given parameters, weight, and data for the bias reduced EM Emax model.

Usage

comp_Q_firth(data, theta, alpha, weight, fit_mis)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

alpha

Parameters of logisitic missing indicator model.

weight

A vector of working weight from the EM iterations.

fit_mis

an object of class "formula": a symbolic description of the model to be fitted.

Details

Calculating the maximization function of bias reduced EM algorithm for each iteration.

Value

A value of function estimation

Estimation of emax parameters in EM algorithm iteration.

Description

Calls Newton-Raphson optimizer MaxNR, for Emax model.

Usage

comp_theta(data, weight, theta)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

weight

A vector of working weight from the EM iterations.

theta

Initial value of parameters for optimization. The order of the variables is (E0,log(ED50),Emax).

Details

Fits the Emax model with defined working weights using MaxNR.

Value

A vector of parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

Cox–Snell bias-corrected estimator (one-step using clinDR MLE)

Description

Starts from the MLE from clinDR fitEmax and applies a user-supplied Cox–Snell bias correction.

Usage

comp_theta_cox_snell(data = NULL, weight = NULL, theta = NULL)

Arguments

data

A data.frame (or list) with y (0/1) and dose.

weight

Numeric vector of case weights.

theta

Numeric(3) initial guess; ignored if clinDR MLE succeeds.

Details

Requires helpers comp_bias() and Comp_I_score() in your package.

Value

A list with:

par

bias-corrected parameter vector c(e0, emax, led50).

vc

Variance-covariance matrix (generalized inverse of information).

See Also

fitEmax

Examples

 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
 colnames(theta_true)<- c('e_0','emax','led_50')
 theta_true <- as.data.frame(theta_true)
 dose_set <- c(0,7.5,22.5,75,225)
 n=355
 data <-sim_data(theta_true,n,dose_set)
 res <- comp_theta_cox_snell(data=data )

Estimation of emax parameters in Jeffery's prior penalized IL algorithm iteration.

Description

Calls Newton-Raphson optimizer MaxNR, for Emax model.

Usage

comp_theta_firth(data, weight, theta)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

weight

A vector of working weight from the EM iterations.

theta

Initial value of parameters for optimization. The order of the variables is (E0,log(ED50),Emax).

Details

Fits the Emax model with defined working weights using MaxNR with bias reduction.

Value

A vector of parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

Firth-corrected estimating equation solution (score-based)

Description

Solves the Firth-corrected score equations for an Emax-type binary-response model. Requires user-provided score and information helpers.

Usage

comp_theta_firth_score(data = NULL, weight = NULL, theta = NULL)

Arguments

data

A data.frame (or list) with at least y (0/1) and dose.

weight

Numeric vector of case weights, same length as data$y.

theta

Numeric vector (length 3): c(e0, emax, led50) for initialization.

Details

This function depends on internal helpers you must supply in the package: Score_e0(), Score_emax(), Score_led50(), Comp_Fish_inf(), and Comp_I_score().

Value

A list with elements:

par

numeric(3) estimated parameters.

Fisher.inf

Observed/expected information matrix at the solution (from Comp_I_score).

vc

Fisher-based variance/covariance (from Comp_Fish_inf).

score

Score vector evaluated at the solution.

See Also

fitEmax

Examples

 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
 colnames(theta_true)<- c('e_0','emax','led_50')
 theta_true <- as.data.frame(theta_true)
 dose_set <- c(0,7.5,22.5,75,225)
 n=355
 data <-sim_data(theta_true,n,dose_set)
 res <- comp_theta_firth_score(data=data )

Jeffreys-penalized likelihood estimator via Newton–Raphson

Description

Maximizes the Jeffrey's prior-penalized log-likelihood for the binary Emax model using maxLik::maxNR.

Usage

comp_theta_jeffrey(data = NULL, weight = NULL, theta = NULL)

Arguments

data

A data.frame (or list) with y (0/1) and dose.

weight

Numeric vector of case weights.

theta

Numeric(3) initial guess c(e0, emax, led50).

Details

Requires helpers Comp_Hess(), Comp_I(), and Comp_Hess_deriv() that compute the (penalized) Hessian, information, and derivatives of the Hessian.

Value

A list with:

par

Estimated parameter vector.

hessian

Final Hessian returned by optimizer.

vc

Fisher-based variance/covariance.

See Also

fitEmax

Examples

 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
 colnames(theta_true)<- c('e_0','emax','led_50')
 theta_true <- as.data.frame(theta_true)
 dose_set <- c(0,7.5,22.5,75,225)
 n=355
 data <-sim_data(theta_true,n,dose_set)
 res <- comp_theta_jeffrey(data=data )

Calculate the variance covariance matrix of estimated parameters by EmaxEM

Description

This provides the estimated VCOV matrix of parameters using IL method by EmaxEM.

Usage

comp_vcov(em.emax.fit)

Arguments

em.emax.fit

an object for result from EmaxEM

Details

Internal function for variance covariance estimation.

Value

A list of two variance covariance matrices, one for Emax model parameter, one for logistic missing model parameter.

See Also

ginv inv.logit

Calculate the variance covariance matrix of estimated parameters by EmaxEM_firth

Description

This provides the estimated VCOV matrix of parameters using FIL method by EmaxEM_firth.

Usage

comp_vcov_firth(em.emax.fit)

Arguments

em.emax.fit

an object for result from EmaxEM_firth

Details

Internal function for variance covariance estimation.

Value

A list of two variance covariance matrices, one for Emax model parameter, one for logistic missing model parameter.

See Also

ginv inv.logit

Estimation of working weight in EM algorithm iteration.

Description

Internal function for estimating the observation weights for each iteration of EM algorithm.

Usage

comp_weight(data, theta, alpha, mis_form)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

alpha

Parameters of logisitic missing indicator model.

mis_form

an object of class "formula": a symbolic description of the model to be fitted.

Details

Calculating the working weights of EM algorithm for each iteration.

Value

A vector of weights for each observation in the data.

Fitting IL method with Emax model and binary response missing data.

Description

This provides the estimates using IL method.

Usage

fitEmaxEM(
  data = NULL,
  theta_0 = NULL,
  alpha_0 = NULL,
  mis_form = as.formula(mis ~ y + dose + x1 + x2)
)

Arguments

data

Dataset contain missing response indicated as'NA', including response variable as 'y' and dose variable as 'dose'.

theta_0

Initial value of Emax model parameters for optimization, Default: NULL

alpha_0

Initial value of logistic missingness model parameters for optimization., Default: NULL

mis_form

an object of class "formula": a symbolic description of the model to be fitted, Default: as.formula(mis ~ y + dose + x1 + x2)

Value

list of fitted values:

theta

the final fitted parameters of Emax model

alpha

the final fitted parameters of logistic missing model

weight

the final fitted weight for each observation in EM

Q

the value of Q function for maximizationfor each iteration of EM

K

the total number of iterations of EM to converge

vcov_theta

the estimated variance covariance matrix of theta

vcov_alpha

the estimated variance covariance matrix of alpha

See Also

fitEmax

Examples


 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
 colnames(theta_true)<- c('e_0','emax','led_50')
 theta_true <- as.data.frame(theta_true)
 dose_set <- c(0,7.5,22.5,75,225)
 n=355
 alpha_true = c(0.5,1,-0.5,0,0) #mis 15 typical
 data <-sim_data(theta_true,n,dose_set,alpha_true)
 res <- fitEmaxEM(data=data$data,mis_form=as.formula(mis~y+dose) )

Fitting bias reduced IL method with Emax model and binary response missing data.

Description

This provides the estimates using FIL method.

Usage

fitEmaxEM_firth(
  data = NULL,
  theta_0 = NULL,
  alpha_0 = NULL,
  mis_form = as.formula(mis ~ y + dose + x1 + x2)
)

Arguments

data

Dataset contain missing response indicated as'NA', including response variable as 'y' and dose variable as 'dose'.

theta_0

Initial value of Emax model parameters for optimization, Default: NULL

alpha_0

Initial value of logistic missingness model parameters for optimization., Default: NULL

mis_form

an object of class "formula": a symbolic description of the model to be fitted, Default: as.formula(mis ~ y + dose + x1 + x2)

Value

list of fitted values:

theta

the final fitted parameters of Emax model

alpha

the final fitted parameters of logistic missing model

weight

the final fitted weight for each observation in EM

Q

the value of Q function for maximizationfor each iteration of EM

K

the total number of iterations of EM to converge

vcov_theta

the estimated variance covariance matrix of theta

vcov_alpha

the estimated variance covariance matrix of alpha

See Also

fitEmax

Examples


 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
 colnames(theta_true)<- c('e_0','emax','led_50')
 theta_true <- as.data.frame(theta_true)
 dose_set <- c(0,7.5,22.5,75,225)
 n=355
 alpha_true = c(0.5,1,-0.5,0,0) #mis 15 typical
 data <-sim_data(theta_true,n,dose_set,alpha_true)
 res <- fitEmaxEM_firth(data=data$data,mis_form=as.formula(mis~y+dose) )

Log likelihood estimation of binary Emax model

Description

Estimate Log likelihood of given parameters with data for binary Emax model.

Usage

log_Emax_i(data, theta)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

theta

Parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

Details

Internal function of calculating the maximization function of EM algorithm.

Value

A vector of log likelihood of each observation.

Log likelihood estimation of logisitic missing indicator model

Description

Estimate Log likelihood of given parameters with data for logisitic missing indicator model.

Usage

log_missing_i(data, alpha, mis_form)

Arguments

data

A complete dataset without missingness, including response variable as 'y' and dose variable as 'dose'.

alpha

Parameters of logisitic missing indicator model.

mis_form

an object of class "formula": a symbolic description of the model to be fitted.

Details

Internal function of calculating the maximization function of EM algorithm.

Value

A vector of log likelihood of each observation.

Simulate dataset for testing Emaxem and Emaxem_firth

Description

FUNCTION_DESCRIPTION

Usage

sim_data(theta, n, dose_set, alpha = NULL)

Arguments

theta

True parameters of Emax model. The order of the variables is (E0,log(ED50),Emax).

n

number of observations.

dose_set

A vector indicate the dose set for the dose-response relationship.

alpha

True parameters of logisitic missing model.

Details

DETAILS

Value

A list of two datasets. One is with missingness on repsonse, and the other is the full complete data.

Examples

 #EXAMPLE1
 theta_true=matrix(c(qlogis(0.1),qlogis(0.8)-qlogis(0.1),log(7.5)),1,3)
colnames(theta_true)<- c('e_0','emax','led_50')
theta_true <- as.data.frame(theta_true)
dose_set <- c(0,7.5,22.5,75,225)
n=355
alpha_true = c(0.5,1,-0.5,0,0) #mis 15 typical
data <-sim_data(theta_true,n,dose_set,alpha_true)