Multivariate and psychometric methods

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

alpha

Estimate raw Cronbach alpha without reversing items

Stata

alpha q1 q2 q3, asis

R Package:psych

psych::alpha(data[c("q1","q2","q3")],check.keys=FALSE)

Differs from Stata asis disables Stata item-direction reversal; check.keys=FALSE makes the R direction rule explicit.

Generate a mean score from a reliability scale

Stata

alpha q1 q2 q3, generate(score) asis

R Package:base

score <- rowMeans(data[c("q1","q2","q3")])

Differs from Stata asis disables Stata item-direction reversal; check.keys=FALSE makes the R direction rule explicit.

Assess reliability when each item is deleted

Stata

alpha q1 q2 q3, item asis

R Package:psych

psych::alpha(data[c("q1","q2","q3")],check.keys=FALSE)$alpha.drop

Differs from Stata asis disables Stata item-direction reversal; check.keys=FALSE makes the R direction rule explicit.

Estimate standardized Cronbach alpha

Stata

alpha q1 q2 q3, std asis

R Package:psych

psych::alpha(data[c("q1","q2","q3")],check.keys=FALSE)$total$std.alpha

Differs from Stata asis disables Stata item-direction reversal; check.keys=FALSE makes the R direction rule explicit.

ca

Perform simple correspondence analysis of a contingency table

Stata

ca g h

R Package:MASS

MASS::corresp(unclass(table(data$g,data$h)),nf=2)

Differs from Stata MASS returns standard coordinates; Stata ca defaults to normalize(symmetric), so plotted coordinates differ in scale. Singular values and inertias are comparable, while axes may reverse signs. The table retains category names.

camat

Perform correspondence analysis of a supplied count matrix

Stata

camat C, dimensions(2)

R Package:MASS

MASS::corresp(C,nf=2)

Differs from Stata MASS returns standard correspondence coordinates; Stata defaults to symmetric normalization. Singular values are comparable, but coordinate scaling and axis signs can differ.

candisc

Estimate canonical linear discriminant directions

Stata

candisc x y z w, group(g)

R Package:MASS

MASS::lda(g ~ x+y+z+w,data=data,prior=rep(1/3,3))

Differs from Stata MASS lda supplies discriminant directions with equal class priors. This is a discriminant-analysis workflow, not the full Stata candisc canonical-test and normalization output.

canon

Estimate canonical correlations between two variable sets

Stata

canon (x y) (z w)

R Package:stats

stats::cancor(data[c("x","y")],data[c("z","w")])

Differs from Stata Returns canonical correlations and coefficient matrices for the two variable blocks. Canonical signs and coefficient normalization are arbitrary; Stata also reports significance tests.

cluster

Fit averagelinkage hierarchical clustering on standardized variables

Stata

cluster averagelinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="average")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

8 more examples and alternatives

Fit centroidlinkage hierarchical clustering on standardized variables

Stata

cluster centroidlinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="centroid")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

Fit completelinkage hierarchical clustering on standardized variables

Stata

cluster completelinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="complete")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

Cut a hierarchical cluster tree into three groups

Stata

cluster generate group = groups(3)

R Package:stats

stats::cutree(fit,k=3)

Differs from Stata Cuts the existing hclust tree at three groups. The result is an integer membership vector; cluster labels are arbitrary and Stata generates a named dataset variable.

Partition standardized observations into three k-means clusters

Stata

cluster kmeans x y z w, k(3)

R Package:stats

stats::kmeans(scale(as.matrix(data[c("x","y","z","w")])),centers=3,nstart=10)

Differs from Stata Algorithms, starts and random streams differ, so cluster numbering and local optima may differ.

Fit medianlinkage hierarchical clustering on standardized variables

Stata

cluster medianlinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="median")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

Fit singlelinkage hierarchical clustering on standardized variables

Stata

cluster singlelinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="single")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

Fit wardslinkage hierarchical clustering on standardized variables

Stata

cluster wardslinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="ward.D2")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

Fit waveragelinkage hierarchical clustering on standardized variables

