Analyzing censored data with the Tobit model

This replicates Table 2.2 in Breen (1996). 

clear

// Simulate data
set obs 2000
generate ui = rnormal(0, 2)
generate xi = rnormal()
generate yi_star = 1 + 2*xi + ui
drop ui

// Fit OLS model
regress yi_star xi
eststo ols1
estadd scalar sigma = e(rmse)

// Truncate variable
generate yi = yi_star if yi_star > 0
replace  yi = 0       if yi_star <= 0
recode   yi (0 = .), gen(yi_h) 

// (A) OLS (using all observations on y
//     including y1 = 0)
regress yi xi
eststo ols2
estadd scalar sigma = e(rmse)

// (B) OLS (yi > 0) 
regress yi xi if yi >0
eststo ols3
estadd scalar sigma = e(rmse)

// (C) Heckman 2-step
heckman yi_h xi, select(xi) twostep 
eststo heckman

// (D) Tobit
tobit yi xi, ll(0)
eststo tobit
estadd scalar sigma = _b[sigma:_cons]

// Table 2.2


esttab ols2 ols3 heckman tobit ols1, b(3) se(3) nostar stat(sigma) ///
       mtitles("(A) OLS incl. yi = 0" ///
               "(B) OLS yi > 0" ///
               "(C) Heckman 2-step" ///
               "(D) Tobit" ///
               "OLS yi_star") ///
       coeflabel(_cons "alpha" xi "beta") ///
       collabels() ///
       drop(mills:lambda sigma:_cons) ///
       order(_cons xi) ///
       unstack ///
       modelwidth(20) nonumber

Reference

Breen, Richard. 1996. Regression Models. Censored, Sample Selected, or Truncated Data. Sage. doi: 10.4135/9781412985611

Selection of regression predictors

clear
set seed 1
set obs 1000

// Generate random variable y
generate y = rnormal()

// Generate 50 random variables x
forvalues i = 01/50 {
 generate x`i' = rnormal()
}

// Model 1: Regress y on the x's
quietly regress y x1-x50
estimates store model1
coefplot model1, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 1") ysize(8) xsize(3) ///
                 drop(_cons) msymbol(o) name(model1, replace)

// Identify variables significant at 10% level
forvalues i = 1/50 {
 local t = _b[x`i'] / _se[x`i']
 local p = 2 * ttail(e(df_r), abs(`t'))
 if `p' <= .10 {
   local significant10 `significant10' x`i'
 }
 di "x`i'" _skip(5) `t' _skip(5) `p' _skip(5) "`significant10'"
}

// Model 2: Regress y on the x's significant at the 10 per cent level
quietly regress y `significant10'
estimates store model2
coefplot model2, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 2") ysize(8) xsize(3) ///
                 drop(_cons) msymbol(o) name(model2, replace)
    

// Identify variables significant at 25% level
estimates restore model1
forvalues i = 1/50 {
 local t = _b[x`i'] / _se[x`i']
 local p = 2 * ttail(e(df_r), abs(`t'))
 if `p' <= .25 {
   local significant25 `significant25' x`i'
 }
 di "x`i'" _skip(5) `t' _skip(5) `p' _skip(5) "`significant25'"
}

