diff --git a/.Rbuildignore b/.Rbuildignore
index b7287de76..dffa79a3b 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -2,3 +2,4 @@
 ^README\.Rmd$
 ^README\.html$
 ^\.github$
+^\.here$
diff --git a/.here b/.here
new file mode 100644
index 000000000..e69de29bb
diff --git a/DESCRIPTION b/DESCRIPTION
index 806c383ac..98eaa69df 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -27,6 +27,7 @@ Suggests:
     fixest,
     ggeffects,
     haven,
+    here,
     ivreg,
     lme4,
     MASS,
diff --git a/tests/testthat/stata/clean.R b/tests/testthat/stata/clean.R
new file mode 100755
index 000000000..df0ce49c7
--- /dev/null
+++ b/tests/testthat/stata/clean.R
@@ -0,0 +1,50 @@
+library(haven)
+library(ivreg)
+library(betareg)
+library(MASS)
+library(data.table)
+library(testthat)
+library(readxl)
+
+results <- list()
+
+#  stats::glm
+tmp <- read.table("results/stats_glm_01.txt", sep = "\t", skip = 4) |>
+    setNames(c("term", "dydxstata", "std.errorstata")) |>
+    head(2)
+results[["stats_glm_01"]] <- tmp
+
+#  stats::lm
+tmp <- read.table("results/stats_lm_01.txt", sep = "\t", skip = 4) |>
+    setNames(c("term", "dydxstata", "std.errorstata")) |>
+    head(2)
+results[["stats_lm_01"]] <- tmp
+
+#  MASS::polr 
+tmp <- read.table("results/MASS_polr_01.txt", sep = "\t", skip = 4) |>
+    setNames(c("group", "dydxstata_x1", "std.errorstata_x1", "dydxstata_x2", "std.errorstata_x2")) |>
+    head(3) |>
+    transform(group = 2:4) |>
+    data.table()
+cols <- names(tmp)
+tmp[, (cols) := lapply(.SD, as.numeric)]
+tmp <- melt(tmp, id.vars = "group")
+tmp[, c("statistic", "term") := tstrsplit(variable, "_")]
+tmp[, variable := NULL]
+tmp <- dcast(tmp, group + term ~ statistic, value.var = "value")
+results[["MASS_polr_01"]] <- tmp
+
+#  ivreg::ivreg
+tmp <- read.table("results/ivreg_ivreg_01.txt", sep = "\t", skip = 4) |>
+    setNames(c("term", "dydxstata", "std.errorstata")) |>
+    head(2)
+results[["ivreg_ivreg_01"]] <- tmp
+
+#  betareg::betareg 
+tmp <- read.table("results/betareg_betareg_01.txt", sep = "\t", skip = 4) |>
+    setNames(c("term", "dydxstata", "std.errorstata")) |>
+    head(1)
+results[["betareg_betareg_01"]] <- tmp
+
+# save
+saveRDS(results, file = "stata.rds")
diff --git a/tests/testthat/stata/data/MASS_polr_01.dta b/tests/testthat/stata/data/MASS_polr_01.dta
new file mode 100755
index 000000000..be1ddefff
Binary files /dev/null and b/tests/testthat/stata/data/MASS_polr_01.dta differ
diff --git a/tests/testthat/stata/data/betareg_betareg_01.dta b/tests/testthat/stata/data/betareg_betareg_01.dta
new file mode 100755
index 000000000..67277b534
Binary files /dev/null and b/tests/testthat/stata/data/betareg_betareg_01.dta differ
diff --git a/tests/testthat/stata/data/fixest_fixest_01.dta b/tests/testthat/stata/data/fixest_fixest_01.dta
new file mode 100755
index 000000000..8eec33397
Binary files /dev/null and b/tests/testthat/stata/data/fixest_fixest_01.dta differ
diff --git a/tests/testthat/stata/data/ivreg_ivreg_01.dta b/tests/testthat/stata/data/ivreg_ivreg_01.dta
new file mode 100755
index 000000000..46d2a378f
Binary files /dev/null and b/tests/testthat/stata/data/ivreg_ivreg_01.dta differ
diff --git a/tests/testthat/stata/data/stats_glm_01.dta b/tests/testthat/stata/data/stats_glm_01.dta
new file mode 100755
index 000000000..83c27c854
Binary files /dev/null and b/tests/testthat/stata/data/stats_glm_01.dta differ
diff --git a/tests/testthat/stata/data/stats_lm_01.dta b/tests/testthat/stata/data/stats_lm_01.dta
new file mode 100755
index 000000000..081ad3e51
Binary files /dev/null and b/tests/testthat/stata/data/stats_lm_01.dta differ
diff --git a/tests/testthat/stata/estimate.do b/tests/testthat/stata/estimate.do
new file mode 100755
index 000000000..3caa228b6
--- /dev/null
+++ b/tests/testthat/stata/estimate.do
@@ -0,0 +1,41 @@
+cd ~/Dropbox/marginsxp_crosscheck
+
+
+* stats::glm
+clear
+use data/stats_glm_01.dta
+quiet logit y c.x1##c.x2
+quiet margins, dydx(x1 x2) post
+outreg2 using "results/stats_glm_01.xls", dec(10) excel replace noaster sideway noparen stats(coef se)
+
+
+* stats::lm
+clear
+use data/stats_lm_01.dta
+quiet reg y c.x1##c.x2
+quiet margins, dydx(x1 x2) post
+outreg2 using "results/stats_lm_01.xls", dec(10) excel replace noaster sideway noparen stats(coef se)
+
+
+* betareg::betareg
+clear
+use data/betareg_betareg_01.dta
+quiet betareg yield i.batch temp
+quiet margins, dydx(temp) post
+outreg2 using "results/betareg_betareg_01.xls", dec(10) excel replace noaster sideway noparen stats(coef se)
+
+
+* MASS::polr
+clear
+use data/MASS_polr_01.dta
+quiet ologit y x1 x2
+quiet margins, dydx(x1 x2) post
+outreg2 using "results/MASS_polr_01.xls", dec(10) excel replace noaster sideway noparen stats(coef se)
+
+
+* ivreg::ivreg
+clear
+use data/ivreg_ivreg_01.dta
+quiet ivregress 2sls Q D (P = D F A)
+quiet margins, dydx(P D) post
+outreg2 using "results/ivreg_ivreg_01.xls", dec(10) excel replace noaster sideway noparen stats(coef se)
diff --git a/tests/testthat/stata/generate.R b/tests/testthat/stata/generate.R
new file mode 100755
index 000000000..008721b30
--- /dev/null
+++ b/tests/testthat/stata/generate.R
@@ -0,0 +1,54 @@
+library(haven)
+library(ivreg)
+library(betareg)
+
+################
+#  stats::glm  #
+################
+set.seed(1024)
+N <- 100
+dat <- data.frame(x1 = rnorm(N), x2 = rnorm(N))
+dat$y <- plogis(dat$x1 + dat$x2 + dat$x1 * dat$x2)
+dat$y <- rbinom(N, 1, dat$y)
+haven::write_dta(dat, path = "data/stats_glm_01.dta")
+
+
+###############
+#  stats::lm  #
+###############
+set.seed(1024)
+N <- 100
+dat <- data.frame(x1 = rnorm(N), x2 = rnorm(N))
+dat$y <- dat$x1 + dat$x2 + dat$x1 * dat$x2 + rnorm(N)
+haven::write_dta(dat, path = "data/stats_lm_01.dta")
+
+
+################
+#  MASS::polr  #
+################
+set.seed(1024)
+N <- 1000
+dat <- data.frame(x1 = rnorm(N), x2 = rnorm(N))
+dat$y <- dat$x1 + dat$x2 + dat$x1 * dat$x2 + rnorm(N)
+dat$y <- cut(dat$y, breaks = 4)
+dat$y <- factor(as.numeric(dat$y))
+haven::write_dta(dat, path = "data/MASS_polr_01.dta")
+
+
+##################
+#  ivreg::ivreg  #
+##################
+haven::write_dta(ivreg::Kmenta, path = "data/ivreg_ivreg_01.dta")
+
+
+######################
+#  betareg::betareg  #
+######################
+data("GasolineYield", package = "betareg")
+haven::write_dta(GasolineYield, path = "data/betareg_betareg_01.dta")
+
+
+####################
+#  fixest::fixest  #
+####################
+haven::write_dta(mtcars, path = "data/fixest_fixest_01.dta")
diff --git a/tests/testthat/stata/stata.rds b/tests/testthat/stata/stata.rds
new file mode 100755
index 000000000..ed86d6f4b
Binary files /dev/null and b/tests/testthat/stata/stata.rds differ
diff --git a/tests/testthat/test-dydx_MASS.R b/tests/testthat/test-dydx_MASS.R
index 7191103d1..bf66b8a67 100644
--- a/tests/testthat/test-dydx_MASS.R
+++ b/tests/testthat/test-dydx_MASS.R
@@ -1,19 +1,33 @@
 skip_if_not_installed("MASS")
 
 library("margins")