Stata

cluster waveragelinkage x y z w

R Package:stats

stats::hclust(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),method="mcquitty")

Differs from Stata Standardize variables before the analogous Stata call. Ward distance and squaring conventions require care.

clustermat

Cluster a precomputed dissimilarity matrix

Stata

clustermat completelinkage D

R Package:stats

stats::hclust(stats::as.dist(D),method="complete")

Differs from Stata Uses complete linkage on the supplied dissimilarity matrix. Dissimilarity definitions and tie breaking affect the tree; R returns an hclust object rather than a Stata named cluster solution.

discrim

Classify observations using five nearest neighbors

Stata

discrim knn x y z w, group(g) k(5)

R Package:class

class::knn(train=as.matrix(data[c("x","y","z","w")]),test=as.matrix(data[c("x","y","z","w")]),cl=data$g,k=5)

Differs from Stata Selects predictors by name and classifies the training observations. class::knn uses unweighted majority voting and random tie breaking; Stata defaults to equal class priors, which can change decisions for unbalanced classes. Use held-out data for predictive accuracy.

Fit LDA discriminant analysis

Stata

discrim lda x y z w, group(g)

R Package:MASS

MASS::lda(g ~ x+y+z+w,data=data,prior=rep(1/3,3))

Differs from Stata Prior class probabilities and covariance estimation conventions should be specified before comparing results.

Classify using leave-one-out linear discriminant analysis

Stata

discrim lda x y z w, group(g) lootable

R Package:MASS

MASS::lda(g ~ x+y+z+w,data=data,CV=TRUE,prior=rep(1/3,3))

Stata lootable requests leave-one-out classification; MASS CV=TRUE performs the corresponding prediction workflow.

Use equal prior class probabilities in LDA

Stata

discrim lda x y z w, group(g) priors(equal)

R Package:MASS

MASS::lda(g ~ x+y+z+w,data=data,prior=rep(1/3,3))

Differs from Stata Equal class priors are explicit where shown. MASS coefficient normalization and class order can differ from Stata; compare posterior classifications after aligning priors and covariance estimates.

Fit QDA discriminant analysis

Stata

discrim qda x y z w, group(g)

R Package:MASS

{ analysis <- droplevels(stats::model.frame(g ~ x + y + z + w, data = data, na.action = stats::na.omit)); groups <- length(unique(analysis$g)); MASS::qda(g ~ x + y + z + w, data = analysis, prior = rep(1 / groups, groups)) }

Differs from Stata Explicit equal priors are computed from classes remaining after complete-case filtering of the response and predictors, matching Stata discrim qda default priors(equal). MASS otherwise uses sample class proportions. Align covariance normalization and class labels before comparing posterior probabilities.

estat

Report confirmatory-factor model fit indices

Stata

estat gof, stats(all)

R Package:lavaan

lavaan::fitMeasures(fit,c("chisq","df","cfi","tli","rmsea","srmr"))

Differs from Stata Selects six lavaan SEM fit indices. Model identification, likelihood scaling, missing-data estimation and baseline-model definitions must match; this is not every statistic in Stata stats(all).

4 more examples and alternatives

Assess factorability using Kaiser-Meyer-Olkin measures

Stata

estat kmo

R Package:psych

psych::KMO(stats::cor(X))

Differs from Stata Computes KMO measures from the supplied correlation matrix. R returns overall and per-variable measures in a list; match the same correlation and missing-data sample used for Stata factor.

Inspect SEM modification indices

Stata

estat mindices

R Package:lavaan

lavaan::modindices(fit)

Differs from Stata Returns lavaan modification indices for admissible freed parameters. Align model constraints, information matrix and missing-data likelihood; candidate parameter sets can differ from Stata.

Inspect residual correlations after SEM

Stata

estat residuals

R Package:lavaan

lavaan::lavResiduals(fit,type="raw")

Differs from Stata Uses raw observed-minus-implied covariance residuals, matching Stata default raw residuals. lavaan also returns standardized summaries with its own scaling, which are not the same Stata output columns.

