diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 85cd7f2..69cfc6a 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -1,12 +1,17 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +# +# NOTE: This workflow is overkill for most R packages and +# check-standard.yaml is likely a better choice. +# usethis::use_github_action("check-standard") will install it. on: push: - branches: - - master + branches: [main, master] pull_request: - branches: - - master -name: R-CMD-check +name: R-CMD-check.yaml + +permissions: read-all jobs: R-CMD-check: @@ -18,72 +23,40 @@ jobs: fail-fast: false matrix: config: - - {os: macOS-latest, r: 'devel'} - - {os: macOS-latest, r: 'release'} - - {os: macOS-latest, r: 'oldrel'} - - {os: windows-latest, r: 'devel'} + - {os: macos-latest, r: 'release'} + - {os: windows-latest, r: 'release'} - - {os: windows-latest, r: 'oldrel'} - - {os: windows-latest, r: '3.5'} - - {os: ubuntu-16.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} - - {os: ubuntu-16.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} - - {os: ubuntu-16.04, r: 'oldrel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} - - {os: ubuntu-16.04, r: '3.5', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} - - {os: ubuntu-16.04, r: '3.4', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} - - {os: ubuntu-16.04, r: '3.3', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"} + # use 4.0 or 4.1 to check with rtools40's older compiler + - {os: windows-latest, r: 'oldrel-4'} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'oldrel-1'} + - {os: ubuntu-latest, r: 'oldrel-2'} + - {os: ubuntu-latest, r: 'oldrel-3'} + - {os: ubuntu-latest, r: 'oldrel-4'} env: - R_REMOTES_NO_ERRORS_FROM_WARNINGS: true - RSPM: ${{ matrix.config.rspm }} + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - - uses: r-lib/actions/setup-r@master + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + use-public-rspm: true - - uses: r-lib/actions/setup-pandoc@master - - - name: Query dependencies - run: | - install.packages('remotes') - saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) - shell: Rscript {0} - - - name: Cache R packages - if: runner.os != 'Windows' - uses: actions/cache@v1 + - uses: r-lib/actions/setup-r-dependencies@v2 with: - path: ${{ env.R_LIBS_USER }} - key: ${{ runner.os }}-r-new-${{ matrix.config.r }}-1-${{ hashFiles('.github/depends.Rds') }} - restore-keys: ${{ runner.os }}-r-${{ matrix.config.r }}-1- - - - name: Install system dependencies - if: runner.os == 'Linux' - env: - RHUB_PLATFORM: linux-x86_64-ubuntu-gcc - run: | - Rscript -e "remotes::install_github('r-hub/sysreqs')" - sysreqs=$(Rscript -e "cat(sysreqs::sysreq_commands('DESCRIPTION'))") - sudo -s eval "$sysreqs" - - - name: Install dependencies - run: | - remotes::install_deps(dependencies = TRUE) - remotes::install_cran("rcmdcheck") - shell: Rscript {0} - - - name: Check - env: - _R_CHECK_CRAN_INCOMING_: false - run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check") - shell: Rscript {0} + extra-packages: any::rcmdcheck + needs: check - - name: Upload check results - if: failure() - uses: actions/upload-artifact@master + - uses: r-lib/actions/check-r-package@v2 with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check + upload-snapshots: true + build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION deleted file mode 100644 index 0080f95..0000000 --- a/CRAN-SUBMISSION +++ /dev/null @@ -1,3 +0,0 @@ -Version: 1.0.2 -Date: 2024-01-17 06:28:29 UTC -SHA: 70b42698119d523bbacc9294c9238b7014715e7a diff --git a/DESCRIPTION b/DESCRIPTION index 7cd9764..bab771a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: estimatr Type: Package Title: Fast Estimators for Design-Based Inference -Version: 1.0.2 +Version: 1.0.6 Authors@R: c(person("Graeme", "Blair", email = "graeme.blair@gmail.com", role = c("aut", "cre")), person("Jasper", "Cooper", email = "jjc2247@columbia.edu", role = c("aut")), person("Alexander", "Coppock", email = "alex.coppock@yale.edu", role = c("aut")), @@ -9,7 +9,8 @@ Authors@R: c(person("Graeme", "Blair", email = "graeme.blair@gmail.com", role = person("Luke", "Sonnet", email = "luke.sonnet@gmail.com", role = c("aut")), person("Neal", "Fultz", email = "nfultz@gmail.com", role = c("ctb")), person("Lily", "Medina", email = "lilymiru@gmail.com", role = c("ctb")), - person("Russell", "Lenth", email = "russell-lenth@uiowa.edu", role = c("ctb"))) + person("Russell", "Lenth", email = "russell-lenth@uiowa.edu", role = c("ctb")), + person("Molly", "Offer-Westort", email = "mollyow@uchicago.edu", role = c("ctb"))) Description: Fast procedures for small set of commonly-used, design-appropriate estimators with robust standard errors and confidence intervals. Includes estimators for linear regression, instrumental variables regression, difference-in-means, Horvitz-Thompson estimation, and regression improving precision of experimental estimates by interacting treatment with centered pre-treatment covariates introduced by Lin (2013) . URL: https://declaredesign.org/r/estimatr/, https://github.com/DeclareDesign/estimatr BugReports: https://github.com/DeclareDesign/estimatr/issues diff --git a/NEWS.md b/NEWS.md index c91ef9c..77de895 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,17 @@ +# estimatr 1.0.6 + +* Allows for prediction with lm_lin() when treatment is a factor and/or multi-valued. +* Adds saved treatment_levels to the returned lm_lin model object. +* Stops prediction for lm_lin if the treatment values in new data are not a subset of treatment_levels. +* Standardizes model fit for lm_lin() models with no intercept. +* Adds tests to ensure identical predictions from lm_lin() models where treatment is either numeric or factorial, and fit with/without an intercept. +* Adds relevant examples to predict and lm_robust and lm_lin documentation. +* Adds Molly Offer-Westort as a contributor. + +# estimatr 1.0.4 + +* Test suite changes for M1 mac stay current on CRAN. + # estimatr 1.0.2 * Minor documentation changes to stay current on CRAN. diff --git a/R/S3_predict.R b/R/S3_predict.R index 2e2bbfb..5aa4740 100644 --- a/R/S3_predict.R +++ b/R/S3_predict.R @@ -61,6 +61,27 @@ #' new_dat$w <- runif(n) #' predict(lm_out, newdata = new_dat, weights = w, interval = "prediction") #' +#' # Works for 'lm_lin' models as well +#' dat$z <- sample(1:3, size = nrow(dat), replace = TRUE) +#' lmlin_out1 <- lm_lin(y ~ z, covariates = ~ x, data = dat) +#' predict(lmlin_out1, newdata = dat, interval = "prediction") +#' +#' # Predictions from Lin models are equivalent with and without an intercept +#' # and for multi-level treatments entered as numeric or factor variables +#' lmlin_out2 <- lm_lin(y ~ z - 1, covariates = ~ x, data = dat) +#' lmlin_out3 <- lm_lin(y ~ factor(z), covariates = ~ x, data = dat) +#' lmlin_out4 <- lm_lin(y ~ factor(z) - 1, covariates = ~ x, data = dat) +#' +#' predict(lmlin_out2, newdata = dat, interval = "prediction") +#' predict(lmlin_out3, newdata = dat, interval = "prediction") +#' predict(lmlin_out4, newdata = dat, interval = "prediction") +#' +#' # In Lin models, predict will stop with an error message if new +#' # treatment levels are supplied in the new data +#' new_dat$z <- sample(0:3, size = nrow(new_dat), replace = TRUE) +#' # predict(lmlin_out, newdata = new_dat) +#' +#' #' @export predict.lm_robust <- function(object, newdata, @@ -74,30 +95,6 @@ predict.lm_robust <- function(object, X <- get_X(object, newdata, na.action) - # lm_lin scaling - if (!is.null(object$scaled_center)) { - demeaned_covars <- - scale( - X[ - , - names(object$scaled_center), - drop = FALSE - ], - center = object$scaled_center, - scale = FALSE - ) - - # Interacted with treatment - treat_name <- attr(object$terms, "term.labels")[1] - interacted_covars <- X[, treat_name] * demeaned_covars - - X <- cbind( - X[, attr(X, "assign") <= 1, drop = FALSE], - demeaned_covars, - interacted_covars - ) - } - # Get coefs coefs <- as.matrix(coef(object)) @@ -224,9 +221,64 @@ get_X <- function(object, newdata, na.action) { X <- model.matrix(rhs_terms, mf, contrasts.arg = object$contrasts) + # lm_lin scaling (moved down from predict.lm_robust) + if (!is.null(object$scaled_center)) { + # Covariates + demeaned_covars <- + scale( + X[ + , + names(object$scaled_center), + drop = FALSE + ], + center = object$scaled_center, + scale = FALSE + ) + + # Handle treatment variable reconstruction + treat_name <- attr(object$terms, "term.labels")[1] + treatment <- mf[, treat_name] + vals <- sort(unique(treatment)) + old_vals <- object$treatment_levels + + # Ensure treatment levels in newdata are subset of those for model fit + if (!all(as.character(vals) %in% as.character(old_vals))) { + stop( + "Levels of treatment variable in `newdata` must be a subset of those ", + "in the model fit." + ) + } + treatment <- model.matrix(~ factor(treatment, levels = old_vals) - 1) + + colnames(treatment) <- paste0(treat_name, "_", old_vals) + # Drop out first group if there is an intercept + if (attr(rhs_terms, "intercept") == 1) treatment <- treatment[, -1, drop = FALSE] + + # Interactions matching original fitting logic + n_treat_cols <- ncol(treatment) + n_covars <- ncol(demeaned_covars) + + interaction_matrix <- matrix(0, nrow = nrow(X), ncol = n_covars * n_treat_cols) + + for (i in 1:n_covars) { + cols <- (i - 1) * n_treat_cols + (1:n_treat_cols) + interaction_matrix[, cols] <- treatment * demeaned_covars[, i] + } + + X <- cbind( + if (attr(rhs_terms, "intercept") == 1) { + matrix(1, nrow = nrow(X), ncol = 1, dimnames = list(NULL, "(Intercept)")) + }, + treatment, + if (attr(rhs_terms, "intercept") == 1 || ncol(treatment) == 1) demeaned_covars, + interaction_matrix + ) + } + return(X) } + add_fes <- function(preds, object, newdata) { # Add factors! diff --git a/R/estimatr_lm_lin.R b/R/estimatr_lm_lin.R index 37a0eab..23ce146 100644 --- a/R/estimatr_lm_lin.R +++ b/R/estimatr_lm_lin.R @@ -83,7 +83,7 @@ #' \item{weighted}{whether or not weights were applied} #' \item{call}{the original function call} #' \item{fitted.values}{the matrix of predicted means} -#' We also return \code{terms} and \code{contrasts}, used by \code{predict}, +#' We also return \code{terms}, \code{contrasts}, and \code{treatment_levels}, used by \code{predict}, #' and \code{scaled_center} (the means of each of the covariates used for centering them). #' #' @seealso \code{\link{lm_robust}} @@ -127,9 +127,12 @@ #' #' lm_lin(y ~ z_clust, covariates = ~ x, data = dat, clusters = clusterID) #' -#' # Works with multi-valued treatments +#' # Works with multi-valued treatments, whether treatment is specified as a +#' # factor or not #' dat$z_multi <- sample(1:3, size = nrow(dat), replace = TRUE) +#' #' lm_lin(y ~ z_multi, covariates = ~ x, data = dat) +#' lm_lin(y ~ factor(z_multi), covariates = ~ x, data = dat) #' #' # Stratified estimator with blocks #' dat$blockID <- rep(1:5, each = 8) @@ -137,11 +140,23 @@ #' #' lm_lin(y ~ z_block, ~ factor(blockID), data = dat) #' +#' # Fitting the model without an intercept provides estimates of mean outcomes +#' # under each respective treatment condition +#' lm_lin(y ~ z_multi - 1, covariates = ~ x, data = dat) +#' +#' # Predictions are the same in equivalent models with and without an intercept +#' lmlin_out3 <- lm_lin(y ~ z_multi - 1, covariates = ~ x, data = dat) +#' lmlin_out4 <- lm_lin(y ~ z_multi, covariates = ~ x, data = dat) +#' +#' predict(lmlin_out3, newdata = dat, se.fit = TRUE, interval = "confidence") +#' predict(lmlin_out4, newdata = dat, se.fit = TRUE, interval = "confidence") +#' #' \dontrun{ #' # Can also use 'margins' package if you have it installed to get #' # marginal effects #' library(margins) -#' lmlout <- lm_lin(y ~ z_block, ~ x, data = dat) +#' # Instruct 'margins' to treat z as a factor +#' lmlout <- lm_lin(y ~ factor(z_block), ~ x, data = dat) #' summary(margins(lmlout)) #' #' # Can output results using 'texreg' @@ -230,7 +245,7 @@ lm_lin <- function(formula, design_mat_treatment <- colnames(design_matrix)[treat_col] # Check case where treatment is not factor and is not binary - if (any(!(treatment %in% c(0, 1)))) { + if (any(!(treatment %in% c(0, 1))) | (!has_intercept&ncol(treatment) ==1) ) { # create dummies for non-factor treatment variable # Drop out first group if there is an intercept @@ -313,20 +328,10 @@ lm_lin <- function(formula, interacted_covars ) } else { - # If no intercept, but treatment is only one column, - # need to add base terms for covariates - if (n_treat_cols == 1) { - X <- cbind( - treatment, - demeaned_covars, - interacted_covars - ) - } else { - X <- cbind( - treatment, - interacted_covars - ) - } + X <- cbind( + treatment, + interacted_covars + ) } # ---------- @@ -360,6 +365,12 @@ lm_lin <- function(formula, return_list[["scaled_center"]] <- center setNames(return_list[["scaled_center"]], original_covar_names) + # Store unique treatment values + if(attr(terms(model_data), "dataClasses")[attr(terms(model_data),"term.labels")[1]] == "factor"){ + return_list[["treatment_levels"]] <- model_data$xlevels[[1]] + } else { + return_list[["treatment_levels"]] <- sort(unique(design_matrix[, design_mat_treatment])) + } return_list[["call"]] <- match.call() diff --git a/cran-comments.md b/cran-comments.md index ca83863..0214e44 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,6 +1,6 @@ ## Submission -Small changes to update based on Kurt's 8/19/23 email and also address C++ warning in R CMD CHECK. +Fixes M1 Mac errors. There are no changes to worse in reverse depends. diff --git a/estimatr.Rproj b/estimatr.Rproj index bae6b6c..1c4d969 100644 --- a/estimatr.Rproj +++ b/estimatr.Rproj @@ -1,5 +1,5 @@ Version: 1.0 -ProjectId: 533c90b7-36c9-4d31-ae3d-2adc22079fe1 +ProjectId: 2690710e-30ff-4c33-a75e-7b4e02172a9d RestoreWorkspace: Default SaveWorkspace: Default diff --git a/man/estimatr.Rd b/man/estimatr.Rd index a7aedcd..7c86faa 100644 --- a/man/estimatr.Rd +++ b/man/estimatr.Rd @@ -33,6 +33,7 @@ Other contributors: \item Neal Fultz \email{nfultz@gmail.com} [contributor] \item Lily Medina \email{lilymiru@gmail.com} [contributor] \item Russell Lenth \email{russell-lenth@uiowa.edu} [contributor] + \item Molly Offer-Westort \email{mollyow@uchicago.edu} [contributor] } } diff --git a/man/lm_lin.Rd b/man/lm_lin.Rd index 14e8c14..8ba0149 100644 --- a/man/lm_lin.Rd +++ b/man/lm_lin.Rd @@ -95,7 +95,7 @@ following components: \item{weighted}{whether or not weights were applied} \item{call}{the original function call} \item{fitted.values}{the matrix of predicted means} -We also return \code{terms} and \code{contrasts}, used by \code{predict}, +We also return \code{terms}, \code{contrasts}, and \code{treatment_levels}, used by \code{predict}, and \code{scaled_center} (the means of each of the covariates used for centering them). } \description{ @@ -149,9 +149,12 @@ dat$z_clust <- cluster_ra(clusters = dat$clusterID) lm_lin(y ~ z_clust, covariates = ~ x, data = dat, clusters = clusterID) -# Works with multi-valued treatments +# Works with multi-valued treatments, whether treatment is specified as a +# factor or not dat$z_multi <- sample(1:3, size = nrow(dat), replace = TRUE) + lm_lin(y ~ z_multi, covariates = ~ x, data = dat) +lm_lin(y ~ factor(z_multi), covariates = ~ x, data = dat) # Stratified estimator with blocks dat$blockID <- rep(1:5, each = 8) @@ -159,11 +162,23 @@ dat$z_block <- block_ra(blocks = dat$blockID) lm_lin(y ~ z_block, ~ factor(blockID), data = dat) +# Fitting the model without an intercept provides estimates of mean outcomes +# under each respective treatment condition +lm_lin(y ~ z_multi - 1, covariates = ~ x, data = dat) + +# Predictions are the same in equivalent models with and without an intercept +lmlin_out3 <- lm_lin(y ~ z_multi - 1, covariates = ~ x, data = dat) +lmlin_out4 <- lm_lin(y ~ z_multi, covariates = ~ x, data = dat) + +predict(lmlin_out3, newdata = dat, se.fit = TRUE, interval = "confidence") +predict(lmlin_out4, newdata = dat, se.fit = TRUE, interval = "confidence") + \dontrun{ # Can also use 'margins' package if you have it installed to get # marginal effects library(margins) - lmlout <- lm_lin(y ~ z_block, ~ x, data = dat) + # Instruct 'margins' to treat z as a factor + lmlout <- lm_lin(y ~ factor(z_block), ~ x, data = dat) summary(margins(lmlout)) # Can output results using 'texreg' diff --git a/man/predict.lm_robust.Rd b/man/predict.lm_robust.Rd index 27792d4..43257bc 100644 --- a/man/predict.lm_robust.Rd +++ b/man/predict.lm_robust.Rd @@ -90,4 +90,25 @@ predict(lm_out, newdata = new_dat) new_dat$w <- runif(n) predict(lm_out, newdata = new_dat, weights = w, interval = "prediction") +# Works for 'lm_lin' models as well +dat$z <- sample(1:3, size = nrow(dat), replace = TRUE) +lmlin_out1 <- lm_lin(y ~ z, covariates = ~ x, data = dat) +predict(lmlin_out1, newdata = dat, interval = "prediction") + +# Predictions from Lin models are equivalent with and without an intercept +# and for multi-level treatments entered as numeric or factor variables +lmlin_out2 <- lm_lin(y ~ z - 1, covariates = ~ x, data = dat) +lmlin_out3 <- lm_lin(y ~ factor(z), covariates = ~ x, data = dat) +lmlin_out4 <- lm_lin(y ~ factor(z) - 1, covariates = ~ x, data = dat) + +predict(lmlin_out2, newdata = dat, interval = "prediction") +predict(lmlin_out3, newdata = dat, interval = "prediction") +predict(lmlin_out4, newdata = dat, interval = "prediction") + +# In Lin models, predict will stop with an error message if new +# treatment levels are supplied in the new data +new_dat$z <- sample(0:3, size = nrow(new_dat), replace = TRUE) +# predict(lmlin_out, newdata = new_dat) + + } diff --git a/tests/testthat/helper-lm-robust-se.R b/tests/testthat/helper-lm-robust-se.R index 4eb98d0..53e6307 100644 --- a/tests/testthat/helper-lm-robust-se.R +++ b/tests/testthat/helper-lm-robust-se.R @@ -164,3 +164,4 @@ BMlmSE <- function(model, clustervar=NULL, ell=NULL, IK=TRUE) { se = se, se.Stata = se.Stata ) } + diff --git a/tests/testthat/test-lm-cluster.R b/tests/testthat/test-lm-cluster.R index a860409..164548c 100644 --- a/tests/testthat/test-lm-cluster.R +++ b/tests/testthat/test-lm-cluster.R @@ -122,7 +122,7 @@ test_that("lm cluster se", { for (se_type in cr_se_types) { lmr_rd <- lm_robust(Y ~ X + Z + X2, data = dat, clusters = J, se_type = se_type) lmr_full <- lm_robust(Y ~ X + Z, data = dat, clusters = J, se_type = se_type) - expect_identical( + expect_equal( tidy(lmr_rd)[1:3, ], tidy(lmr_full) ) diff --git a/tests/testthat/test-lm-lin.R b/tests/testthat/test-lm-lin.R index 555f6ed..6a97e33 100644 --- a/tests/testthat/test-lm-lin.R +++ b/tests/testthat/test-lm-lin.R @@ -308,10 +308,19 @@ test_that("Test LM Lin", { lmlo <- lm_lin(Y ~ z + 0, ~X1, data = dat) expect_equal( lmlo$term, - c("z", "X1_c", "z:X1_c") + c("z0", "z1", "z0:X1_c", "z1:X1_c") + ) + + # works with a multi-value treatment with no intercept + dat$z <- factor(rbinom(nrow(dat), 2, 0.5)) + lmlo <- lm_lin(Y ~ z + 0, ~X1, data = dat) + expect_equal( + lmlo$term, + c("z0", "z1", "z2", "z0:X1_c", "z1:X1_c", "z2:X1_c") ) # behaves correctly with "odd" covariates (see issue #283) + dat$z <- rbinom(nrow(dat), 1, 0.5) lmlo <- lm_lin(Y ~ z, ~ is.na(X1), data = dat) expect_equal( lmlo$term, diff --git a/tests/testthat/test-lm-robust_margins.R b/tests/testthat/test-lm-robust_margins.R index 433e6c3..76140dd 100644 --- a/tests/testthat/test-lm-robust_margins.R +++ b/tests/testthat/test-lm-robust_margins.R @@ -97,17 +97,19 @@ test_that("lm robust + cluster can work with margins", { test_that("lm lin can work with margins", { skip_if_not_installed("margins") data("alo_star_men") - lml <- lm_lin(GPA_year1 ~ ssp, ~ gpa0, data = alo_star_men, se_type = "classical") + # instruct margins to treat treatment as a factor + lml <- lm_lin(GPA_year1 ~ factor(ssp), ~ gpa0, data = alo_star_men, se_type = "classical") alo_star_men$gpa0_tilde <- alo_star_men$gpa0 - mean(alo_star_men$gpa0) - lmo <- lm(GPA_year1 ~ ssp * gpa0_tilde, data = alo_star_men) + lmo <- lm(GPA_year1 ~ factor(ssp) * gpa0_tilde, data = alo_star_men) lml_sum <- margins:::summary.margins(margins::margins(lml, vce = "delta")) lmo_sum <- margins:::summary.margins(margins::margins(lmo, vce = "delta")) expect_equal( - round(lml_sum[, 4], 5), - round(lmo_sum[, 4], 5) + lml_sum[, 4], + lmo_sum[, 4], + tolerance = 0.000001 ) }) diff --git a/tests/testthat/test-modelsummary.R b/tests/testthat/test-modelsummary.R index 451f457..ada0553 100644 --- a/tests/testthat/test-modelsummary.R +++ b/tests/testthat/test-modelsummary.R @@ -3,10 +3,12 @@ context("S3 - modelsummary works") test_that("modelsummary works with glance", { - # skip_if_not_installed("modelsummary") - skip("modelsummary non reproducible errors") + skip_if_not_installed("modelsummary") library(modelsummary) + + set.seed(5) + model1 <- lm_robust(mpg ~ am, mtcars) model2 <- lm_robust(mpg ~ am, mtcars, clusters = cyl) model3 <- lm_lin(mpg ~ am, ~ cyl, mtcars) @@ -14,9 +16,9 @@ test_that("modelsummary works with glance", { mso <- modelsummary(list(model1, model2, model3), output = "data.frame") expect_equal(colnames(mso), c("part", "term", "statistic", - "Model 1", "Model 2", "Model 3")) + "(1)", "(2)", "(3)")) - expect_equal(nrow(mso), 12L) + expect_equal(nrow(mso), 15L) expect_equal(ncol(mso), 6L) @@ -28,7 +30,7 @@ test_that("modelsummary works with glance", { gof_omit = c("N|[sS]tatistic|p.value|p{1}"), output = "data.frame") - expect_equal(nrow(mso), 6) + expect_equal(nrow(mso), 10) expect_equal(ncol(mso), 5) diff --git a/tests/testthat/test-s3-methods.R b/tests/testthat/test-s3-methods.R index 4ae25ab..58bfd49 100644 --- a/tests/testthat/test-s3-methods.R +++ b/tests/testthat/test-s3-methods.R @@ -644,6 +644,59 @@ test_that("predict works", { predict(lm_int_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] ) + # With and without factorial treatment, intercept, predictions should be the same + lml_out <- lm_lin(y ~ x, covariates = ~ z + I(z^2) +cl, data = dat, se_type = "classical") + lml_0_out <- lm_lin(y ~ x + 0, covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + lml_f_out <- lm_lin(y ~ factor(x), covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + lm_f0_out <- lm_lin(y ~ factor(x) + 0, covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + + expect_equivalent( + predict(lml_0_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + ) + + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lm_f0_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + + # For multi-level treatments, with and without factorial treatment, intercept, predictions should be the same + dat$x_multi <- rep(0:2, times = 7)[1:20] + new_dat$x_multi <- rep(0:2, times = 7) + + lml_out <- lm_lin(y ~ x_multi, covariates = ~ z + I(z^2) +cl, data = dat, se_type = "classical") + lml_0_out <- lm_lin(y ~ x_multi + 0, covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + lml_f_out <- lm_lin(y ~ factor(x_multi), covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + lm_f0_out <- lm_lin(y ~ factor(x_multi) + 0, covariates = ~ z + I(z^2) + cl, data = dat, se_type = "classical") + + expect_equivalent( + predict(lml_0_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + ) + + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lml_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + + expect_equivalent( + predict(lml_f_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)], + predict(lm_f0_out, new_dat, se.fit = TRUE, interval = "confidence")[c(1, 2)] + ) + # working with rank deficient X head(dat) dat$z2 <- dat$z