+library("haven")
+library("data.table")
 
-test_that("polr: test against margins", {
+test_that("polr vs. margins", {
     skip("no idea why this fails")
     tmp <- data.frame(mtcars)
     tmp$carb <- as.factor(tmp$carb)
     mod <- MASS::polr(carb ~ hp + am + mpg, data = tmp) 
     res <- marginsxp(mod, variance = NULL)
     mar <- margins(mod)
-
     expect_s3_class(res, "data.frame")
     expect_equal(dim(res), c(480, 8))
-
     # TODO: not supported yet
     expect_error(marginsxp(mod, variance = NULL), regexp = "group_name")
     expect_error(marginsxp(mod, group_names = "1"), regexp = "not yet supported")
 })
+
+test_that("polr vs. Stata", {
+    stata <- readRDS(test_path("stata/stata.rds"))[["MASS_polr_01"]]
+    dat <- read_dta(test_path("stata/data/MASS_polr_01.dta"))
+    mod <- MASS::polr(factor(y) ~ x1 + x2, data = dat)
+    mfx <- marginsxp(mod, 
+                     variance = NULL, 
+                     prediction_type = "probs")
+    mfx <- data.table(mfx)
+    ame <- mfx[, list(dydx = mean(dydx)), by = c("group", "term")][
+               , group := as.numeric(group)]
+    ame <- merge(ame, stata, by = c("group", "term"))
+    expect_equal(ame$dydx, ame$dydxstata, tolerance = 0.001)
+})
diff --git a/tests/testthat/test-dydx_betareg.R b/tests/testthat/test-dydx_betareg.R
index 98705a12d..da4a111a5 100644
--- a/tests/testthat/test-dydx_betareg.R
+++ b/tests/testthat/test-dydx_betareg.R
@@ -1,6 +1,9 @@
 skip_if_not_installed("betareg")
 
 library("margins")
+library("haven")
+library("betareg")
+library("data.table")
 
 test_that("betareg", {
     data("GasolineYield", package = "betareg")
@@ -11,3 +14,14 @@ test_that("betareg", {
     })
     expect_true(cor(as.numeric(mar$temp), res$temp, use = "complete.obs") > .99999)
 })
+
+
+test_that("betareg vs. Stata", {
+    stata <- readRDS(test_path("stata/stata.rds"))[["betareg_betareg_01"]]
+    dat <- read_dta(test_path("stata/data/betareg_betareg_01.dta"))
+    mod <- betareg::betareg(yield ~ factor(batch) + temp, data = dat)
+    mfx <- suppressWarnings(data.table(marginsxp(mod, variance = NULL)))
+    ame <- mfx[, list(dydx = mean(dydx)), by = "term"]
+    ame <- merge(ame, stata)
+    expect_equal(ame$dydx, ame$dydxstata, tolerance = 0.00001)
+})
diff --git a/tests/testthat/test-dydx_ivreg.R b/tests/testthat/test-dydx_ivreg.R
index dd91e50e6..add523129 100644
--- a/tests/testthat/test-dydx_ivreg.R
+++ b/tests/testthat/test-dydx_ivreg.R
@@ -1,6 +1,9 @@
 skip_if_not_installed("ivreg")
 
 library("margins")
+library("haven")
+library("data.table")
+library("ivreg")
 
 test_that("ivreg: vs. margins", {
     data(Kmenta, package = "ivreg")
@@ -9,3 +12,14 @@ test_that("ivreg: vs. margins", {
     mar <- data.frame(margins(mod, unit_ses = TRUE))
     marginsxp:::test_against_margins(res, mar, tolerance = .01)
 })    
+
+test_that("ivreg: vs. Stata", {
+    dat <- read_dta(test_path("stata/data/ivreg_ivreg_01.dta"))
+    stata <- readRDS(test_path("stata/stata.rds"))[["ivreg_ivreg_01"]]
+    mod <- ivreg::ivreg(Q ~ P + D | D + F + A, data = dat)
+    mfx <- marginsxp(mod)
+    mfx <- data.table(mfx)
+    ame <- mfx[, list(dydx = mean(dydx), std.error = mean(std.error)), by = "term"]
+    ame <- merge(ame, stata, by = "term")
+    expect_equal(ame$dydx, ame$dydxstata, tolerance = 0.0001)
+})
diff --git a/tests/testthat/test-dydx_stats.R b/tests/testthat/test-dydx_stats.R
index cbfc52be6..98507f4ab 100644
--- a/tests/testthat/test-dydx_stats.R
+++ b/tests/testthat/test-dydx_stats.R
@@ -1,4 +1,6 @@
 library("margins")
+library("haven")
+library("data.table")
 
 test_that("glm", {
     set.seed(1024)
@@ -17,6 +19,28 @@ test_that("glm", {
 })
 
 
+test_that("glm vs. Stata", {
+    stata <- readRDS(test_path("stata/stata.rds"))[["stats_glm_01"]]
+    dat <- read_dta(test_path("stata/data/stats_glm_01.dta"))
+    mod <- glm(y ~ x1 * x2, family = binomial, data = dat)
+    mfx <- data.table(marginsxp(mod))
+    ame <- mfx[, list(dydx = mean(dydx), std.error = mean(std.error)), by = "term"]
+    ame <- merge(ame, stata)
+    expect_equal(ame$dydx, ame$dydxstata, tolerance = 0.00001)
+})
+
+
+test_that("lm vs. Stata", {
+    stata <- readRDS(test_path("stata/stata.rds"))[["stats_lm_01"]]
+    dat <- read_dta(test_path("stata/data/stats_lm_01.dta"))
+    mod <- lm(y ~ x1 * x2, data = dat)
+    mfx <- data.table(marginsxp(mod))
+    ame <- mfx[, list(dydx = mean(dydx), std.error = mean(std.error)), by = "term"]
+    ame <- merge(ame, stata)
+    expect_equal(ame$dydx, ame$dydxstata, tolerance = 0.00001)
+})
+
+
 test_that("lm with interactions", {
     counterfactuals <- expand.grid(hp = 100, am = 0:1)
     mod <- lm(mpg ~ hp * am, data = mtcars)
@@ -27,14 +51,14 @@ test_that("lm with interactions", {
 })
 
 
-test_that("TODO: loess vcov error is raised too early to catch", {
+test_that("vcov(loess) does not exist", {
     mod <- loess(mpg ~ wt, data = mtcars)
     expect_error(marginsxp(mod), regexp = "not yet supported")
 })
 
 
 test_that("loess error", {
-    skip("not sure why I get different results for loess")
+    skip("loess produces different results under margins and marginsxp")
     mod <- loess(mpg ~ wt, data = mtcars)
     res <- marginsxp(mod, variance = NULL)
     mar <- data.frame(margins(mod))