// Model 3: Regress y on the x's significant at the 25 per cent level
quietly regress y `significant25'
estimates store model3
coefplot model3, xline(0) xscale(alt) ylabel(, labsize(*.7)) ///
                 xtitle("Regression weights and 90% CI's", size(*.8)) levels(90) ///
                 xlabel(, format(%6.2f) labsize(*.8)) ///
                 title("Model 3") ysize(8) xsize(3)  ///
                 drop(_cons) msymbol(o) name(model3, replace)

graph combine model1 model2 model3, row(1) xcommon

Interactions with and without main effects

set scheme s1mono
clear

set obs 10000
set seed 1

gen a = rnormal()
gen b = rnormal()
gen e = 2*rnormal()

gen y = 2 + 3*a + 4*b + 5*a*b + e

summarize y a b

estimates clear
regress y c.a##c.b
estimates store M1
regress y c.a  c.b
estimates store M2
regress y      c.b c.a#c.b
estimates store M3
regress y c.a      c.a#c.b
estimates store M4
regress y          c.a#c.b
estimates store M5

estimates restore M1
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M1, replace) noci ///
             title("M1: Main terms and interaction") legend(col(1))

estimates restore M2
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M2, replace) noci ///
             title("M2: Main terms, no interaction") legend(off)

estimates restore M3
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M3, replace) noci ///
             title("M3: Interaction, one main term missing") legend(off)

estimates restore M4
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M4, replace) noci ///
             title("M4: Interaction, other main term missing") legend(off)

estimates restore M5
margins, at(a=(-1 0 1) b=(-1 0 1))
marginsplot, name(M5, replace) noci ///
             title("M5: Interaction term, both main terms missing") legend(off)

grc1leg M1 M2 M3 M4 M5, col(2) ysize(11) xsize(8) ring(0) pos(5) ycommon
estimates table M1 M2 M3 M4 M5, b(%9.2f) stat(F r2_a)

Calculate inflection point of curvilinear relationships

// Generate data
clear
set seed 1
set obs 500
generate  e = rnormal(0,10)     
generate  x = rnormal(0,2)     
generate  y  = (-2*x + 2*(x*x) + e)
regress y c.x##c.x

matrix list e(b) // Find out names of coefficients

// Plot regression coefficients
coefplot, xline(0) ///
xtitle(" " "Regression coefficients and 95% CI's") ///
drop(_cons) scheme(s1mono) ///
coeflabels(x = "x" c.x#c.x = "x{char 178}") ///
ciopts(recast(rcap)) ///
headings(x = "U-shaped effect of x") ///
name(coeffs, replace)

// Calculate bend
local bend = round(-_b[x]/(2*_b[c.x#c.x]), .01)
display "b1 = " round(_b[x], .01) _newline ///
        "b2 = " round(_b[c.x#c.x], .01) _newline ///
        "-b1/(2*b2) = " `bend'

predict yhat  // Generate predicted values

