| Title: | Computes the Variance-Covariance Matrix of Multidimensional Parameters Using M-Estimation |
|---|---|
| Description: | Provides a flexible framework for estimating the variance-covariance matrix of estimated parameters. Estimation relies on unbiased estimating functions to compute the empirical sandwich variance. (i.e., M-estimation in the vein of Tsiatis et al. (2019) <doi:10.1201/9780429192692>. |
| Authors: | François Grolleau |
| Maintainer: | François Grolleau <[email protected]> |
| License: | MIT + file LICENCE |
| Version: | 0.2.1 |
| Built: | 2024-11-07 06:37:24 UTC |
| Source: | CRAN |
Provides a flexible framework for estimating the variance-covariance matrix of a multidimensional parameter. Estimation relies on providing unbiased estimating functions to compute the empirical sandwich variance. (i.e., M-estimation in the vein of Tsiatis et al. (2019) <doi:10.1201/9780429192692>).
get_vcov(data, thetas, M)get_vcov(data, thetas, M)
data |
a dataframe with proper variable (i.e., column) names. |
thetas |
a list of the (properly named) estimated parameters. |
M |
a list of expressions detailing the unbiased estimating functions with the same ordering as |
A list with elements vcov, se, and jacob.
vcov |
pxp matrix of the estimated asymptotic variance-covariance matrix of the estimated parameters in |
se |
length p vector of the estimated asymptotic standard error for the estimated parameters in |
jacob |
a list of lists containing expressions for computing the jacobian matrix. |
François Grolleau <[email protected]>
Stefanski, LA. and Boos DD. (2002)
The Calculus of M-Estimation, The American Statistician,
doi:10.1198/000313002753631330.
Tsiatis, A. A., Davidian, M., Holloway, S. T. and Laber, E. B (2019) Dynamic Treatment Regimes: Statistical Methods for Precision Medicine, CRC Press, doi:10.1201/9780429192692.
#### ## Simulate data #### set.seed(123) n <- 10000 # number of simulated iid observations x_1 <- rnorm(n); x_2 <- rnorm(n) # generate x_1 and x_2 true_thetas <- c(2,3) # generate true parameters X <- model.matrix(~-1+x_1+x_2) # build the design matrix y <- rbinom(n, 1, 1/(1 + exp(-X %*% true_thetas)) ) # generate Y from X and true_thetas dat <- data.frame(x_1=x_1, x_2=x_2, y=y) # build a simulated dataset #### ## Fit a LR model (estimated parameters solve unbiased estimating equations) #### mod <- glm(y~-1 + x_1 + x_2, data=dat, family = "binomial") #### ## Get variance covariance matrix for all parameters solving unbiased estimating equations #### # Put estimated parameters in a list thetas_hat <- list(theta_1=coef(mod)[1], theta_2=coef(mod)[2]) # Build a list of unbiased estimating functions # NB: parameters' names must be consistent with the previous list psi_1 <- expression( ((1/(1+exp( -( theta_1 * x_1 + theta_2 * x_2 ) ))) - y ) * x_1 ) psi_2 <- expression( ((1/(1+exp( -( theta_1 * x_1 + theta_2 * x_2 ) ))) - y ) * x_2 ) est_functions <- list(psi_1, psi_2) ## Pass arguments and run get_vcov res <- get_vcov(data=dat, thetas=thetas_hat, M=est_functions) # Estimted variance covariance matrix is similar to that obtain from glm res$vcov vcov(mod) # So are the standard errors for the estimated parameters res$se summary(mod)$coefficients[,2]#### ## Simulate data #### set.seed(123) n <- 10000 # number of simulated iid observations x_1 <- rnorm(n); x_2 <- rnorm(n) # generate x_1 and x_2 true_thetas <- c(2,3) # generate true parameters X <- model.matrix(~-1+x_1+x_2) # build the design matrix y <- rbinom(n, 1, 1/(1 + exp(-X %*% true_thetas)) ) # generate Y from X and true_thetas dat <- data.frame(x_1=x_1, x_2=x_2, y=y) # build a simulated dataset #### ## Fit a LR model (estimated parameters solve unbiased estimating equations) #### mod <- glm(y~-1 + x_1 + x_2, data=dat, family = "binomial") #### ## Get variance covariance matrix for all parameters solving unbiased estimating equations #### # Put estimated parameters in a list thetas_hat <- list(theta_1=coef(mod)[1], theta_2=coef(mod)[2]) # Build a list of unbiased estimating functions # NB: parameters' names must be consistent with the previous list psi_1 <- expression( ((1/(1+exp( -( theta_1 * x_1 + theta_2 * x_2 ) ))) - y ) * x_1 ) psi_2 <- expression( ((1/(1+exp( -( theta_1 * x_1 + theta_2 * x_2 ) ))) - y ) * x_2 ) est_functions <- list(psi_1, psi_2) ## Pass arguments and run get_vcov res <- get_vcov(data=dat, thetas=thetas_hat, M=est_functions) # Estimted variance covariance matrix is similar to that obtain from glm res$vcov vcov(mod) # So are the standard errors for the estimated parameters res$se summary(mod)$coefficients[,2]