Compute squared multiple correlations for scale items

Stata

estat smc

R Package:stats

1-1/diag(solve(stats::cor(X)))

Differs from Stata Computes squared multiple correlations as one minus the reciprocal diagonal of the inverse correlation matrix. Requires a nonsingular correlation matrix and returns a vector rather than Stata display tables.

factor

Fit a one-factor maximum-likelihood factor model

Stata

factor q1 q2 q3 q4 q5 q6, factors(1) ml

R Package:stats

stats::factanal(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),factors=1,rotation="none")

Differs from Stata Matches the stated ML or principal-component extraction and factor count. Loading signs and rotations are indeterminate; psych principal-component scores and scaling need not reproduce Stata factor postestimation.

3 more examples and alternatives

Fit iterated principal-axis factors

Stata

factor q1 q2 q3 q4 q5 q6, factors(2) ipf

R Package:psych

psych::fa(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),nfactors=2,fm="pa",rotate="none")

Differs from Stata Principal-axis algorithms and convergence tolerances vary; do not assume identical loadings.

Fit a two-factor maximum-likelihood factor model

Stata

factor q1 q2 q3 q4 q5 q6, factors(2) ml

R Package:stats

stats::factanal(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),factors=2,rotation="none")

Differs from Stata Matches the stated ML or principal-component extraction and factor count. Loading signs and rotations are indeterminate; psych principal-component scores and scaling need not reproduce Stata factor postestimation.

Fit a principal-component factor solution

Stata

factor q1 q2 q3 q4 q5 q6, factors(2) pcf

R Package:psych

psych::principal(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),nfactors=2,rotate="none")

Differs from Stata Matches the stated ML or principal-component extraction and factor count. Loading signs and rotations are indeterminate; psych principal-component scores and scaling need not reproduce Stata factor postestimation.

factormat

Fit maximum-likelihood factors from a correlation matrix

Stata

correlate q1 q2 q3 q4 q5 q6
matrix R = r(C)
factormat R, n(300) factors(2) ml

R Package:stats

stats::factanal(covmat=R,n.obs=300,factors=2,rotation="none")

Differs from Stata Fixed workflow. R is a fixed correlation-matrix object; sample size supplies likelihood-based fit statistics.

hotelling

Test a vector of means against zero with Hotelling T squared

Stata

hotelling x y

R Package:stats

{ n <- nrow(data); k <- 2; m <- colMeans(data[c("x","y")]); T2 <- n*drop(t(m)%*%solve(stats::cov(data[c("x","y")]),m)); F <- (n-k)/(k*(n-1))*T2; c(T2=T2,F=F,p.value=stats::pf(F,k,n-k,lower.tail=FALSE)) }

Differs from Stata Computes the one-sample Hotelling T-squared and its F transformation for these two variables. Assumes a complete sample and invertible covariance; Stata additionally displays mean summaries.

kap

Estimate unweighted Cohen kappa for two raters

Stata

kap rater1 rater2

R Package:psych

psych::cohen.kappa(table(data$rater1,data$rater2))$kappa

Differs from Stata Returns unweighted Cohen kappa from the cross-tabulation. It omits Stata agreement summaries, standard errors and confidence intervals; align missing-value and category handling.

manova

Test multivariate group differences using Pillai

Stata

manova y z = g

R Package:stats

summary(stats::manova(cbind(y,z)~g,data=data),test="Pillai")

Differs from Stata Stata displays all four MANOVA statistics; this R example selects the named statistic.

Fit multivariate analysis of covariance

Stata

manova y z = g c.x

R Package:stats

stats::manova(cbind(y,z)~g+x,data=data)

Differs from Stata Fits the stated multivariate linear model. R summary defaults to Pillai while Stata displays several multivariate criteria; choose the same criterion, hypothesis and factor contrasts.

3 more examples and alternatives

Test multivariate group differences using Wilks

Stata

manova y z = g

R Package:stats

summary(stats::manova(cbind(y,z)~g,data=data),test="Wilks")

Differs from Stata Stata displays all four MANOVA statistics; this R example selects the named statistic.

