[Bolade Banougnin]

Thank you so much Andrew for your brilliance. I appreciate your time. Sure I will reach out to Ben Jann on this.

[Bolade Banougnin]

Please post images in .PNG format as recommended in the FAQs. I am not able to open the .GPH image created by Stata 17 in Stata 16. That being the case, I ran the code in #14 and three p-values were missing. [ATTACH=CONFIG]n1679353[/ATTACH]



On further investigation, I note that these values are not missing in the stored results (highlighted below).

. est replay SexLastUnsafeConsist0

--------------------------------------------------------------------------------------
Model SexLastUnsafeConsist0
--------------------------------------------------------------------------------------

Mixed-effects logistic regression Number of obs = 1,174
Group variable: ID Number of groups = 587

Obs per group:
min = 2
avg = 2.0
max = 2

Integration method: mvaghermite Integration pts. = 7

Wald chi2(.) = .
Log likelihood = . Prob > chi2 = .
-------------------------------------------------------------------------------------
SexLastUnsafeCons~t | Coef. Std. Err. t P>|t| [95% Conf. Interval]
--------------------+----------------------------------------------------------------
XC | -.6959997 .6446828 -1.08  0.280 -1.959555 .5675554
XM | .3778831 .3858202 0.98 0.327 -.3783106 1.134077
Rural | -.6248963 .5701184 -1.10 0.273 -1.742308 .4925153
AgeC | .4053572 .1053465 3.85 0.000 .1988818 .6118326
HomeTypeInformal | -.5606224 .7523278 -0.75 0.456 -2.035158 .9139131
NecessitiesAll_r | .5613208 .474195 1.18 0.237 -.3680843 1.490726
ARTever | -.8648473 .479 -1.81 0.071 -1.80367 .0739755
Relationship | 2.024283 .5228796 3.87 0.000 .999458 3.049108
SchEnrol | -2.499274 1.336178 -1.87 0.061 -5.118134 .1195862
j | .8185097 .2963594 2.76 0.006 .2376559 1.399363
_cons | -2.661133 1.376298 -1.93 0.053 -5.358626 .0363608
--------------------+----------------------------------------------------------------
var(_cons[ID])| 8.865928 3.155539 4.413333 17.81073
-------------------------------------------------------------------------------------
LR test vs. logistic model: chi2(.) = . Prob > chi2 = .

Note: LR test is conservative and provided only for reference.

. est replay SexLastUnsafeConsist1

--------------------------------------------------------------------------------------
Model SexLastUnsafeConsist1
--------------------------------------------------------------------------------------

Mixed-effects logistic regression Number of obs = 1,532
Group variable: ID Number of groups = 766

Obs per group:
min = 2
avg = 2.0
max = 2

Integration method: mvaghermite Integration pts. = 7

Wald chi2(.) = .
Log likelihood = . Prob > chi2 = .
-------------------------------------------------------------------------------------
SexLastUnsafeCons~t | Coef. Std. Err. t P>|t| [95% Conf. Interval]
--------------------+----------------------------------------------------------------
XC | -.6728635 .3554466 -1.89 0.058 -1.369526 .023799
XM | .3563452 .215419 1.65 0.098 -.0658684 .7785587
Rural | .2718676 .235966 1.15 0.249 -.1906173 .7343525
AgeC | .5020056 .061067 8.22 0.000 .3823165 .6216948
HomeTypeInformal | .1759705 .2726909 0.65 0.519 -.3584938 .7104348
NecessitiesAll_r | -.0812407 .2505907 -0.32 0.746 -.5723895 .4099081
ARTever | -.670258 .2284569 -2.93 0.003 -1.118025 -.2224908
Relationship | 1.025744 .2646242 3.88 0.000 .5070905 1.544398
SchEnrol | .1936175 .3525943 0.55 0.583 -.4974547 .8846896
j | 1.00031 .1776544 5.63 0.000 .652114 1.348506
_cons | -3.3551 .4931144 -6.80 0.000 -4.321587 -2.388614
--------------------+----------------------------------------------------------------
var(_cons[ID])| 2.449458 .6544243 1.450954 4.135104
-------------------------------------------------------------------------------------
LR test vs. logistic model: chi2(.) = . Prob > chi2 = .

Note: LR test is conservative and provided only for reference.

Furthermore, if I asked for the coefficients to be printed, they are available. Since you have 1500+ observations, my workaround involves approximating the p-values using the normal distribution in these cases where they are missing. I can only conclude that this is a bug in coefplot. I would recommend that you send a dataset and some code that reproduces the issue to Ben Jann so that he may have a look at it. The approximation is highlighted in the code below.

use "C:\Users\amus\Downloads\20220825 sample data.dta"

mi set mlong // mi set data for imputation

mi reshape wide SexDrunkDrugs SexLastUnsafeConsist SexPartnerFreq2More InequitSex Access XC XM Use, i(ID) j(j)


* Identifying which variables in the imputation model have missing information.
mi register imputed NecessitiesAll_r HomeTypeInformal SexDrunkDrugs0 SexDrunkDrugs1 InequitSex0 InequitSex1

*** Mutiple imputation
*MI using chained equations/MICE (also known as the fully conditional specification or sequential generalized regression) for binary outcomes. In simulation studies (Lee & Carlin, 2010; Van Buuren, 2007), MICE has been show to produce estimates that are comparable to MVN method.

