Twopm- post estimation

[Talar Markossian]

I am interested to calculate predicted mean expenditures by condition group (conditions 1-4). The outcome is total expenditures including some 0 values, so I fit the two-part model. Main interest is to use margins to estimate predicted mean for each condition category. One of the conditions (condition 3) has 0 (no) observations with 0 expenditures and that is why the logit model coefficient is 0, and no predicted margin was calculated for condition 3. Can I still use the twopm, or any suggestions for this situation since the main interest is predicted mean for conditions?

svy: twopm TOTEXP22 i.diag2 if exclusion==1, firstpart(logit, nolog) secondpart(glm, family(gamma) link(log) nolog)
(running twopm on estimation sample)

Survey data analysis

Number of strata = 105 Number of obs = 4,023
Number of PSUs = 348 Population size = 49,972,532
Design df = 243
F(2, 242) = 5.23
Prob > F = 0.0059

------------------------------------------------------------------------------
| Linearized
TOTEXP22 | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
logit |
diag2 |
condition2 | -.2460358 .2919973 -0.84 0.400 -.8212046 .329133
condition3 | 0 (empty)
condition4 | 1.22004 .4498981 2.71 0.007 .3338425 2.106238
|
_cons | 3.447114 .1897759 18.16 0.000 3.073298 3.820929
-------------+----------------------------------------------------------------
glm |
diag2 |
condition2 | .1630623 .0944394 1.73 0.086 -.0229621 .3490866
condition3 | .8484066 .051607 16.44 0.000 .7467524 .9500608
condition4 | .6767057 .0865197 7.82 0.000 .5062814 .8471299
|
_cons | 9.467698 .051607 183.46 0.000 9.366044 9.569353
------------------------------------------------------------------------------

. margins i.diag2

Adjusted predictions

Number of strata = 105 Number of obs = 4,023
Number of PSUs = 348 Population size = 49,972,532
Subpop. no. obs = 3,973
Subpop. size = .
Model VCE: Linearized Design df = 243

Expression: twopm combined expected values, predict()

------------------------------------------------------------------------------
| Delta-method
| Margin std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
diag2 |
condition1 | 12535.97 661.786 18.94 0.000 11232.4 13839.54
condition2 | 14630.29 1188.604 12.31 0.000 12289.01 16971.57
condition4 | 25211.33 1478.899 17.05 0.000 22298.24 28124.43
------------------------------------------------------------------------------

[George Ford]

Been noodling on this a minute. See if this helps.

clear all
webuse womenwk , clear
replace wage = 0 if wage==.
g wagem = wage > 0