Test multivariate group differences using Hotelling-Lawley

Stata

manova y z = g

R Package:stats

summary(stats::manova(cbind(y,z)~g,data=data),test="Hotelling-Lawley")

Differs from Stata Stata displays all four MANOVA statistics; this R example selects the named statistic.

Test multivariate group differences using Roy

Stata

manova y z = g

R Package:stats

summary(stats::manova(cbind(y,z)~g,data=data),test="Roy")

Differs from Stata Stata displays all four MANOVA statistics; this R example selects the named statistic.

mca

Analyze multiple categorical variables through indicator MCA

Stata

mca g h j, method(indicator)

R Package:MASS

MASS::mca(data, nf=2)

Differs from Stata MASS uses an indicator-matrix MCA; inertia normalization differs from adjusted and joint MCA.

mds

Fit classical metric multidimensional scaling

Stata

mds x y z w, id(id) dimension(2) noplot

R Package:stats

stats::cmdscale(stats::dist(as.matrix(data[c("x","y","z","w")])),k=2,eig=TRUE)

Differs from Stata Uses classical metric scaling of Euclidean distances. Stata mds defaults and stress criteria can differ; coordinates are comparable only after aligning distances, method, rotation and reflection.

Fit nonmetric multidimensional scaling

Stata

mds x y z w, id(id) method(nonmetric) std noplot

R Package:MASS

MASS::isoMDS(stats::dist(scale(as.matrix(data[c("x","y","z","w")]))),k=2,trace=FALSE)

Differs from Stata Nonmetric optimization, stress definitions and starting configurations can differ.

mdsmat

Scale an existing dissimilarity matrix

Stata

mdsmat D, dimension(2) noplot

R Package:stats

stats::cmdscale(stats::as.dist(D),k=2)

Differs from Stata Applies classical metric scaling to the supplied dissimilarity matrix. Match Stata metric method and eigenvalue treatment; coordinate rotations and reflections are arbitrary.

mvreg

Fit multivariate multiple regression

Stata

mvreg y z = x w

R Package:stats

stats::lm(cbind(y,z) ~ x+w,data=data)

Differs from Stata Returns a multivariate lm fit on the common complete-case sample. Stata mvreg displays separate equations and cross-equation summaries; request the corresponding R summaries explicitly.

mvtest

Compute Mardia multivariate skewness and kurtosis tests

Stata

mvtest normality x y z w, stats(skewness kurtosis)

R Package:psych

psych::mardia(X,plot=FALSE)

Differs from Stata Finite-sample corrections may differ; this selects Mardia components rather than Stata default Doornik-Hansen.

pca

Fit correlation-matrix PCA

Stata

pca x y z w

R Package:stats

stats::prcomp(as.matrix(data[c("x","y","z","w")]),center=TRUE,scale.=TRUE)

Differs from Stata Selects correlation or covariance PCA explicitly. Component signs are arbitrary; prcomp returns singular-value-based objects and scores whose scaling differs from Stata displays.

2 more examples and alternatives

Retain two principal components

Stata

pca x y z w, components(2)

R Package:stats

stats::prcomp(as.matrix(data[c("x","y","z","w")]),scale.=TRUE,rank.=2)

Differs from Stata Selects correlation or covariance PCA explicitly. Component signs are arbitrary; prcomp returns singular-value-based objects and scores whose scaling differs from Stata displays.

Fit covariance-matrix PCA

Stata

pca x y z w, covariance

R Package:stats

stats::prcomp(as.matrix(data[c("x","y","z","w")]),center=TRUE,scale.=FALSE)

Differs from Stata Selects correlation or covariance PCA explicitly. Component signs are arbitrary; prcomp returns singular-value-based objects and scores whose scaling differs from Stata displays.

pcamat

Perform PCA directly on a correlation matrix

Stata

correlate x y z w
matrix R = r(C)
pcamat R, n(50)

R Package:stats

stats::princomp(covmat=R)

Differs from Stata Fixed workflow.

predict

