Hacky workaround for prtest and pwcorr following MI - are these equivalent?

[Paul Sanfilippo]

I am helping a student with an analysis of a dataset with a fair bit of missing data and am exploring using multiple imputation.

There's nothing too elaborate in the analysis - the main components are to test for differences in proportions and calculate correlations. To that end prtest (one sample) and pwcorr don't work with mi estimate that I can see.

For the former is there an equivalence to running an intercept only logistic model and for the latter, could I effectively run a regress on the standardised versions of the variables (to obtain the same results as prtest and pwcorr, but on multiply imputed datasets)? Like in the following (just on the one dataset to illustrate):

sysuse auto, clear

* prtest vs logit-margins
prtest foreign == 0.5
logit foreign // p-value is test of difference of log-odds from 0?
margins // p-value is test of difference of probability from 0.5?

* pwcorr vs regress using z-scores
pwcorr weight length
egen weight_z = std(weight)
egen length_z = std(length)
reg weight_z length_z
reg length_z weight_z

Thanks,

Paul

[Felix Bittmann]

For the prtest I am not sure whether this gives the same results. For the correlation you could easily use mibeta.

clear all
version 16.1
set seed 1234
sysuse auto


pwcorr mpg price
reg mpg price, beta

*Create missing data*
gen random = runiform()
sort random
replace mpg = . in 1/10
replace price = . in 20/30


mi set flong
mi register imputed mpg price
mi impute chained (regress) mpg price, add(10)

mibeta mpg price    //Standardized coefficient = -0.455

[Joro Kolev]

I think for the -prtest- logit is not equivalent.

For one thing, -prtest- is exact, -logit- is approximate, so they should not be equivalent.

For another, they are not equivalent in the example Original Poster suggests:

. sysuse auto, clear
(1978 Automobile Data)

. prtest foreign==.5

One-sample test of proportion                   Number of obs      =        74

------------------------------------------------------------------------------
    Variable |       Mean   Std. Err.                     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     foreign |   .2972973   .0531331                      .1931583    .4014363
------------------------------------------------------------------------------
    p = proportion(foreign)                                       z =  -3.4874
Ho: p = 0.5

     Ha: p < 0.5                 Ha: p != 0.5                   Ha: p > 0.5
 Pr(Z < z) = 0.0002         Pr(|Z| > |z|) = 0.0005          Pr(Z > z) = 0.9998

. logit foreign, noheader

Iteration 0:   log likelihood =  -45.03321  
Iteration 1:   log likelihood =  -45.03321  
------------------------------------------------------------------------------
     foreign |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |  -.8602013   .2543331    -3.38   0.001    -1.358685   -.3617176
------------------------------------------------------------------------------

In logit we are modelling log[p/(1-p)] = _cons, also log(1)= 0, therefore the prtest above is equivalent to the test that the intercept is 0.

The p-values of the two test are different Pr(|Z| > |z|) = 0.0005 vs p-value logit = 0.001.

[Joro Kolev]

For the linear case (my machine does not recognise the -mibeta- command Felix mentions), bivariate regression with beta weights is equivalent to pwcorr, e.g.,

. pwcorr price turn, sig

             |    price     turn
-------------+------------------
       price |   1.0000 
             |
             |
        turn |   0.3096   1.0000 
             |   0.0073
             |

.  reg price turn, noheader beta
------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|                     Beta
-------------+----------------------------------------------------------------
        turn |   207.5794   75.12938     2.76   0.007                 .3096174
       _cons |  -2064.987   2996.813    -0.69   0.493                        .
------------------------------------------------------------------------------

[Felix Bittmann]

I somehow thought that since version 16 mibeta is included in the main installation, maybe I am wrong. However, it is an official ado by Stata Corp.

findit mitbeta

[daniel klein]

I somehow thought that since version 16 mibeta is included in the main installation, maybe I am wrong. However, it is an official ado by Stata Corp.

Actually, anything that is not included/distributed as part of the Stata installation (or update) is not officially approved by StataCorp. Yulia Marchenko, who wrote mibeta, happens to work for StataCorp, but that does not make mibeta part of official Stata. I am not stressing this to be a smart ass; I am stressing this because I believe there is a reason why mibeta is not part of official Stata. The reason being that there simply is not an agreed-on way how to get standardized coefficients from multiply imputed data (see, e.g., Van Ginkel, 2020).



Van Ginkel, J.R. 2020. Standardized Regression Coefficients and Newly Proposed Estimators for R2 in Multiply Imputed Data. Psychometrika, 85, 128--205.

[Felix Bittmann]

To say it with Lothar Matthäus: again what learned. Thanks Daniel for the clarification

[Federico Tedeschi]

For the linear case (my machine does not recognise the -mibeta- command Felix mentions), bivariate regression with beta weights is equivalent to pwcorr, e.g.,

. pwcorr price turn, sig

| price turn
-------------+------------------
price | 1.0000
|
|
turn | 0.3096 1.0000
| 0.0073
|

. reg price turn, noheader beta
------------------------------------------------------------------------------
price | Coef. Std. Err. t P>|t| Beta
-------------+----------------------------------------------------------------
turn | 207.5794 75.12938 2.76  0.007  .3096174
_cons | -2064.987 2996.813 -0.69 0.493 .
------------------------------------------------------------------------------

I think this is an alternative to "pwcorr" to get the point estimate. To get the level of statistical significance, however, one should do something different, because the multiply-imputed nature of the dataset must be taken into account (any command line not specifying that the data are multiply imputed will lead Stata to consider observations as independent, thus inflating sample size, leading to an over-estimation of statistical significance).