Statalist - Forums

wrong SE in dsregress with cluster() option

Adrien Bouguen — Sat, 03 Oct 2026 15:54:00 GMT

Hi statalisters,

it looks like when you run the following model:

Code:

dsregress Y T, controls(X1 X2...) cluster (clustervar)

dsregress does not generate an error and simply provide the double lasso estimator without correcting for clustering.

The correct code should be (as you probably know):

Code:

dsregress Y T, controls(X1 X2...) vce(cluster clustervar)

Yet, why does stata authorize the first coding? In a classical regression reg Y T, cluster(clustervar) and regression reg Y T, vce(cluster clustervar) provide the same SE.

I think this should not be allowed since it is a source of error.

What do you think?

Help with using Monte Carlo Simulation to do sensitivity power analysis on a mixed multilevel model

Apostolos Ng — Fri, 02 Oct 2026 02:42:07 GMT

Hi,
I am a psychology honours student who is using a mixed multilevel model for one of my research questions. I am considering using Monte Carlo simulation to do a sensitivity power analysis. However, we have not learn this at uni, and I am not sure if I am doing it correctly. Would someone be able to check if my coding is correct? It seems to run fine on my version of Stata 17, but as I am unfamiliar with it I just wanted to check if I have done the coding correctly. My variables are
pid: participant ID (level 2)
agegrp: age group, 1 = older, 2 = younger (older is reference level)
cval: cue valence, 1 = negative, 2 = neutral, 3 = positive
avg_rating_cue: affective appraisal rating (DV), averaged within each cue valence condition
prop_spec_val_c: proportion of specific memories in each valence condition (predictor, grand centred)

Step 1: Estimate inputs from the real data
* SD components of the predictor (valence main effect removed)
mixed prop_spec_val_c i.cval || pid:, reml
scalar sd_spec_b = exp(_b[lns1_1_1:_cons]) // between-person SD
scalar sd_spec_w = exp(_b[lnsig_e:_cons]) // within-person SD
scalar sd_spec_tot = sqrt(sd_spec_b^2 + sd_spec_w^2)

* Outcome model (observed)
mixed avg_rating_cue i.agegrp##c.prop_spec_val_c i.cval || pid:, reml
scalar sd_intercept = exp(_b[lns1_1_1:_cons]) // random-intercept SD
scalar sd_resid = exp(_b[lnsig_e:_cons]) // residual SD
scalar b_cons = _b[_cons]
scalar b_agegrp = _b[2.agegrp]
scalar b_spec = _b[c.prop_spec_val_c] // slope in the reference group (older)