****************************************************************************************************
** PROGRAM TO GET MARGINS
****************************************************************************************************
capture program drop my_twopart_margin
program define my_twopart_margin, rclass
    args dv dum edlev
    
    // Check for zeros in this education group (in original data)
    qui count if educ==`edlev' & `dv'==0
    local has_zeros = (r(N) > 0)
    
    if `has_zeros' {
        // Two-part model
        qui logit `dum' i.educ age married children
        // Don't replace! Use margins at() instead
        qui margins, at(educ=`edlev') nose
        matrix b1 = r(b)
        local pr = b1[1,1]
        
        qui glm `dv' i.educ age married children if `dum'==1, family(gamma) link(log)
        qui margins, at(educ=`edlev') nose
        matrix b2 = r(b)
        local mu = b2[1,1]
        
        return scalar margin = `pr' * `mu'
    }
    else {
        // No zeros: just second stage
        qui glm `dv' i.educ age married children if `dum'==1, family(gamma) link(log)
        qui margins, at(educ=`edlev') nose
        matrix b2 = r(b)
        
       return scalar margin = b2[1,1]
    }
end

****************************************************************************************************
** DATA AS IS
****************************************************************************************************
** TWOOM
twopm wage i.educ age married children , firstpart(logit) secondpart(glm, family(gamma) link(log))
estimates store e1
margins i.educ ,  post

// Run for each education level
foreach lev in 10 12 16 20 {
    bootstrap margin=r(margin), reps(50) seed(123): my_twopart_margin wage wagem `lev'
    di _n "Education `lev':"
    estat bootstrap, all
}

****************************************************************************************************
** CHANGE DATA TO EXCLUDE EDUC == 16
****************************************************************************************************
g wage2 = wage
replace wage2 = 0.001 if educ==16 & wage==0
g wagem2 = wage2 > 0

** TWOPM
twopm wage2 i.educ age married children , firstpart(logit) secondpart(glm, family(gamma) link(log))
estimates store e1
margins i.educ ,  post

// Run for each education level
foreach lev in 10 12 16 20 {
    bootstrap margin=r(margin), reps(50) seed(123): my_twopart_margin wage2 wagem2 `lev'
    di _n "Education `lev':"
    estat bootstrap, all
}

[George Ford]

Delta Method (just use simple margins from no 0 model)

** First Stage
logit wagem i.educ age married children , r
estimates store e2

** Second Stage  
glm wage i.educ age married children if wagem==1, family(gamma) link(log) vce(robust)
estimates store e3

** Delta method approach
foreach lev in 10 12 16 20 {
    
    // Check for zeros in this education group
    qui count if educ==`lev' & wage==0
    local has_zeros = (r(N) > 0)
    
    if `has_zeros' {
                  
        // First stage - logit
        qui estimates restore e2
        qui margins, at(educ=`lev') vce(unconditional) post
        matrix b1 = e(b)
        matrix V1 = e(V)
        local pr = b1[1,1]
        local var_pr = V1[1,1]
        
        // Second stage - GLM
        qui estimates restore e3
        qui margins, at(educ=`lev') vce(unconditional) post
        matrix b2 = e(b)
        matrix V2 = e(V)
        local mu = b2[1,1]
        local var_mu = V2[1,1]
        
        // Delta method margins: Var(pr*mu) ≈ mu^2*Var(pr) + pr^2*Var(mu)
        local margin = `pr' * `mu'
        local var_margin = (`mu'^2 * `var_pr') + (`pr'^2 * `var_mu')
        local se_margin = sqrt(`var_margin')
        
        // z-stat and p-value
        local z_stat = `margin' / `se_margin'
        local p_value = 2*(1 - normal(abs(`z_stat')))
        
        // 95% CI
        local ci_lower = `margin' - 1.96 * `se_margin'
        local ci_upper = `margin' + 1.96 * `se_margin'
        
        if `lev' == 10 {
        di "Educ Level" _col(15) "Margins" _col(32) "SE" _col(42) "Z" _col(50) "p-val" _col(70) "95% CI"
        di `lev' _col(15) %6.3f `margin' _col(30) %6.3f `se_margin' _col(40) %6.3f `z_stat' _col(50) %6.4f `p_value' _col(65) "[" %6.3f `ci_lower' ", " %6.3f `ci_upper' "]"
        } 
        else {
        di `lev' _col(15) %6.3f `margin' _col(30) %6.3f `se_margin' _col(40) %6.3f `z_stat' _col(50) %6.4f `p_value' _col(65) "[" %6.3f `ci_lower' ", " %6.3f `ci_upper' "]"
        }
        
    }
    else {
        // No zeros: just use second stage margins 
        di "Education `lev' (no zeros - second stage only):"
        qui estimates restore e3
        margins, at(educ=`lev') vce(unconditional)
    }
}

[Talar Markossian]

Hi, this is extremely helpful and thank you for addressing my question, I appreciate the help. I ran the Delta Method (second approach) and it works perfectly. Do you have any thoughts whether the Delta Method (second approach) would still apply after running the logit and glm models with svy: (Survey design) - any reasons why not? I did quick search and it looks like it's okay. Thank you.

[George Ford]

I'd study that a bit. I think you'd need vce(unconditional) in the margins command to get the right SE. I don't work much with svy data. You can pretty close to the right result by using weights and clustering on PSU, depending on how complex the survey design is.

[Talar Markossian]

Ok, yes, thank you for the help!