Predict posterior probabilities for all discriminant classes

Stata

predict p1 p2 p3, pr

R Package:MASS

data[c("p1","p2","p3")] <- stats::predict(fit)$posterior

Differs from Stata Uses the existing discriminant fit and its class order. Returned labels or posterior-probability columns must match the fitted classes; R factor order need not match Stata category codes.

4 more examples and alternatives

Compute the first two principal-component scores

Stata

predict pc1 pc2, score

R Package:stats

data[c("pc1","pc2")] <- stats::predict(fit)[,1:2,drop=FALSE]

Differs from Stata Returns scores from the existing PCA fit. R score scaling and arbitrary component signs can differ from Stata predict, score; align normalization before comparing values.

Predict discriminant-analysis class membership

Stata

predict predicted, classification

R Package:MASS

data$predicted <- stats::predict(fit)$class

Differs from Stata Uses the existing discriminant fit and its class order. Returned labels or posterior-probability columns must match the fitted classes; R factor order need not match Stata category codes.

Compute Bartlett-method factor scores

Stata

predict score1 score2, bartlett

R Package:stats

data[c("score1","score2")] <- stats::factanal(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),factors=2,rotation="none",scores="Bartlett")$scores

Differs from Stata Fixed measurement model with indicators q1-q6; output score names can be changed.

Compute regression-method factor scores

Stata

predict score1 score2, regression

R Package:stats

data[c("score1","score2")] <- stats::factanal(as.matrix(data[c("q1","q2","q3","q4","q5","q6")]),factors=2,rotation="none",scores="regression")$scores

Differs from Stata Fixed measurement model with indicators q1-q6; output score names can be changed.

rotate

Use the oblique oblimin loading-rotation criterion

Stata

rotate, oblique oblimin

R Package:GPArotation

GPArotation::oblimin(unclass(fit$loadings))

Differs from Stata Rotates the supplied loading matrix only. Loading signs, axis order and normalization can differ from Stata; correlated-factor solutions also require comparing the factor-correlation matrix.

Use the orthogonal quartimax loading-rotation criterion

Stata

rotate, orthogonal quartimax

R Package:GPArotation

GPArotation::quartimax(unclass(fit$loadings))

Differs from Stata Rotates the supplied loading matrix only. Loading signs, axis order and normalization can differ from Stata; correlated-factor solutions also require comparing the factor-correlation matrix.

Rotate factor loadings using promax

Stata

rotate, promax

R Package:stats

stats::promax(fit$loadings,m=3)

Differs from Stata Loadings are identified only up to sign and ordering; rotation normalization conventions may differ.

Rotate factor loadings using varimax

Stata

rotate, varimax

R Package:stats

stats::varimax(fit$loadings,normalize=FALSE)

Differs from Stata Loadings are identified only up to sign and ordering; rotation normalization conventions may differ.

rotatemat

Rotate an externally supplied loading matrix

Stata

rotatemat L, varimax

R Package:stats

stats::varimax(L)

Differs from Stata Applies varimax to the supplied matrix. Match normalization and column ordering; equivalent rotated solutions may reverse column signs.

sem

Fit a three-factor confirmatory measurement model

Stata

sem (visual -> x1 x2 x3) (textual -> x4 x5 x6) (speed -> x7 x8 x9), latent(visual textual speed)

R Package:lavaan

lavaan::cfa("visual =~ x1+x2+x3; textual =~ x4+x5+x6; speed =~ x7+x8+x9",data=data)

Differs from Stata Latent-variable scale constraints and mean structures must be aligned when comparing estimates.

1 more examples and alternatives

Report standardized SEM parameters

Stata

sem, standardized

R Package:lavaan

lavaan::standardizedSolution(fit)

Differs from Stata Fixed workflow.

tetrachoric

Estimate a tetrachoric correlation by maximum likelihood

Stata

tetrachoric a b

R Package:polycor

polycor::polychor(data$a,data$b,ML=TRUE,std.err=TRUE)

Differs from Stata Each binary variable must have two observed levels; inference assumes latent bivariate normality.