Step 2: simulate one dataset and test the interaction
capture program drop sim_rq3
program sim_rq3, rclass
version 16
syntax, delta(real) n1(integer) n2(integer) nobs(integer) [alpha(real 0.05)]
local ntot = `n1' + `n2'
local reject = .
local conv = 0
quietly {
drop _all
set obs `ntot'
generate long pid = _n
generate byte agegrp = cond(_n <= `n1', 1, 2) // 1 = Older, 2 = Younger

generate double u_i = rnormal(0, sd_intercept) // random intercept
generate double s_i = rnormal(0, sd_spec_b) // person's mean specificity
expand `nobs' // 3 observations per person
generate double spec = s_i + rnormal(0, sd_spec_w)

* delta = true slope difference (added to group 1)
generate double slope = cond(agegrp==1, b_spec + `delta', b_spec)
generate double rating = b_cons + b_agegrp*(agegrp==2) ///
+ slope*spec + u_i + rnormal(0, sd_resid)

capture mixed rating i.agegrp##c.spec || pid:, reml iter(200)
if _rc == 0 {
if e(converged) == 1 {
local conv = 1
testparm i.agegrp#c.spec
local p = chi2tail(r(df), r(chi2)) // Wald test
local reject = (`p' < `alpha')
}
}
}
return scalar reject = `reject'
return scalar conv = `conv'
end

Step 3: power curve
local n_older = 34
local n_younger = 36
local nobs = 3
local reps = 1000
local seed = 12345

tempname results
postfile `results' double delta double delta_per_sd double power double conv_rate ///
using rq3_power_curve.dta, replace

foreach d of numlist 0 1.0(0.1)1.4 {
quietly simulate reject=r(reject) conv=r(conv), reps(`reps') seed(`seed') nodots: ///
sim_rq3, delta(`d') n1(`n_older') n2(`n_younger') nobs(`nobs')
quietly summarize reject
local pwr = r(mean)
quietly summarize conv
local cr = r(mean)
post `results' (`d') (`d'*sd_spec_tot) (`pwr') (`cr')
}
postclose `results'

Thank you for your help.

get AMEs by combining -margins- with some sort of parallel process?

Jonathan Horowitz — Thu, 01 Oct 2026 16:55:51 GMT

Dear Statalist:

I am running a massive multinomial logistic regression using restricted data on an organization's servers. I'm trying to get average marginal effects.

This is not my machine, I have to connect to it remotely and use their virtual machine. This means that upgrading to a Stata / MP that can use 16 (or 32, or 64) cores is not an option (unfortunately).

I had an idea to use the user-written package -parallel- by Yon and Quistorff to try and speed it up. I do not know if this is possible, but if it is theoretically possible I would like to try.

Some context:
This multinomial logistic regression has 60 outcomes, and the sample I'm running it on has over six million respondents. The multinomial logistic regression usually finishes in about 12-15 hours, and if I get predicted values using -margins- with no standard errors at all (-nose- option) it usually takes about 45 hours. I am currently trying to get predicted values using margins with standard errors and it is definitely using every last drop of RAM to do it--currently about about 110 gigabytes of RAM. I suppose it might even take more if I had the RAM to spare. I'm also not sure when it will finish. Maybe tomorrow, maybe Monday, maybe never.

However, something very different happened when I tried to get average marginal effects (AMEs) using -margins- and getting standard errors. The process started off by using up huge chunks of processing power (between 15% and 50% of the available CPU) and about 3 gigabytes of RAM for about a day. It then dropped down to about 5% to 20% of CPU, but increased to using about 6 gigabytes of RAM. And it kept doing this for about a week before I finally killed the process.

I suppose I can just set it up to run and come back and check it every week until it's done, but I'm not sure it will ever finish that way. I'm also not sure why it's not using more available RAM. I'm hoping that maybe there's a way to parallelize the AMEs further, chopping up my RAM into 6 gigabyte chunks across multiple instances. But that depends very much on how -margins- calculates the AME, and it is possible I am misunderstanding that part.

Has anyone tried doing something like this with -parallel-? Or anything else? If I'm patient enough eventually I can get predicted values without the SEs, but getting SEs for the AMEs seems to require something stronger.

Thanks for your time,
Jonathan

SSC activity, September 2026

KitBaum — Thu, 01 Oct 2026 02:23:19 GMT

During the last month, 28 packages have been added to the SSC Archive, and 42 packages have been updated, for a total of 4,014 packages.

At this week's Stata Conference in Boston, I will present details of the SSC-Next Generation project, which will transform the archive. One aspect of that project is already implemented: the ssc2 package (sec describe ssc2) provides access to earlier versions of community-contributed packages to enhance reproducibility. The Stata Conference RePEc series will share aspects of the SSC-NG project in the next few days.

Stata Sociology Virtual Symposium 2026

Samuel Mossberg — Wed, 30 Sep 2026 08:21:37 GMT

The 2026 Stata Sociology Virtual Symposium is a meeting of researchers in sociology from around the world discussing current theory and applied methods using Stata. The proceedings consist of invited talks by top Stata users in a virtual platform that allows you to experience this one-day event from wherever you are. The symposium is organized by StataNordic and Metrika Consulting.

Array

Program and Registration available at: https://www.statanordic.com/stata-so...mposium26.html

Winsorization of variables using winsor2

Bose Sudipta — Tue, 29 Sep 2026 12:04:57 GMT

I winsorized 25 variables simultaneously using the winsor2 command in Stata. However, I noticed that the first variable in the variable list does not appear to have been winsorized correctly and contains some unusual values after running the command. In contrast, the remaining variables appear to have been winsorized as expected.

Is this a known issue with winsor2, or could it be related to the way the command is specified? Has anyone encountered a situation where the first variable in a winsor2 command is treated differently from the others?

Any guidance on how to diagnose or resolve this issue would be greatly appreciated.

hmerge: a prototype hash-based -merge- (m:1, 1:1) that skips the master sort, benchmarks and request for testers

CJ Libassi — Fri, 25 Sep 2026 03:00:35 GMT

It occurred to me today that once upon a time I remember gtools having a long to-do list of speedups the author wanted to implement but never did, and that there were probably plenty of speedups we might get if we turned agents onto the problem of extending the work of gtools or finding other speedups. When I asked it to take a look at everything, it came up with a speedup to merge. What follows below is its writeup of those results, and I have confirmed the code runs on my own computer.

I have a prototype command, hmerge, that does merge m:1 and merge 1:1 with a hash table instead of sorting the master data. With the same syntax as merge and identical results, it runs about 8 times faster than merge on 10 million master observations with one integer key and up to 100,000 using observations (5.4 to 13.6 times across the key types I tried), when you don't need the result re-sorted by the merge key. The gain is smaller with very large using files. I'd love feedback, and especially people willing to try to break it.

Setup for everything below: Stata/MP 17.0 (2-core licence) on a MacBook Pro (Apple M5 Pro, 24 GB, macOS 26.5). hmerge is a C plugin plus an ado wrapper, so for now it only exists as a macOS arm64 build. Times are seconds, median of 3 runs with a warm file cache, end to end (master in memory, using file on disk, merged data in memory).

Why it's faster

merge sorts the whole master dataset by the key, moving every variable, before a linear sort-merge join. That sort is nearly all of the cost. With a master already sorted by the key, merge m:1 on 10 million observations takes 0.283s instead of 2.015s. hmerge instead reads the using data into a scratch frame, builds an index on its keys (a direct lookup table for one integer key over a compact range, otherwise a hash table), then looks up each master row where it already sits. The attached figure (hmerge_sortmerge_vs_hashjoin.png) shows the two approaches side by side. Hash collisions can't produce wrong matches, because a match is always confirmed by comparing the full key values.

Headline numbers

merge m:1 on one integer key, 10 million master observations in random order, 3 using variables:

Code:

using obs (J)    merge   hmerge   hmerge,sort   fmerge   frames   merge/hmerge
          100    1.924    0.231         1.929    1.018    7.070          8.3x
      100,000    2.015    0.259         2.028    1.002    7.447          7.8x
    5,000,000    2.299    1.204         3.370    2.910   11.449          1.9x

fmerge is from ftools (SSC), and frames means frame create + use + frlink + frget. At 100,000 using observations the speedup is 5.4x to 13.6x across key types (13.6x for a str12 key, 9.8x for a mixed string and numeric key, which fmerge doesn't accept).

Where it doesn't help

Code:

case (10 million obs, J = 100,000)                          merge   hmerge   ratio
master already sorted by the key                            0.283    0.239    1.2x
merge, then bysort on the key right after                   2.138    2.270    0.94x
panel sorted by person-year that has to keep that order     4.373    0.272    16.1x
1:1 on a unique integer key (J = 10 million)                2.377    1.602    1.5x

So the gain is skipping a sort you don't need. If you want the result sorted by the merge key anyway, merge is as fast or slightly faster. One thing about order: after merge m:1 the master observations are in key order, but Stata no longer treats the data as sorted, so a -by key:- afterwards needs a sort with either command. hmerge leaves the master observations in their original order. Both commands add using-only observations at the end, in key order.

A public benchmark

To check that this isn't an artifact of my own simulated data, I ran the join task from the H2O.ai/DuckDB Labs "database-like ops" benchmark (https://github.com/duckdblabs/db-benchmark, commit 3b074bc, with results for other software at https://duckdblabs.github.io/db-benchmark/). Its generator makes a 10 million observation table and three using tables. I mapped its five questions to merge with keepusing(v2) and added R's data.table on the same machine, limited to 2 threads to match Stata/MP:

Code:

question                         merge   hmerge   fmerge   data.table(2 thr)
q1 small inner (J=10)            1.681    0.323    0.878               0.104
q2 medium inner (J=10,000)       1.959    0.354    0.919               0.178
q3 medium outer (keep master)    1.861    0.163    0.719               0.222
q4 medium inner, string key      3.282    0.485    3.495               0.151
q5 big inner, 1:1 (J=10M)        4.182    1.621    6.635               1.103

All implementations return the same number of observations, and hmerge's result is identical to merge's on all five. The Stata times include reading the using .dta from disk, while data.table already has its tables in memory.

Correctness

- 130 differential test cases compare hmerge with merge on identical inputs: values, _merge, variable order, storage types (including merge's type promotion), formats, variable and value labels, error codes, the data left in memory after an error, and whether any sort flag left on the data is true.
- The tests also pass with the plugin compiled in a checking mode (UndefinedBehaviorSanitizer) that halts on any undefined behavior in the C code.

Other differences from merge: variable notes aren't copied, and r() holds r(path), which records how hmerge did the join, instead of what merge leaves there. For 1:m and m:m merges, the update, replace and force options, and strL variables, hmerge hands the job to merge and says so.

Limitations

It's only been run on Stata 17 with 2 cores on a Mac. Stata 18/19, Windows, Linux and more cores are all untested. To me the shakiest part is that the design relies on plugin memory persisting between two plugin calls, which works but isn't documented in the plugin interface, so a future Stata could break it (it is guarded so a mismatch errors out rather than using stale data). At a very large using file hmerge uses somewhat more memory than merge (594 vs 498 MB above the data itself with 5 million using observations).

What would help most

1. If you can build and run it on Windows or Linux, or on Stata 18/19 or more cores, and report timings.
2. Any case where hmerge and merge give different results.

Code, tests, benchmark do-files and a longer write-up: https://github.com/clibassi/hmerge (install: net install hmerge, from("https://raw.githubusercontent.com/clibassi/hmerge/main/")). A short demo do-file that builds a 10 million observation person-year panel and times both commands is in demo/hmerge_speedup_demo.do.

The index idea for integer keys follows gtools by Mauricio Cáceres Bravo (SSC, https://github.com/mcaceresb/stata-gtools), and ftools by Sergio Correia (SSC, https://github.com/sergiocorreia/ftools) was the obvious comparison. Claude built this with a little nudging from me.

If you see mistakes or notice any crazy choices, that would be super helpful. Thanks!

CJ Libassi

Array

Reducing data set

Euslaner — Thu, 24 Sep 2026 20:22:51 GMT

Is there a way to exclude cases from a dataset if they do not have sufficient observations?

I wrote:

keep house if _n(VCF0900c >= 45)

but it didn't work Is there a solution to this problem? Thanks for any help.

New Stata command: xtie2 – Interactive panel data estimation with automatic pairwise interactions

Zehra Yalniz — Thu, 24 Sep 2026 05:44:01 GMT

Dear Statalist members,

I would like to introduce xtie2, a new user-written Stata command for fixed-effects panel estimation with automatic pairwise interactions.

The command automatically generates all pairwise interactions among the explanatory variables and performs post-estimation analysis for the primary interaction defined by the first two explanatory variables. It reports the interaction coefficient and significance test, calculates marginal effects, evaluates them at selected moderator values, and computes marginal-effect zero-crossing critical values. It also checks whether these critical values fall within the observed sample range.

The reported critical values should be interpreted as marginal-effect zero-crossing values, not as Hansen-type panel threshold estimates.

The package is currently available from GitHub and can be installed in Stata using:

net install xtie2, from("https://raw.githubusercontent.com/zehrayalnizz41-cmd/xtie2/main/")

After installation:

help xtie2

Example:

xtie2 unemployment energy_dependency energy_price_index inflation gdp_growth, robust

GitHub repository:
https://github.com/zehrayalnizz41-cmd/xtie2

I would be grateful for comments, suggestions, and bug reports from Stata users.

Best regards,
Dr. Zehra Yalnız
Kocaeli University, Türkiye
ORCID: 0000-0003-2633-2022

relative survival poisson model

Matt Piercy — Wed, 23 Sep 2026 23:07:18 GMT

HI Statalist,
I am trying to model rs using a poisson model as per Royston and Lambert's book Flexible Paremetric Survival Analysis using Stata (2011). I have found the mean/variance is not equal. Does this matter? Is there a way of correcting for this in rs models?
Thanks,
Matt

Repeated cross-sectional LCA with multi-level data

Evie Gates — Wed, 23 Sep 2026 13:18:18 GMT

Hi there,

I want to conduct a latent class analysis of repeated cross-sectional data. Participants in the dataset were sampled and recruited at the school level; therefore, I would like to include a random effect for school clustering in the LCA; however, I am struggling to find a method to conduct a multilevel LCA in Stata. Does anybody have any recommendations on how to do this?

Many Thanks

PSA: Replacing string values does not entirely remove the old values from the saved .dta file

Nils Enevoldsen — Tue, 22 Sep 2026 22:49:04 GMT

I want to alert people who de-identify, anonymize, or share data to a behavior of Stata's .dta writer. I reported it to StataCorp Technical Support on 7 September. Details of their reply are at the end.

In short: replace name = "" does not remove names from a saved dataset. Neither does replacing names with shorter pseudonyms. The old values stay in the file, minus their first byte. Stata does not show them, but anyone can read them with a text editor or grep.

1. String data

When a str# value is replaced by a shorter one, Stata writes the new value and a null terminator, and leaves the rest of the fixed-width field as it was. The rest of the old value is then saved to disk. Here is a reproduction on StataNow 19.5:

Code:

version 19
clear
input str30 id
"Maria Fernanda Oliveira"
"João Carlos Pereira"
"Ana Beatriz Santos"
end
replace id = string(_n)
save pseudonymised.dta, replace

Stata shows id as 1, 2 and 3. On disk, the three 30-byte fields are

Code:

1\0ria Fernanda Oliveira\0\0\0\0\0\0\0
2\0ão Carlos Pereira\0\0\0\0\0\0\0\0\0\0
3\0a Beatriz Santos\0\0\0\0\0\0\0\0\0\0\0\0

Blanking the names with replace id = "" gives the same result. So does -saveold- to formats 117 and 115. -help dta- documents this at the level of the file format: a short string is followed by "random bytes". But those bytes are not random. They are the value you just overwrote.

Dropping the variable removes its contents from that save. However, freed string memory gets reused. Any string variable you create later in the same session, including one in a new dataset after -clear-, can pick up the dropped values and write them into the next file you save.

I scanned a few hundred published datasets from public research archives. About a quarter of those with a str# variable wider than some of its values have readable text in the padding that appears in no visible value in the file. The text includes addresses, place names and free-text responses. In one file, the padding held the full names of several individuals who appear nowhere in the dataset's visible values. The files most at risk are those with wide string variables that hold short values.

2. Metadata

The fixed-width metadata slots (variable names, formats, value label names, variable labels, sort list) are also written without being cleared first. The padding contains whatever memory Stata last used for them. That can include names and labels from a different dataset saved earlier in the same session. Using -clear-, -clear all- or -label drop _all- in between does not help. The padding also survives a load and re-save in a fresh session. About half of the published files I scanned have fragments like this. You can see it in the auto dataset Stata ships with: the slot for -turn- contains "turn\0ng_circle".

Code:

mata: fh = fopen(findfile("auto.dta"), "r"); s = fread(fh, 100000); fclose(fh); strpos(s, "ng_circle")

StataCorp's response

Technical Support told me that StataCorp does not consider this a bug, but they said it is prompting internal discussion. One option they are considering is zeroing the unused space in strings on -save-, either automatically or through an option, because doing it on every save could slow Stata down.

In my view, the current behavior is a disclosure risk for personally identifiable information (PII).

Updated egfr package on SSC

Phil Clayton — Mon, 21 Sep 2026 21:16:12 GMT

Thanks as always to Kit Baum, a new version of the egfr program is available on SSC. This program calculates eGFR from serum creatinine (and/or cystatin C) and is useful for people analysing datasets that include kidney function.

The main update is that it can now calculate the 2021 race-free CKD-EPI eGFR (with or without cystatin C). I also removed the default formula which was the old 4-variable MDRD formula, which is no longer appropriate as a default. Old code expecting that default can be updated if needed to specify that formula directly and results will be unchanged.

To install this package run:

Code:

ssc install egfr

If you already have it installed you can update it with:

Code:

adoupdate egfr, update

New package on SSC: bayesinits (chain-specific initial values for bayesmh)

Ben A. Dwamena — Mon, 21 Sep 2026 20:55:52 GMT

Thanks as always to Kit Baum, a new package -bayesinits- is now
available from the SSC archive.

To install:

. ssc install bayesinits
. help bayesinits

-bayesinits- (version 1.4.1) generates domain-aware, randomized,
reproducible initial values for one or more Markov chains, for use
with -bayesmh-'s init#() options. Running multiple chains from
dispersed, feasible starting points is the basis of the Gelman-Rubin
convergence diagnostic (-bayesstats grubin-), but constructing such
values by hand is tedious and error-prone when parameters are
constrained: a carelessly drawn variance may be nonpositive, a
correlation may fall outside (-1, 1), and mixture weights may fail to
sum to one.

-bayesinits- samples on domain-appropriate working scales and
transforms back, so every chain starts at a valid point. Six parameter
domains are supported:

unconstrained() normal draws on the raw scale
positive() log-scale draws, exponentiated
corr() Fisher-z draws, mapped back via tanh()
fisherz() draws directly on the z scale
phi() angular parameters in (0, pi)
simplex() normalized log-normal weights (positive, sum 1)

Chain-specific strings are returned in r(init1), ..., r(initK),
together with composite strings r(init_all) and r(init_all_wrap) that
can be passed to -bayesmh- verbatim. Default centers and scales may be
overridden per parameter via center() and scale(), which accept
numeric expressions. The mle option derives centers and scales
automatically from the estimation results in memory: centers at the
point estimates, scales as inflate() (default 2) times the
delta-method standard error on the working scale, with map() linking
Bayesian parameter names to e(b) coefficient names. seed() and
rngstream() provide reproducibility and independent streams for
parallel chains.

A typical workflow:

. bayesinits, nchains(3) unconstrained(b0 b1) positive(sig2)
. bayesmh y = ({b0} + {b1}*x), likelihood(normal({sig2})) ///
prior({b0}, normal(0,100)) prior({b1}, normal(0,100)) ///
prior({sig2}, igamma(0.01,0.01)) ///
nchains(3) `r(init_all)' rseed(12345)
. bayesstats grubin

