Models and tests

Start with the expanded example for each command. Open the variants for options and alternative R approaches. Read the qualifications before adapting the code.

ameans

Compare arithmetic, geometric and harmonic means

Stata

ameans y

R stats

c(arithmetic=mean(data$y),geometric=exp(mean(log(data$y))),harmonic=1/mean(1/data$y))

Differs from Stata Returns the three point estimates for positive nonmissing values; Stata also reports their confidence intervals.

Stata reference · R reference

biprobit

Fit jointly correlated binary probit equations

Stata

biprobit y1 y2 x

R VGAM

VGAM::vglm(cbind(y1,y2)~x,family=VGAM::binom2.rho(),data=data)

Differs from Stata Both binary equations share x and use a constant latent correlation; rho link parameterizations differ but fitted joint probabilities have the same model meaning.

Stata reference · R reference

centile

Estimate a median and bootstrap its confidence interval

Stata

centile y, centile(50)

R boot

{ fit <- boot::boot(data,function(d,i) unname(stats::quantile(d$y[i],0.5,type=6)),R=499); list(median=fit$t0,interval=boot::boot.ci(fit,type="perc")$percent[4:5]) }

Point-estimate correspondence; this uses a percentile bootstrap interval, whereas Stata centile defaults to a binomial order-statistic interval.

Stata reference · R reference

contrast

Form group contrasts against the reference category

Stata

contrast r.g

R multcomp

multcomp::glht(fit,linfct=multcomp::mcp(g="Dunnett"))

Differs from Stata Returns contrast estimates and covariance. Stata contrast defaults to unadjusted tests; summary(glht) defaults to single-step multiplicity adjustment. Use summary(result, test=multcomp::adjusted(“none”)) for unadjusted p-values.

Stata reference · R reference

cpoisson

Fit a right-censored Poisson count model

Stata

cpoisson y x, ul(4)

R VGAM

VGAM::vglm(VGAM::SurvS4(y,status)~x,family=VGAM::cens.poisson(),data=data)

Differs from Stata Fixed right-censoring limit 4; status is one for exact counts and zero for counts known to be at least 4. Censoring must not be treated as truncation.

Stata reference · R reference

elasticnet

Fit cross-validated elasticnet linear prediction

Stata

elasticnet linear y x z, alpha(0.5)

R glmnet

glmnet::cv.glmnet(as.matrix(data[c("x","z")]),data$y,alpha=0.5,nfolds=10)

Prediction workflow: folds, standardization, penalty grids, loss scaling and lambda selection differ from Stata defaults; coefficients need not match.

Stata reference · R reference

esize

Estimate Cohen standardized mean differences for two independent groups

Stata

esize twosample y, by(g) cohensd

R effectsize

effectsize::cohens_d(y~g,data=data,pooled_sd=TRUE)

Differs from Stata Uses pooled SD and first-group minus second-group direction; reverse the sign if the substantive contrast is reversed.

Stata reference · R reference

esizei

Compute a pooled-SD effect size from group summaries

Stata

esizei 30 10 2 40 12 3, cohensd

R stats

(10-12)/sqrt(((30-1)*2^2+(40-1)*3^2)/(30+40-2))

Fixed summary-statistic workflow returns Cohen d only; confidence intervals and Hedges correction are not included.

Stata reference · R reference

icc

Estimate a one-way individual-measure intraclass correlation

Stata

icc y id

R psych

{ wide <- stats::reshape(data,idvar="id",timevar="rater",direction="wide"); psych::ICC(wide[,-1],lmer=FALSE)$results[1,,drop=FALSE] }

Differs from Stata Fixed balanced design has three measurements per target. R starts from wide repeated-measure columns; ICC1 is the one-way absolute-agreement individual-measure form.

Stata reference · R reference

lasso

Fit cross-validated lasso linear prediction

Stata

lasso linear y x z

R glmnet

glmnet::cv.glmnet(as.matrix(data[c("x","z")]),data$y,alpha=1,nfolds=10)

Prediction workflow: folds, standardization, penalty grids, loss scaling and lambda selection differ from Stata defaults; coefficients need not match.

Stata reference · R reference

manovatest

Test a factor term in a fitted multivariate model

Stata

manovatest g

R stats

summary(fit,test="Pillai")$stats["g",,drop=FALSE]

Differs from Stata After fitting the corresponding MANOVA model. This R row extracts Pillai; Stata manovatest reports four multivariate statistics.

Stata reference · R reference

nestreg

Compare successively added regression blocks

Stata

nestreg: regress y (x) (z)

R stats

{ reduced <- stats::lm(y~x,data=data); full <- stats::lm(y~x+z,data=data); stats::anova(reduced,full) }

Differs from Stata Compares these two nested OLS fits with a partial F test; it does not implement the general Stata prefix that sequences arbitrary model blocks.

Stata reference · R reference

orthog

Orthogonalize centered predictor columns in their stated order

Stata

orthog x z, generate(ox oz)

R stats

{ Q <- qr.Q(qr(cbind(1,as.matrix(data[c("x","z")])))); data[c("ox","oz")] <- Q[,-1,drop=FALSE]*sqrt(nrow(data)); data[c("ox","oz")] }

Differs from Stata Unweighted QR gives an orthogonal basis with mean square one; signs can differ from Stata modified Gram-Schmidt output.

Stata reference · R reference

orthpoly

Construct three orthogonal polynomial columns