// Plot results
twoway (line yhat x, sort(x)) ///
       (scatter y x) ///
   , legend(off) ///
    xline(`bend') ///
    ytitle("y") ///
    xtitle("x") ///
    name(scatter, replace) ///
    title("Inflection point at `bend'")
    

Random graphs (46): Plotting the contours of interactions

Next to the classical layout for interaction plots, it is also possible to make use of contour plots for visualizing interactions.

// Simulate data
clear
set seed 1
set obs 100

generate  e  = 0 + (100000 - 0) * runiform()  // To generate random variates over the
generate  x1 = 0 + (600 - 0) * runiform()     // interval [a,b), a+(b-a)*runiform()
generate  x2 = 0 + (600 - 0) * runiform()  
generate  y  = (x1 + x2 + (x1*x2) + e) / 1000

// Calculations for Aiken & West (1991) style plot
qui regress y c.x1##c.x2

qui sum x1
local x1_minsd = r(mean) - r(sd)
local x1_mean  = r(mean)
local x1_plusd = r(mean) + r(sd)

qui sum x2
local x2_minsd = r(mean) - r(sd)
local x2_mean  = r(mean)
local x2_plusd = r(mean) + r(sd)

// Calculate predicted values for plot
margins, at(c.x2 = (`x2_minsd' `x2_mean' `x2_plusd') ///
            c.x1 = (`x1_minsd' `x1_mean' `x1_plusd')) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
   xlabel(`x2_minsd' "-1 SD" `x2_mean' "Average x2" `x2_plusd' "+1 SD") ///
   title("Aiken & West (1991)-style plot for interactions") ///
   ytitle("Predicted y") ///
   xtitle("") ///
   plotopts(msymbol(none)) ///        // Turn off markers
   plot1opts(lpattern(longdash)) ///  // Define line types here
   plot2opts(lpattern(solid)) ///
   plot3opts(lpattern(shortdash)) ///
   legend(subtitle(x1) ///
   order(4 "- 1 SD" 5 "Average x1" 6 "+ 1 SD" 2 "95% CI")) ///
   name(lines, replace)
        // For some reason, marginsplot ignores the -label option-, therefore
        // -order- is used here

// Calculations for contour plot
qui regress y c.x1##c.x2
predict y_hat
// Plot
qui twoway contour y_hat x1 x2, ///
       title("Contour plot for interactions") ///
    ztitle("Predicted y") ///
    name(contours, replace)
// Combine figures
graph combine lines contours, col(1) ysize(8)

Random graphs (27): Between and within regression

// Case A
clear
set seed 1
set obs 300

generate  e = 0 + (200 - 0) * runiform()  // To generate random variates over the
generate  x = 0 + (600 - 0) * runiform()  // interval [a,b), a+(b-a)*runiform()

generate  y = e // Relationship between X and Y random
generate id = 1
replace   y = (e + 200) if x > 200 & x <= 400
replace  id = 2         if x > 200 & x <= 400
replace   y = (e + 400) if x > 400
replace  id = 3         if x > 400

bys id: egen ym = mean(y)  // Generate cluster means for Y and X
bys id: egen xm = mean(x)

twoway (scatter y x if id == 1, msymbol(oh)) ///
       (scatter y x if id == 2, msymbol(dh)) ///
       (scatter y x if id == 3, msymbol(sh)) ///
       (scatter ym xm,          msymbol(S)) ///
       (lfit ym xm,             lpattern(solid)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
      , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/  legend(off) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case A:" "Regional correlation," "individual orthogonality") ///
        xlabel(none) ylabel(none) name(casea, replace)

// Case B
clear
set seed 1

set obs 4 // Generate four clusters

generate id = _n
generate u_i = 300 if id == 1  // Generate cluster-specific intercepts
replace u_i = 0 if id == 2
replace u_i = 0 if id == 3
replace u_i = -300 if id == 4

expand 100  // Generate units in clusters

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0
drop if id == 2 & x >= 300
drop if id == 3 & x  < 300

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter ym xm, msymbol(S)) ///
       (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter y  x  if id == 4, msymbol(th)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit y x if id == 4,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
      , legend(order(2 "A" 3 "B" 4 "C" 5 "D" 1 "Means") size(small) title(Region:, size(small)) row(1)) ///
        xtitle("Component 1") ytitle("Component 2") ///
        title("Case B:" "Regional orthogonality," "individual correlation") ///
        xlabel(none) ylabel(none) name(caseb, replace)


// Case C
clear
set seed 1
set obs 3

generate  id = _n
generate u_i = 300 if id == 1
replace u_i = 0 if id == 2
replace u_i = -300 if id == 3

expand 150

bysort id: generate x = -100 + (600 + 100) * runiform()
generate e_ij = rnormal(0,70)                // Generate level-1 error term

generate y = x + u_i + e_ij
drop if y < 0 | x < 0

bys id: egen ym = mean(y)
bys id: egen xm = mean(x)

twoway (scatter y  x  if id == 1, msymbol(oh)) ///
       (scatter y  x  if id == 2, msymbol(dh)) ///
       (scatter y  x  if id == 3, msymbol(sh)) ///
       (scatter ym xm, msymbol(S)) ///
       (lfit y x if id == 1,    lpattern(dash)) ///
       (lfit y x if id == 2,    lpattern(dash)) ///
       (lfit y x if id == 3,    lpattern(dash)) ///
       (lfit ym xm, lpattern(solid)) ///
       , /*legend(order(1 "A" 2 "B" 3 "C") title(Region))*/ legend(off) ///
       xtitle("Component 1") ytitle("Component 2") ///
       title("Case C:" "Negative regional correlation," "positive individual correlation") ///
       xlabel(none) ylabel(none) ///
       name(casec, replace)

grc1leg casea caseb casec, legendfrom(caseb) span pos(6) row(1) name(combined, replace)