Or, anchored at a preliminary fit with threefold inflation:

. regress y x
. bayesinits, nchains(4) unconstrained(b0 b1) positive(sig2) ///
mle map(b0=_cons b1=x) inflate(3) center(sig2=e(rmse)^2)

-bayesinits- requires Stata 15 or later. It is particularly useful for
hierarchical models with correlation structures (for example,
bivariate diagnostic test accuracy meta-analysis), mixture models, and
any setting where dispersed multi-chain runs and credible convergence
assessment matter.

Comments, bug reports, and suggestions are welcome.

Ben Adarkwa Dwamena, MD
Division of Nuclear Medicine and Molecular Imaging
Department of Radiology, University of Michigan
bdwamena@umich.edu

Stata module to download BCRP (Central Reserve Bank of Peru) data: dbcrp

Anthony Huaynate — Mon, 21 Sep 2026 19:52:28 GMT

Dear Statalist community,

I would like to announce the release of dbcrp, a simple Stata package I built to download macroeconomic time series directly from the API of the Central Reserve Bank of Peru (BCRP). As an economics student, I often needed a quick and clean way to get BCRP data into Stata without messy workarounds. This package streamlines the process.

To install:

Code:

ssc install dbcrp

Quick example using monthly inflation (PN01271PM):

Code:

clear all
dbcrp PN01271PM, names(ipc)
tsline ipc

Further options include date ranges, frequency selection, and batch downloads of multiple series (the API supports up to 20 series per request).

GitHub: https://github.com/Anthony-maxed/STATA-dbcrp
LinkedIn: https://www.linkedin.com/in/anthony-...fre-1a3897263/

Feedback and suggestions are very welcome.

Best regards,
Anthony Huaynate Onofre
Universidad ESAN