Stata

orthpoly x, generate(p1 p2 p3) degree(3)

R stats

data[c("p1","p2","p3")] <- stats::poly(data$x,degree=3)*sqrt(nrow(data))

Differs from Stata Rescales the R orthonormal basis to mean square one; polynomial signs and stored coefficient conventions may differ.

Stata reference · R reference

predictnl

Apply a nonlinear prediction transform with delta-method standard errors

Stata

predictnl p = invlogit(predict(xb)), se(se)

R stats

{ pred <- stats::predict(fit,type="link",se.fit=TRUE); data$p <- stats::plogis(pred$fit); data$se <- data$p*(1-data$p)*pred$se.fit; data[c("p","se")] }

Differs from Stata Run after a logit GLM. This uses the derivative p*(1-p) and the fitted linear-index standard error.

Stata reference · R reference

pwcompare

Compare all group pairs with Tukey simultaneous intervals

Stata

pwcompare g, mcompare(tukey)

R stats

stats::TukeyHSD(stats::aov(y~g,data=data))$g

Differs from Stata After the one-way equal-variance model y on g; R TukeyHSD is not a general replacement for adjusted margins in arbitrary models.

Stata reference · R reference

pwmean

Compare group means using Bonferroni-adjusted pooled tests

Stata

pwmean y, over(g) mcompare(bonferroni)

R stats

stats::pairwise.t.test(data$y,data$g,p.adjust.method="bonferroni",pool.sd=TRUE)

Differs from Stata Uses a common pooled within-group variance and Bonferroni-adjusted pairwise t tests. The returned p-value matrix does not include Stata pairwise mean-difference confidence intervals.

Stata reference · R reference

robvar

Test equal group variances using median-centered Levene inference

Stata

robvar y, by(g)

R car

car::leveneTest(y~factor(g),data=data,center=stats::median)

Differs from Stata Corresponds to robvar W50; Stata additionally reports mean-centered W0 and 10-percent-trimmed W10.

Stata reference · R reference

runtest

Test whether a binary threshold sequence has random runs

Stata

runtest y, threshold(0)

R stats

{ above <- data$y>0; n1 <- sum(above); n0 <- sum(!above); runs <- 1+sum(diff(above)!=0); expected <- 1+2*n1*n0/(n1+n0); variance <- 2*n1*n0*(2*n1*n0-n1-n0)/((n1+n0)^2*(n1+n0-1)); z <- (runs-expected)/sqrt(variance); c(runs=runs,z=z,p.value=2*stats::pnorm(-abs(z))) }

Differs from Stata Uses the large-sample uncorrected runs statistic, matching Stata default continuity convention; missing values must be removed before this workflow.

Stata reference · R reference

slogit

Fit a rank-one stereotype multinomial logit model

Stata

slogit y x z, dimension(1) baseoutcome(3)

R VGAM

VGAM::rrvglm(y~x+z,family=VGAM::multinomial(refLevel="3"),data=data,Rank=1)

Differs from Stata Uses the reduced-rank multinomial-logit form of the one-dimensional Anderson stereotype model. Category score identification, scale and sign conventions can differ.

Stata reference · R reference

symmi

Test matched-table symmetry directly from cell counts

Stata

symmi 10 8 \ 3 10

R stats

stats::mcnemar.test(matrix(c(10,8,3,10),2,byrow=TRUE),correct=FALSE)

Differs from Stata For this 2-by-2 table the symmetry test is the uncorrected McNemar test. Stata also reports table summaries; this R result is an htest object.

Stata reference · R reference

tnbreg

Fit zero-truncated NB2 regression

Stata

tnbreg y x, ll(0)

R VGAM

VGAM::vglm(y~x,family=VGAM::posnegbinomial(),data=data)

Differs from Stata Zero-truncated NB2 with constant size; VGAM size equals 1/alpha. Information-matrix and optimizer conventions can differ.

Stata reference · R reference

tpoisson

Fit zero-truncated Poisson regression

Stata

tpoisson y x, ll(0)

R VGAM

VGAM::vglm(y~x,family=VGAM::pospoisson(),data=data)

Differs from Stata Zero-truncation only: fitted values are conditional positive-count means, while log-link coefficients parameterize the latent Poisson lambda.

Stata reference · R reference

truncreg

Fit a zero-left-truncated normal regression

Stata

truncreg y x, ll(0)

R VGAM

VGAM::vglm(y~x,family=VGAM::posnormal(),data=data)

Differs from Stata VGAM posnormal represents the special case of normal observations selected above zero; this example does not cover varying or upper truncation limits.

Stata reference · R reference

ztest

Test a mean using a known population standard deviation

Stata

ztest y == 20, sd(3)

R stats

{ z <- (mean(data$y)-20)/(3/sqrt(nrow(data))); c(z=z,p.value=2*stats::pnorm(-abs(z))) }

Differs from Stata The value 3 is assumed known as a population SD, not estimated from the sample.

Stata reference · R reference

ztesti

Compare independent means from known-variance summaries

Stata

ztesti 30 10 2 40 12 3

R stats

{ z <- (10-12)/sqrt(2^2/30+3^2/40); c(z=z,p.value=2*stats::pnorm(-abs(z))) }

Differs from Stata Treats the supplied standard deviations as known and uses a normal reference distribution. Returns the z statistic and two-sided p-value, without Stata confidence limits.

Stata reference · R reference