* Clear any xtset that is ongoing
mi xtset, clear
mi /*mimpt*/ impute chained (logit) NecessitiesAll_r HomeTypeInformal SexDrunkDrugs0 ///
SexDrunkDrugs1 InequitSex0 InequitSex1 = Access0 Access1 ///
XC0 XM0 XC1 XM1 Sex ARTever AgeC Rural Relationship ///
SchEnrol, add(10) rseed(20210101) /*skipnonconvergence(5)*/ savetrace(tracedata, replace)

* Bring data back to long format
mi reshape long SexDrunkDrugs SexLastUnsafeConsist SexPartnerFreq2More InequitSex Access XC XM Use, i(ID) j(j)

* We can check the data and see if there is no problems in our imputation
summ SexDrunkDrugs SexLastUnsafeConsist SexPartnerFreq2More InequitSex HomeTypeInformal NecessitiesAll_r
local y SexDrunkDrugs SexLastUnsafeConsist SexPartnerFreq2More InequitSex
local x0 Access
local x1 XC XM
local x Rural AgeC Sex HomeTypeInformal NecessitiesAll_r ARTever Relationship SchEnrol j
local xsex Rural AgeC HomeTypeInformal NecessitiesAll_r ARTever Relationship SchEnrol j
foreach ys of local y {
forvalues f = 0/1 {
cap noi mi estimate, post saving(`ys'`f', replace) cmdok : melogit `ys' `x1' `xsex' if Sex == `f' || ID:, or cov(un)
cap noi estimate store `ys'`f'
}
}
*set trace on
set scheme s1color
cap drop xpos
gen xpos = 3.5
coefplot (SexLastUnsafeConsist0, label("Boy") msymbol(S) mcolor(black) mfcolor(white)) ///
(SexLastUnsafeConsist1, label("Girl") msymbol(O) mcolor(black) mfcolor(black)), ///
bylabel("Infrequent condom use") mlabel(cond(@pval<0.001, "{it:p} < 0.001", ///
cond(@pval>0.001 &@pval<., "{it:p} = " + string(@pval,"%9.3f"), ///
 cond(missing(@pval) & 2*normal(@t)<0.001, "{it:p} < 0.001", "{it:p} = " + string(2*normal(@t),"%9.3f")) ))) ///
mlabc(none) addplot(scatter @at xpos, ms(i) mlab(@mlbl) mlabcolor(black) ///
mlabpos(9) mlabgap(-10) mlabsize(small)) graphr(margin(r=10)) xsc(r(-10 5))

[ATTACH=CONFIG]n1679354[/ATTACH]

Hi again Andrew,

In the event where the t-statistics is equal to a value says 0.92, we have a p-value that is higher than 1. Maybe they should be another condition for if 2*normal(@t)>1 where the p-value could be evaluated to 2*normal(-@t). I tried Laplace integration method and was confronted to that specific problem. The updated code is:

cap drop xpos
gen xpos = 8.5
coefplot    (SexDrunkDrugs0x0 SexDrunkDrugs0, label("Boy") msymbol(S) mcolor(black) mfcolor(white))     ///
            (SexDrunkDrugs1x0 SexDrunkDrugs1, label("Girl") msymbol(O) mcolor(black) mfcolor(black)),     ///
            bylabel("Sex after substance use")                                                             ///
      ||    (SexLastUnsafeConsist0x0 SexLastUnsafeConsist0)                                                ///
            (SexLastUnsafeConsist1x0 SexLastUnsafeConsist1)                                            ,    ///
            bylabel("Infrequent condom use")                                                             ///
      ||    (SexPartnerFreq2More0x0 SexPartnerFreq2More0)                                                ///
            (SexPartnerFreq2More1x0 SexPartnerFreq2More1)                                            ,     ///
            bylabel("Multiple sexual partnership")                                                        ///
      ||     (InequitSex0x0 InequitSex0)                                                                    ///
            (InequitSex1x0 InequitSex1)                                                                ,     ///
            bylabel("Inequitable sexual partnership")                                                    ///
      ||  , eform keep(Access XC XM) ciopts(lcolor(gs8) lwidth(thick)) xtitle("Adjusted odds ratios (ORs)") xlab(0(1)10,format(%3.0fc)) xline(1, lcolor(red)) ylab(1 `" "Access" "({it:vs. } no access)" "' 2 `" "Health content" "use ({it:vs. } no access)" "' 3 `" "Social media" "use ({it:vs. } no access)" "', labsize(medsmall)) byopts(xrescale) subtitle(, fcolor(none) lstyle(none)) legend(order(2 "ORs for boys" 4 "ORs for girls" 3 "95% CI") row(1) ring(0) pos(10) region(lstyle(none)))             ///
      mlabel(cond(@pval<0.001, "{it:p}<0.001", ///
                cond(@pval>0.001 &@pval<., "{it:p}=" + string(@pval,"%9.3f"), ///
                    cond(missing(@pval) & 2*normal(@t)>1, "{it:p}=" + string(2*normal(-@t),"%9.3f"),    ///
                        cond(missing(@pval) & 2*normal(@t)<0.001, "{it:p}<0.001", "{it:p}=" + string(2*normal(@t),"%9.3f")) )))) ///
  mlabc(none) addplot(scatter @at xpos, ms(i) mlab(@mlbl) mlabcolor(black) ///
  mlabpos(9) mlabgap(-10) mlabsize(small)) graphr(margin(r=10)) xsc(r(0 10))

This yields the p-values in the regression table outputs.

[Andrew Musau]

In the event where the t-statistics is equal to a value says 0.92, we have a p-value that is higher than 1.

Good point! Another way to do it is

2*normal(-abs(@t))