Parmby gives 'normal' instead of 'logit' confidence intervals when using proportion command

[Warren Holroyd]

Afternoon all,

I have a question about the excellent parmby module, provided by Roger Newson.

As a simple example, I'm running the following command:

clear all
frames reset
sysuse auto
proportion rep78

As expected, the citype on proportion defaults to logit and yields the following results:

Proportion estimation             Number of obs   =         69

--------------------------------------------------------------
             |                                   Logit
             | Proportion   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       rep78 |
          1  |   .0289855   .0201966      .0070794    .1110924
          2  |    .115942   .0385422       .058317    .2173648
          3  |   .4347826   .0596787      .3214848    .5553295
          4  |   .2608696   .0528625      .1695907    .3788629
          5  |   .1594203   .0440694      .0895793     .267702
--------------------------------------------------------------

Running the same code in parmby...

parmby "proportion rep78", frame(output1, replace)
frame output1: list, noobs clean

...yields a different set of confidence intervals:

    parmseq      parm    estimate      stderr   dof           t           p        min95       max95  
          1   1.rep78   .02898551   .02019662    68    1.435166   .15582355   -.01131623   .06928724  
          2   2.rep78   .11594203   .03854218    68   3.0081856   .00368254    .03903231   .19285175  
          3   3.rep78   .43478261   .05967869    68   7.2853911   4.333e-10    .31569563   .55386958  
          4   4.rep78   .26086957    .0528625    68   4.9348699   5.465e-06    .15538409   .36635504  
          5   5.rep78   .15942029   .04406936    68   3.6174863    .0005668    .07148126   .24735932  

Parmby returns confidence intervals for the proportion as if I had run the proportion command with the argument citype(normal).

I have attempted the following which have not solved the problem:

• explicitly added citype(logit) to the parmby command
• adjusted the available parmby confidence interval options
• run commands with set trace on to ensure that no unseen error is occurring
• read Mr Newson's paper on confidence intervals to see if there's anything obvious I'm missing

So, specific questions:

• Has anyone experienced this or can explain why it's happening?
• Is this expected behaviour for parmby?
• Does anyone know of a way of "forcing" parmby to pick up the logit transformed CIs, rather than the normal/t-distribution CIs?

Thanks for considering my questions.

Warren

[Warren Holroyd]

Any takers for this one?

I know it's a fairly niche area, given that not everybody is a user or parmby.

Forgot to mention that I'm using Stata16 and the latest version of parmby/parmest.

[Carlo Lazzaro]

Warren:
let's hope that Roger Newson will chime in.

[Warren Holroyd]

For the purposes of closure and helping future users, I should mention that I think I found a roundabout solution after being prompted by a colleague on a related matter.

Essentially, I manually calculate the logit transferred CIs, based on the examples provided by Leonardo Guizzetti in his post in the following thread:

Standard errors and 95% Confidence Intervals for Proportions - Differences between different Stata versions - Statalist

So, using an amendment to my original example:

clear all
frames reset
sysuse auto
proportion rep78

parmby "proportion rep78", frame(output1, replace)

frame output1 {
    gen min95_lt = invlogit(logit(estimate) - invt(dof, 0.975) * stderr/(estimate*(1-estimate)))
    gen max95_lt = invlogit(logit(estimate) + invt(dof, 0.975) * stderr/(estimate*(1-estimate)))
    }
frame output1: list, noobs clean

Which yields:

   parmseq      parm    estimate      stderr   dof           t           p        min95       max95   min95_lt   max95_lt  
          1   1.rep78   .02898551   .02019662    68    1.435166   .15582355   -.01131623   .06928724   .0070794   .1110924  
          2   2.rep78   .11594203   .03854218    68   3.0081856   .00368254    .03903231   .19285175    .058317   .2173648  
          3   3.rep78   .43478261   .05967869    68   7.2853911   4.333e-10    .31569563   .55386958   .3214848   .5553295  
          4   4.rep78   .26086957    .0528625    68   4.9348699   5.465e-06    .15538409   .36635504   .1695907   .3788629  
          5   5.rep78   .15942029   .04406936    68   3.6174863    .0005668    .07148126   .24735932   .0895793    .267702  

This is assumes that:

• you're using Stata v14+ (or wish to emulate it's method)
• your parmby output uses the default resultsset variable names; and
• that you're after 95% confidence.

You'd need to adjust accordingly otherwise.

In my case, we're collecting proportions and other calculations in a batch script, meaning that not all CIs will need to be logit transferred.

In such a case, one could also choose to replace the existing CI variables by enabling the command option in parmby, then conditioning a replacement with a substring match, like:

clear all
frames reset
sysuse auto
proportion rep78

parmby "proportion rep78", frame(output1, replace) command

frame output1 {
    replace min95 = invlogit(logit(estimate) - invt(dof, 0.975) * stderr/(estimate*(1-estimate))) if index(command,"proportion")
    replace max95 = invlogit(logit(estimate) + invt(dof, 0.975) * stderr/(estimate*(1-estimate))) if index(command,"proportion")
    }
frame output1: list, noobs clean

Yielding:

    parmseq            command      parm    estimate      stderr   dof           t           p       min95       max95  
          1   proportion rep78   1.rep78   .02898551   .02019662    68    1.435166   .15582355   .00707941   .11109242  
          2   proportion rep78   2.rep78   .11594203   .03854218    68   3.0081856   .00368254   .05831701   .21736479  
          3   proportion rep78   3.rep78   .43478261   .05967869    68   7.2853911   4.333e-10   .32148479   .55532951  
          4   proportion rep78   4.rep78   .26086957    .0528625    68   4.9348699   5.465e-06   .16959074   .37886295  
          5   proportion rep78   5.rep78   .15942029   .04406936    68   3.6174863    .0005668   .08957931   .26770201  

If you have variable names with the word "proportion" in them, you might need to rename them.

Hope it helps someone.

Regards,
Warren