Anirban Basu

The University of Chic

The 

HERC CyberSe

June 16
Review parametric/semi parametric/semi-parametri

parametri

statistical models for health care 

expenditures

Intent to study how mean costs resp

to shifts in covariate levels to shifts in covariate levels

Will assume “selection of observab

throughout

Main focus is on bias, will look at so

efficiency issues too efficiency issues too.

Slide 

1
Slide 2
 A substantial fraction of people with no health 

substantial fraction of people with no health 

costs

 Handful of people with enormous yearly expe Often top 1% accounts for 25% of the total m

sometimes the top 10% accounts for half of a

 Any attempt to study how the mean of such a 

distribution responds to a covariate would redistribution responds to a covariate would re

careful attention to how different parts of the distribution moves in response to the change

of that covariate 

Slide 

3
 Classes of models: models:

• Single equation models• Multiple equation models

x Two-part modelsx Multi-part modelsx Conditional density estimation

• Mi t x ure mod l e s• Propensity scores• Doubly robust models 

Slide 

4
 Single-equation models: relate mean 

costs to covariates using a single funcform

• Linear models• Transformation models

x Log-OLSx Box-Cox

• Generalized Linear models

 Application to datasets of ~1000 to 10 Discuss goodness-of-fit tests

Slide 

5
Y = Costs = Costs

X = Covariates = (XT, X-T)XT = A binary treatment variable = A binary treatment variableβ = regression coefficients μ = E(Y); = E(Y); μ(x) = E(Y|X=x) (x) = E(Y|X=x)Marginal effect: dμ(X)/dXIncremental effect: Incremental effect: μ(X=1) - μ(X=0)

Slide 

6
Compare Compare the 2-year cost of care un

cost of care un

aggressive treatments (radiation or 

surgery) versus watchful waiting surgery) versus watchful waiting

following diagnosis of prostate can

U SEER se SEER-M di li k d d t b Medicare linked databasData: diagnosis between 1995-2002Follow-up through 2004.No censoring issues No censoring issuesIgnore issues with death

Slide 

7
Main parameters parameters of interest: interest:

• μ(x) = E(y|x) y = Expenditures• Incremental effect of aggressive treatmen

WW on total expenditures (ΔT)

• Avg. Costs for TRT – Avg. Costs for WW

ΔT T = = −= E ( | 1) X X { μ μ ( 1, ) ( 0

Xx X

T T 

• The interpretation of regression coefficie

{ μ μ 

−

Δ ( | 1) = − 

( 1, ) ( 0

T T E X T T X Xx X

T T 

p g 

of secondary concern.

Slide 

8
Y=Xβ + ε

Density

E(Y|Xj = x1) E(Y|Xj = x2) Y

Slide 

9
Y=Xβ + ε

Assumes parallel shift for the entire 

distribution across levels of any 

 one 

and for all X’s

Biased

• the non-linearity in the response, • the instability caused by skewness and k

Inefficient

• common failure to deal with the 

heteroscedasticity (variance increases wmean)

Slide 10
Ln(Y) = Xβ + ε

Density

E(Y|Xj = x1) E(Y|Xj = x2) Y

Slide 11
Ln(Y) = Xβ + ε

Assumes proportional shift for the e

distribution across levels of any one distribution across levels of any one 

and for all X’s

Advantages:

• Usually overcomes skewness issues• Allows for additive effects in the log Allows for additive effects in the log-scal

scal

translates to non-additive effects in the ra

• Reduces problems with heteroscedastici

kurtosis

Slide 12
Scale of estimation of estimation ≠ scale of in

of in

• If modeling $, results in log-$

 Model for geometric mean, not arith

mean

• Models E{log(y)} and not log(E{y})

Slide 13
Ln(Y) = Xβ + ε

Disadvantages:

• Creates a model for log Creates a model for log-$ not $ $, not $.• E(lnY) ≠ ln(E(Y))• Require retransformation, E(Y) Require retransformation, E(Y) = exp(Xβ 

= exp(Xβ

• Duan’s smearing estimator provides an e

for s

• s-hat = Mean{exp(lnY – X β-hat)} 

Slide 14
Ln(Y) = Xβ + ε

Disadvantages:

• E(lnY) = Xβ• lnE(Y) = Xβ + ln{s}• If s = s(X), 

dl En (Y)/dX = β + dln{ (X)}/ {s(X)}/dX

• So, common interpretation of β no longer • If ε~Normal(0, σ2), ln(s) = 0.5 σ2• With log-scale heteroscedasticity, 

Bias = 0.5dσ2(X)/dX

Slide 15
% change in costs between Free & 95%:

From log-OLS regression: 6.349 – 5.998 = 0.351 (or 35.1%) 

M

Truth (assuming log-normality) = 6.349 – 5.998 + 0.5*(2.083 – 

= 0.221 (or 22.1%)

Slide 16
 Plan specific smearing estimate can add

specific smearing estimate can add

 Difficult to come up with appropriate sm

estimator when there are multiple X’s thaheteroscedasticity

 Park Test – test to determine heterosced

log-scale

• Use squared log-scale residuals as an estimator • Study p if X’s predict this estimated variance• reg lny x, robust• predict r, res • gen r2= r^2• glm r2 x, link(log) family(gamma) rob

Slide 17
Slide 18
MODELLING HEALTHCARE EXPENDITURES 

---------------------------------------------------------. use week2_data1, clear . summ totexp,d  totexp --------------------------------------------------------- Percentiles Smallest  1% 78.32 2.78  5% 907.92 3 10% 2702.61 3.63 Obs 25% 8912.18 4.18 Sum of Wgt. 50% 14536.54 Mean 2026 Largest Std. Dev. 229375% 23375.21 279487.5 90% 41139.5 307933.6 Variance 5.2695% 59268.49 316789.1 Skewness 4.5399% 112434.3 393908.4 Kurtosis 39.7Slide 19
. kdensity totexp, saving(grp1, replace) 

0 .00001 .00002 .00003 .00004

Density

0 100000 200000 300000 400000

totexp

kernel = epanechnikov, bandwidth = 1636.38

Kernel density estimate

Slide 20
. bysort aggr_trt: summ totexp, d 

-> aggr_trt = 0  totexp --------------------------------------------------------- Percentiles Smallest  1% 50.8 3.63  5% 648.2 4.18 10% 1423.23 4.57 Obs 25% 4738.09 5 Sum of Wgt. 50% 12317.91 Mean 1920 Largest Std. Dev. 241075% 23087.59 166289 90% 43853.15 204529.9 Variance 5.8195% 65831.95 236276.2 Skewness 3.5299% 114372.3 316789.1 Kurtosis 24.2---------------------------------------------------------Slide 21
-> aggr_trt = 1 

 totexp --------------------------------------------------------- Percentiles Smallest  1% 99.18 2.78  5% 1172.94 3 10% 5141.575 4.37 Obs 25% 10071 4.55 Sum of Wgt. 50% 15343.07 Mean 2074 Largest Std. Dev. 223775% 23522.77 263867.3 90% 39927.79 279487.5 Variance 5.0195% 57092.22 307933.6 Skewness 5.1199% 110871.9 393908.4 Kurtosis 48.8. di (20745.77 - 19201.76)/ 19201.76 .08040982 Slide 22
. /* OLS */ 

. reg totexp aggr_trt, robust Linear regression Number of obs =  F( 1, 7127) =  Prob > F =  R-squared =  Root MSE = ----------------------------------------------------------------- | Robust  totexp | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- aggr_trt | 1544.008 603.106 2.56 0.010 361.7408  _cons | 19201.76 511.6135 37.53 0.000 18198.85 -----------------------------------------------------------------Slide 23
. reg totexp aggr_trt agec white black single marr

mod stage1 stage2 lnc, robust ----------------------------------------------------------------- | Robust  totexp | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- aggr_trt | 3490.596 636.6122 5.48 0.000 2242.646  agec | 205.2538 52.4853 3.91 0.000 102.367  white | 1269.85 1027.947 1.24 0.217 -745.2315  black | 3706.038 1426.357 2.60 0.009 909.9544  single | 946.6497 1241.496 0.76 0.446 -1487.051  married | -1048.995 696.5203 -1.51 0.132 -2414.382  well | -4717.357 1261.522 -3.74 0.000 -7190.314  mod | -2995.753 745.5633 -4.02 0.000 -4457.279  stage1 | -304.7167 602.6593 -0.51 0.613 -1486.108  stage2 | -1803.271 804.8815 -2.24 0.025 -3381.078  lnc | 7871.044 502.2525 15.67 0.000 6886.48  _cons | 15268.15 1349.965 11.31 0.000 12621.82 ------------------------------------------------------------------

. . 

Slide 24
. /* Log-OLS */ 

. capture drop lny . gen lny = ln(totexp) . reg lny aggr_trt, robust Linear regression Number of obs  F( 1, 7127)  Prob > F  R-squared  Root MSE ----------------------------------------------------------------- | Robust  lny | Coef. Std. Err. t P>|t| [95% Con-------------+--------------------------------------------------- aggr_trt | .3431855 .035495 9.67 0.000 .2736047  _cons | 9.147605 .0310711 294.41 0.000 9.086697 -----------------------------------------------------------------. predict xb, xb . gen mu = exp(xb) . summ mu totexp  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 7129 12040.35 1780.561 9391.919  totexp | 7129 20265.17 22937.29 2.78 Slide 25
. **** Overall smearing 

. drop mu . gen smr = exp(lny - xb) . summ smr  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- smr | 7129 1.715782 2.016464 .00021 . gen smear = r(mean) . gen mu = exp(xb)*smear . summ mu totexp  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 7129 20658.62 3055.054 16114.49  totexp | 7129 20265.17 22937.29 2.78 Slide 26
. bysort aggr_trt: summ mu totexp 

-------------------------------------------------------------> aggr_trt = 0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 2219 16114.49 0 16114.49  totexp | 2219 19201.76 24102.25 3.63 -------------------------------------------------------------> aggr_trt = 1  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 4910 22712.27 0 22712.27  totexp | 4910 20745.77 22376.88 2.78 /* Notice difference in log scale Variance */ . bysort aggr_trt: summ lny -------------------------------------------------------------> aggr_trt = 0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- lny | 2219 9.147605 1.463769 1.289233 -------------------------------------------------------------> aggr_trt = 1 -------------+---------------------------------------------- lny | 4910 9.490791 1.202411 1.022451 Slide 27
. gr twoway (kdensity lny if aggr_trt==0) (kdensit

aggr_trt==1), saving 

. 

0 .2 .4 .6

kdensity lny

0

5 10 15

x

kdensity lny kdensity lny

Slide 28
. **** Trt specific smearing 

. drop mu . summ smr if aggr_trt ==0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- smr | 2219 2.044498 2.566275 .0003865 . gen smear0=r(mean) . summ smr if aggr_trt ==1  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- smr | 4910 1.567224 1.690446 .00021 . gen smear1=r(mean) . gen mu =exp(xb)*smear0 if aggr_trt ==0 . replace mu =exp(xb)*smear1 if aggr_trt ==1 . summ mu totexp  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 7129 20265.17 714.9419 19201.76  totexp | 7129 20265.17 22937.29 2.78 Slide 29
. bysort aggr_trt: summ mu totexp 

------------------------------------------------------------aggr_trt = 0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 2219 19201.76 0 19201.76  totexp | 2219 19201.76 24102.25 3.63 -------------------------------------------------------------> aggr_trt = 1  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 4910 20745.77 0 20745.77  totexp | 4910 20745.77 22376.88 2.78 Slide 30
. **** ln-OLS regression 

. drop mu xb . reg lny aggr_trt agec white black single married stage1 stage2 lnc, robust -------------------------------------------------------------- | Robust  lny | Coef. Std. Err. t P>|t| [95% Conf. -------------+------------------------------------------------ aggr_trt | .4609178 .0372838 12.36 0.000 .3878303  agec | .0166473 .0027592 6.03 0.000 .0112385  white | .2875066 .0703636 4.09 0.000 .1495731  black | .239347 .0876897 2.73 0.006 .0674492  single | .0702465 .066109 1.06 0.288 -.0593468  married | -.0016466 .0381762 -0.04 0.966 -.0764833  well | -.3312179 .0710317 -4.66 0.000 -.4704613  mod | -.1428405 .0378877 -3.77 0.000 -.2171118  stage1 | -.0147308 .0337002 -0.44 0.662 -.0807933  stage2 | -.0830997 .0480434 -1.73 0.084 -.177279  lnc | .5280045 .0233166 22.65 0.000 .4822971  _cons | 8.627953 .0881616 97.87 0.000 8.455131 --------------------------------------------------------------Slide 31
. ** Park test - detect heteroscedasticity in log-

. predict xb, xb . gen res2 = (lny- xb)^2 . glm res2 aggr_trt agec white black single marrmod stage1 stage2 lnc, family(gamma) link(log) rob-------------------------------------------------------------- | Robust  res2 | Coef. Std. Err. z P>|z| [95% Con-------------+------------------------------------------------ aggr_trt | -.5202051 .0652508 -7.97 0.000 -.6480943  agec | -.0253315 .0054215 -4.67 0.000 -.0359575  white | -.4874537 .1051322 -4.64 0.000 -.6935091  black | -.1256915 .1388546 -0.91 0.365 -.3978414  single | -.0352572 .1358728 -0.26 0.795 -.3015629  married | -.1367799 .073717 -1.86 0.064 -.2812624  well | .100519 .1393567 0.72 0.471 -.172615  mod | -.0869734 .0785356 -1.11 0.268 -.2409005  stage1 | -.057802 .0720468 -0.80 0.422 -.1990111  stage2 | -.0547015 .1094026 -0.50 0.617 -.2691267  lnc | -.5347174 .046929 -11.39 0.000 -.6266966  _cons | 1.597288 .1376163 11.61 0.000 1.327565 --------------------------------------------------------------. test  chi2( 11) = 283.67  Prob > chi2 = 0.0000 Slide 32
( ) λ Y1 X λ βε +

 λ Æ0, Box-Cox Æ log-OLS

( ) λ Y1 X − = λ βε +

, g

 If λ is known apriori, can use Duan’s sme

estimator

 If λ is estimated, use different estimator

λ σ λ βλ ⎧ − ⎫ =+ + ⎨ ⎬ 2 1/ (1 ) E(Y|X) (X 1) 1

 All concerns about retransformation per

βλ βλ =+ + ⎨ ⎬ ⎩ ⎭ + 2 (Y|X) (X 1) 1 2(1 X ) E

p

Slide 33
 g{E(Y|X)} g{E(Y|X)} = Xβ, g{ } is a link function .} is a link function

 No retransformation problem as link fun

applies to E(Y) and NOT Y. applies to E(Y) and NOTY.

 Estimation of GLM 

• Parametric – FIML• Semi-parametric – Quasi-likelihood

 FIML:

• Gamma, Exponential, Inverse Gaussian, Pois

Neg. Binomial

• glm y x, link(log) family(gamma) 

Slide 34
 Quasi-likelihood likelihood approach approach

1. Relate μ(X) to the linear predictor Xβ using a 

function. g(μ(X)) = Xβ

2. V(Y) ~ h(μ); e.g. V(Y) = фμ2 => Gamma Varia3. Does not make any distributional assumptio

1 t d2 d t 1st and 2nd moments

4. Estimate coefficients by solving quasi-score 

∑ ∑ 1

N N

i

• glm y x, link(log) fam(gamma) robu

β μ μβ −

= =

∑ ∑= − ∂∂ = 1

1 1

( ) ( )0 j i i ii i j

i i

G yV

• glm y x, link(log) fam(gamma) robu

Slide 35
 Example: Example: Gamma with log link Gamma with log link

• ln(μ) = Xβ or μ = exp(Xβ)• V= θ1 μ2• If mean and variance functions are correctly 

β1 provides a consistent estimator of ∂ln(E{y

• If the empirical distribution of Y is truly Gam

FIML will generate the MVUE estimator for β

• If the empirical distribution of Y is deviates f

If the empirical distribution ofY is deviates f

Gamma, but mean and variance functions arappropriately specified, quasi-likelihood wil

robust estimator.

Slide 36
Slide 37
. /* Generalized Linear Models */ 

. glm totexp aggr_trt, link(log) family(gamma) robIteration 3: log pseudolikelihood = -77820.337 Generalized linear models No. of obs = Optimization : ML Residual df =  Scale parameter = Deviance = 7586.035219 1/df) Deviance = Pearson = 9205.852281 1/df) Pearson = Variance function: V(u) = u^2 [GammaLink function : g(u) = ln(u) [Log]  AIC = Log pseudolikelihood = -77820.3371 BIC = -------------------------------------------------------------- | Robust  totexp | Coef. Std. Err. z P>|z| [95% Conf. -------------+------------------------------------------------ aggr_trt | .0773403 .0307692 2.51 0.012 .0170338  _cons | 9.862757 .0266422 370.19 0.000 9.810539 --------------------------------------------------------------Slide 38
. predict mu, mu 

. summ mu totexp  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 7129 20265.17 714.9419 19201.76  totexp | 7129 20265.17 22937.29 2.78 . bysort aggr_trt: summ mu totexp -------------------------------------------------------------> aggr_trt = 0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 2219 19201.76 0 19201.76  totexp | 2219 19201.76 24102.25 3.63 -------------------------------------------------------------> aggr_trt = 1  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu | 4910 20745.77 0 20745.77  totexp | 4910 20745.77 22376.88 2.78 Slide 39
. ** GAMMA FIML 

. glm totexp aggr_trt agec white black single mawell mod stage1 stage2 lnc, link(log) family(gammaIteration 3: log pseudolikelihood = -77563.068 ----------------------------------------------------------------- | Robust  totexp | Coef. Std. Err. z P>|z| [95% Conf-------------+--------------------------------------------------- aggr_trt | .1794719 .0320742 5.60 0.000 .1166076  agec | .010508 .0025216 4.17 0.000 .0055658  white | .0854313 .0532399 1.60 0.109 -.0189169  black | .1795773 .067756 2.65 0.008 .0467779  single | .0333528 .0559583 0.60 0.551 -.0763235  married | -.0542054 .0347019 -1.56 0.118 -.1222198  well | -.2484179 .0642366 -3.87 0.000 -.3743194  mod | -.1447476 .0333845 -4.34 0.000 -.21018 - stage1 | -.0276048 .0291643 -0.95 0.344 -.0847657  stage2 | -.091044 .0429937 -2.12 0.034 -.1753101  lnc | .3698705 .0214508 17.24 0.000 .3278276  _cons | 9.628452 .0680018 141.59 0.000 9.495171 ------------------------------------------------------------------ Slide 40
. ** GAMMA QMLE 

. glm totexp aggr_trt agec white black single mawell mod stage1 stage2 lnc, link(log) family(gammairls Iteration 6: deviance = 7071.498 ----------------------------------------------------------------- | Semi-Robust  totexp | Coef. Std. Err. z P>|z| [95% Conf-------------+--------------------------------------------------- aggr_trt | .1794717 .0327351 5.48 0.000 .1153121  agec | .010508 .0025526 4.12 0.000 .005505  white | .0854311 .0535442 1.60 0.111 -.0195136  black | .1795774 .0681181 2.64 0.008 .0460683  single | .033353 .0559614 0.60 0.551 -.0763292  married | -.0542054 .0349642 -1.55 0.121 -.1227339  well | -.2484182 .0648256 -3.83 0.000 -.375474  mod | -.1447475 .0337717 -4.29 0.000 -.2109388  stage1 | -.0276047 .0292326 -0.94 0.345 -.0848996  stage2 | -.0910438 .0436959 -2.08 0.037 -.1766862  lnc | .3698705 .0216682 17.07 0.000 .3274016  _cons | 9.628452 .0676074 142.42 0.000 9.495944 -----------------------------------------------------------------Slide 41
In order to calculate effect in the ad

order to calculate effect in the ad

scale

Estimate parameters in Estimate parameters in

exp(β0 + β1*TRT + β2*X)

Predict after turning on & off the TR

Predict after turning on & off the TR

mu1 = exp(β0 + β1*1 + β2*X)mu0 = exp( mu0 = exp(β0 + β1 0* + β2 X* )Summ mu1 , mu0

Slide 42
. ** THE METHOD OF RECYCLED PREDICTIONS - TO ESTIMATE EFF

ADDITIVE-SCALE . gen trtold=aggr_trt . replace aggr_trt =1 (2219 real changes made) . predict mu1, mu . replace aggr_trt=0 (7129 real changes made) . predict mu0, mu . replace aggr_trt = trtold (4910 real changes made) . summ mu1 mu0  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu1 | 7129 21447.39 6183.882 12332.27  mu0 | 7129 17923.83 5167.942 10306.23 Slide 43
 Limitations

• Inefficient when empirical distribution of Y h

(kurtotic)

• Log Link is the most common link function us

assumption that log(μ) = Xβ may be violated 

• A i t li k i f ti t k Appropriate link or variance functions not kn

priori

• Various tests exist, not always sufficient to ov

Various tests exist, not always sufficient to ov

these problems

Slide 44
 Proposed Proposed by Manning & Mullahy (200

by Manning & Mullahy (200

• Var(Y|X) = θ1μθ2• For Gamma we assume For Gamma, we assume θ1=1, θ2=2• Test: ln(Var(Y|X)) = ln(θ1) + θ2*ln(μ)• Use Residual ^2 as an estimate for vari

Use Residual ^2 as an estimate for vari

•glm r2 lnmu, link(log) family

robust

Slide 45
 /* Modified Park test */ 

. predict mu, mu . gen lnmu = ln(mu) . gen r2 = (totexp - mu)^2 . glm r2 lnmu, link(log) family(gamma) robust ----------------------------------------------------------------- | Robust  r2 | Coef. Std. Err. z P>|z| [95% Conf-------------+--------------------------------------------------- lnmu | 1.468981 .2492923 5.89 0.000 .9803774  _cons | 5.419501 2.497738 2.17 0.030 .5240253 -----------------------------------------------------------------. test lnmu=1  chi2( 1) = 3.54  Prob > chi2 = 0.0599 . test lnmu=2  chi2( 1) = 4.54  Prob > chi2 = 0.0332 . test lnmu=3  chi2( 1) = 37.72  Prob > chi2 = 0.0000 . drop mu lnmu r2 Slide 46
 Generalized Generalized Gamma FIML Gamma FIML

 Three parameter distribution

• Mean is a function of all three parameters Mean is a function of all three parameters• Only of them is related to Xβ, implicit log-link• Special cases are Gamma, Exponential, Weibull, 

normal

• Except for log-normal case, no retransformation • stset y stset y• gengam2 x

Slide 47
. /* Generalized Gamma Regression */ 

. stset totexp . ****** Download -gengam2- command from http://home.uchicago.edu/~abasu/index_files/Pagec317.htm 

in C:\ado\Personal\ . 

. gengam2 aggr_trt agec white black single marrimod stage1 stage2 lnc, robust ----------------------------------------------------------------- | Robust  _t | Coef. Std. Err. z P>|z| [95% Conf-------------+--------------------------------------------------- aggr_trt | .193495 .0333561 5.80 0.000 .1281183  agec | .0104401 .002439 4.28 0.000 .0056598  white | .0913009 .0521438 1.75 0.080 -.0108992  black | .1798185 .0662471 2.71 0.007 .0499766  single | .0355918 .0546193 0.65 0.515 -.07146  married | -.0508413 .0332982 -1.53 0.127 -.1161046  well | -.2535959 .0625437 -4.05 0.000 -.3761792  mod | -.1453943 .0321754 -4.52 0.000 -.2084568  stage1 | -.0261988 .0281348 -0.93 0.352 -.0813419  stage2 | -.0900821 .0410532 -2.19 0.028 -.1705449  lnc | .3746245 .0203089 18.45 0.000 .3348198  _cons | 9.578895 .0675292 141.85 0.000 9.44654 -------------+---------------------------------------------------Slide 48
 /ln_sig | -.053114 .0147078 -3.61 0.000 -.0819408 

 /kappa | .8766324 .0314132 27.91 0.000 .8150637 -------------+--------------------------------------------------- sigma | .9482719 .013947 .9213265 -----------------------------------------------------------------Tests for identifying distributions -----------------------------------------------------------------Distributions | chi2 df Prob>c---------------------------------|-------------------------------Std. Gamma (kappa = sigma) | 3.27 1 0Log Normal (kappa = 0) | 778.77 1 0Weibull (kappa = 1) | 15.42 1 0Exponential (kappa = sigma = 1) | 50.89 2 0-----------------------------------------------------------------New Variables for Schoenfield Residuals created: _h*, _c* . drop smr smear smear0 smear1 res2 _* Slide 49
Pearson Pearson Correlation Correlation

• Correlation between raw-scale predictio

residuals

Pregibon’s Link Test / Reset Test

• Run same model with Xβ and Xβ^2 as cov

Hosmer-Lemeshow Test

• Plot mean residuals across deciles of Xβ.• Perform joint test that all means are zero

Copa’s Test

Slide 50
. /* GOODNESS OF FIT TESTS */ 

. ** For OLS

. reg totexp aggr_trt agec white black single marmod stage1 stage2 lnc, robust . predict xb, xb . gen mu_ols= xb . gen res_ols = totexp - mu_ols . ** Pearson Corr . pwcorr res_ols mu_ols, sig 

 | res_ols mu_ols -------------+------------------  res_ols | 1.0000  |  mu_ols | -0.0000 1.0000  | 1.0000 . ** Link Test . gen xb2=xb^2 . reg totexp xb xb2, robust  | Robust  totexp | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xb | .2645214 .3449794 0.77 0.443 -.4117406  xb2 | .0000171 8.43e-06 2.03 0.043 5.46e-07  _cons | 7353.987 3357.482 2.19 0.029 772.3251 -----------------------------------------------------------------Slide 51
. ** H-L test 

. xtile xbtile =xb, nq(10) . qui tab xbtile, gen(xbt) . reg res_ols xbt1 xbt2 xbt3 xbt4 xbt5 xbtxbt8 xbt9 xbt10, nocons robust ----------------------------------------------------------------- | Robust  res_ols | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xbt1 | 857.5334 614.4044 1.40 0.163 -346.8819  xbt2 | 681.1362 674.0399 1.01 0.312 -640.1824  xbt3 | -525.1787 642.1524 -0.82 0.413 -1783.988  xbt4 | -142.7663 699.025 -0.20 0.838 -1513.063  xbt5 | 938.962 933.5517 1.01 0.315 -891.0768  xbt6 | -638.9313 716.1985 -0.89 0.372 -2042.893  xbt7 | -1635.909 696.5204 -2.35 0.019 -3001.296  xbt8 | -2250.963 710.2455 -3.17 0.002 -3643.255  xbt9 | 707.1719 1138.647 0.62 0.535 -1524.915  xbt10 | 2007.334 1231.866 1.63 0.103 -407.489 -----------------------------------------------------------------.test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt9  F( 10, 7119) = 2.41  Prob > F = 0.0074  . drop xb xb2 xbt* xbtile Slide 52
. ** For lnOLS

. reg lny aggr_trt agec white black single marrimod stage1 stage2 lnc, robust . predict xb, xb . gen smr = exp(lny- xb) . qui summ smr . gen smear = r(mean) . gen mu_lols= exp(xb)*smear . gen res_lols = totexp - mu_lols . ** Pearson Corr . pwcorr res_lols mu_lols, sig  | res_lols mu_lols -------------+------------------  res_lols | 1.0000  |  mu_lols | -0.1776 1.0000  | 0.0000 . ** Link Test . gen xb2=xb^2 . reg lny xb xb2, robust  | Robust  lny | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xb | 7.586948 1.443073 5.26 0.000 4.758097  xb2 | -.3492898 .0757058 -4.61 0.000 -.4976956 Slide 53
 _cons | -30.9968 6.869053 -4.51 0.000 -44.46218 

-----------------------------------------------------------------

. ** H-L test 

. xtile xbtile =xb, nq(10) . qui tab xbtile, gen(xbt) . reg res_lols xbt1 xbt2 xbt3 xbt4 xbt5 xbxbt8 xbt9 xbt10, nocons robust 

 | Robust  res_lols | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xbt1 | 2479.074 635.594 3.90 0.000 1233.121  xbt2 | 2585.675 755.5887 3.42 0.001 1104.497  xbt3 | 608.5663 694.4084 0.88 0.381 -752.6807  xbt4 | 573.4585 733.1452 0.78 0.434 -863.724  xbt5 | 1438.745 775.6945 1.85 0.064 -81.84715  xbt6 | -215.099 826.1292 -0.26 0.795 -1834.558  xbt7 | -2386.189 757.4182 -3.15 0.002 -3870.954  xbt8 | -3506.212 839.8621 -4.17 0.000 -5152.592  xbt9 | -3966.24 1005.915 -3.94 0.000 -5938.133  xbt10 | -10547.35 1207.501 -8.73 0.000 -12914.42 -----------------------------------------------------------------. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt F( 10, 7119) = 15.10  Prob > F = 0.0000 . drop xb xb2 xbt* xbtile Slide 54
. ** For log-GLM 

. glm totexp aggr_trt agec white black single mawell mod stage1 stage2 lnc, link(log) family(gamma. predict xb, xb . gen mu_glm = exp(xb) . gen res_glm = totexp - mu_glm . ** Pearson Corr . pwcorr res_glm mu_glm, sig  | res_glm mu_glm -------------+------------------  res_glm | 1.0000  |  mu_glm | -0.0062 1.0000  | 0.6022 . ** Link Test . gen xb2=xb^2 . glm totexp xb xb2, link(log) family(gamm----------------------------------------------------------------- | Robust  totexp | Coef. Std. Err. z P>|z| [95% Conf-------------+--------------------------------------------------- xb | 2.601821 3.09824 0.84 0.401 -3.470618  xb2 | -.0806185 .1555559 -0.52 0.604 -.3855025  _cons | -7.950714 15.41964 -0.52 0.606 -38.17266 -----------------------------------------------------------------Slide 55
. ** H-L test 

. xtile xbtile =xb, nq(10) . qui tab xbtile, gen(xbt) . reg res_glm xbt1 xbt2 xbt3 xbt4 xbt5 xbtxbt8 xbt9 xbt10, nocons robust ----------------------------------------------------------------- | Robust  res_glm | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xbt1 | -781.7378 554.103 -1.41 0.158 -1867.944  xbt2 | 473.2616 736.601 0.64 0.521 -970.6954  xbt3 | -323.851 626.755 -0.52 0.605 -1552.477  xbt4 | 49.60152 702.7822 0.07 0.944 -1328.061  xbt5 | 1288.524 797.3756 1.62 0.106 -274.5689  xbt6 | 207.3119 861.6822 0.24 0.810 -1481.841  xbt7 | -558.6503 690.9381 -0.81 0.419 -1913.094  xbt8 | -1955.442 718.1396 -2.72 0.006 -3363.209  xbt9 | 1364.908 1148.222 1.19 0.235 -885.9489  xbt10 | 156.9733 1227.114 0.13 0.898 -2248.535 -----------------------------------------------------------------. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt F( 10, 7119) = 1.48  Prob > F = 0.1382 . drop xb xb2 xbt* xbtile Slide 56
Flexible Flexible Mean Model Mean Model

( ) 1 X if 0

λ μ

β

λ

λ

−

= ≠

log( ) X if 0

β

μ βλ = =

Flexible Variance Model

V(Y) = 

фμθθ =1 => Poisson 

θ =2 => Gamma θ

=3 IG => 

Inv. 

Gauss

Slide 57
 Goal: estimateestimate γ = (βT λ ф θ)  Goal: T estimateestimate γ = (β , λ, ф, θ)  EEE estimator estimates:

Mean model parameters (β, λ) andVar. mo e p ee s( del parameters (ф, θ)

simultaneously using estimating equ

 Extended Estimating Equations (EEE) 

Slide 58
DEL

i 1 G ( )V ( ) − ∂ ∂β

N MOD

i 1

i ii i j

j G (

y )V ( )

β = − μ ∂μ ∂β

i 1 ( )( )

MEAN L

i 1

G i ii i

(

y

)

V ( )

−

λ = −

μ

∂

μ ∂λ

MODEL

i 2 2 G (y ) V V ( V 

i i ii i

−

φ = − ⎡ ⎤ μ − ∂∂

φ ⎣ ⎦

VAR. M

i 2 2 G (y ) V V ( V 

i i ii i

−

θ = − ⎡ ⎤ μ − ∂ ∂θ ⎣ ⎦

V

Slide 59
Latest Stata code available from Stata code available from

http://home.uchicago.edu/~abasu

. summarize totcost meanonly summarize totcost, meanonly. global scale=r(mean)

. generate y=totcost/$scale

. pgl tt m yr x

. pglm y trt x [pw=marsupwt]. pglmpredict y_hat, mu scale($s

. pglmpredict ietrt, ie(trt) sca. pglmpredict me_x, me(x) scale(

Slide 60
. /* EEE - EXTENDED ESTIMATING EQUATIONS */ 

. ****** Download -pglm- command from http://home.uchicago.edu/~abasu/index_files/Pagec317.htm 

in C:\ado\Personal\ 

. qui sum totexp, meanonly . gen y = totexp/r(mean) . global scale = r(mean) . summ y  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- y | 7129 1 1.131858 .0001372 . . pglm y aggr_trt agec white black single marriemod stage1 stage2 lnc Extended GEE with Power Variance Function No of obs Optimization: Fisher's Scoring Residual df Variance: (theta1*mu^theta2) Link: (mu^lambda - 1)/lambda Std Errors: Robust -----------------------------------------------------------------Slide 61
 y | Coef. Std. Err. z P>|z| [95% Conf

-------------+---------------------------------------------------y |  aggr_trt | .176795 .0319912 5.53 0.000 .1140935  agec | .010301 .0025097 4.10 0.000 .005382  white | .0772856 .0524606 1.47 0.141 -.0255354  black | .1771791 .0671802 2.64 0.008 .0455084  single | .0369893 .0563494 0.66 0.512 -.0734535  married | -.053268 .0341707 -1.56 0.119 -.1202413  well | -.2421367 .063453 -3.82 0.000 -.3665023  mod | -.144254 .0339966 -4.24 0.000 -.2108861  stage1 | -.0237874 .0289403 -0.82 0.411 -.0805093  stage2 | -.0913624 .0418396 -2.18 0.029 -.1733666  lnc | .3719294 .0221651 16.78 0.000 .3284867  _cons | -.2769968 .0682525 -4.06 0.000 -.4107693 -------------+---------------------------------------------------lambda |  _cons | .1928841 .303307 0.64 0.525 -.4015866 -------------+---------------------------------------------------theta1 |  _cons | 1.166449 .0810133 14.40 0.000 1.007666 -------------+---------------------------------------------------theta2 |  _cons | 1.457585 .2222653 6.56 0.000 1.021953 -----------------------------------------------------------------. predict xb, xb . pglmpredict mu_pglm, mu scale($scale) Slide 62
. summ mu_pglm totexp 

 Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- mu_pglm | 7129 20265.42 5696.584 10770.24  totexp | 7129 20265.17 22937.29 2.78 . pglmpredict trteffect, ie(aggr_trt) scale($scale. summ trteffect  Variable | Obs Mean Std. Dev. Min -------------+---------------------------------------------- trteffect | 7129 3470.224 790.525 2190.088 . gen res_pglm = totexp - mu_pglm GOODNESS OF FIT

. 

. ** Pearson Corr . pwcorr res_pglm mu_pglm, sig  | res_pglm mu_pglm -------------+------------------  res_pglm | 1.0000  |  mu_pglm | -0.0013 1.0000  | 0.9158 Slide 63
. ** H-L test 

. xtile xbtile =xb, nq(10) . qui tab xbtile, gen(xbt) . reg res_pglm xbt1 xbt2 xbt3 xbt4 xbt5 xbxbt8 xbt9 xbt10, nocons robust  | Robust  res_pglm | Coef. Std. Err. t P>|t| [95% Conf-------------+--------------------------------------------------- xbt1 | -609.5808 553.8617 -1.10 0.271 -1695.314  xbt2 | 646.0331 742.0791 0.87 0.384 -808.6624  xbt3 | -471.0719 618.8141 -0.76 0.447 -1684.131  xbt4 | -84.50918 702.9003 -0.12 0.904 -1462.403  xbt5 | 1116.715 791.1745 1.41 0.158 -434.2219  xbt6 | 365.4751 872.6898 0.42 0.675 -1345.256  xbt7 | -1006.109 683.5575 -1.47 0.141 -2346.085  xbt8 | -1908.943 719.4691 -2.65 0.008 -3319.317  xbt9 | 1258.806 1145.722 1.10 0.272 -987.151  xbt10 | 691.741 1227.973 0.56 0.573 -1715.451 -----------------------------------------------------------------. test xbt1 xbt2 xbt3 xbt4 xbt5 xbt6 xbt7 xbt8 xbt F( 10, 7119) = 1.55  Prob > F = 0.1164 . drop xb xbt* xbtile Slide 64
 Empirical Empirical example example

• Cost attributable to Heart Failure (HF) post m

infarction (MI)

• Use Medstat data, extract all claims between 

and December 31, 2000 

Prior period

6 months

Post period

12 months

Index MI

HF = I(Patient who develop HF in post-period)

MI

Date

Slide 65
Exclusion Criteria Numbe

Starting cohort  No claims in pre-MI or enrollment information 1279  Patients with pre-MI HF only 1696  Patients with pre-MI MI only 806 

P ti t ith Patients with pre-MI HF d MI MI HF and MI 269

 Patients with post-MI MI 575  Post MI HF patients with intermediate coronary 

syndrome and coronary occlusion without MI syndrome and coronary occlusion without MI 18

 Patients with capitation type insurance plans 2258  Patients with VSD 16  Patients with mitral regurgitation 46  Patients with cost of index hospitalization < $500 9 Patients with LOS>100 and cost < $10,000 1  Patients with age<18 years or missing age 3 

Patients with missing plan information Patients with missing plan information 20

 Patient with index date not in (1Dec1998 – 1Dec 2000) 4633 Eligible cohort 

Patients Patients without full 1-year follow year follow-up data. up data. 3427

 Observations with missing covariate values 3 Final sample for this study 

Slide 66
 With Heart Failure

Without Heart Failure

ensity

Without Heart Failure

ensity

De

De

0 20000 40000 60000 80000 100000

Total expenditures

-4 -2 0 

Log-scale residuals for exp

With HF: N=2259 With HF: N = 2259 Mean (sd) = $45,070 (57,366) ; Sk. = 6.0; Kurt. = 73

Sk. = 0.13; Kurt. = 3.2

Without HF: N=5169 Without HF: N = 5169Mean (sd) = $31,708 (30,110) ;Sk. = 4.4; Kurt. = 50

Sk.= -0.03; Kurt. = 3.4

Slide 67
 Final Sample: Sample: 7428 patients patients

2259 (30.4%) with HF

 Dependent variable: Total expenditure Independent variables: HF, age, age-sq

insurance type, Medicare status, proceperformed year of MI type of performed, year of MI, type of MI and

indicators for 30 categories of co-morb

 Compare: OLS, log-OLS(with HF-specif

smearing), GLM Gamma-Log link.

Slide 68
π

ˆ (π π ˆ ˆ Row Col − ) 

M ( td )* Mean (std.err.)

Estimator 

Mean (std.err.)

Log-OLS 

(Hom.) 

Log-OLS 

(Het.) 

Gamma 

(log link)

OLS 15578 (1172) 959 1913 312 

(768) -1913 (627) -312 (566) 

Log-OLS (H )

14619 ( 901) -2872 -1271 

(Hom.) (599) (574) Log-OLS (Het.) 

17491 (1159) 1601 

(239)( ) 

Gamma (log link) 

15890 (1104) 

EEE 14722 (1135)

Slide 69
15000

O

60009000

12000

)

LL

0

3000

6000

esiduals ($ G

E

9000

-6000-3000

Mean Re

-15000-12000

-9000

1 2 3 4 5 6 7 8 9 10

Deciles of Linear Predictor Basu 

Slide 70
 Results Results from the EEE Model from the EEE Model

• Mean Model (Link parameter)• = 0.201 (95%CI: 0.04, 0.37)

ˆ

λ

• Appropriate link ~ μ0.2 or μ0.25• Variance Model:V(y) = фμθ • Variance Model:V(y) = фμ• = 0.80 (95% CI: 0.72, 0.89); = 1.79 (95% C• Appropriate variance ~ μ2

ˆ

φ

ˆ

θ

 Suggest fifth or fourth-root link with G

variance.

Slide 71
π

ˆ (π π ˆ ˆ Row Col − ) 

M ( td )* Mean (std.err.)

Estimator 

Mean (std.err.)

Log-OLS 

(Hom.) 

Log-OLS 

(Het.) 

Gamma 

(log link)

OLS 15578 (1172) 959 1913 312 

(768) -1913 (627) -312 (566) 

Log-OLS (H )

14619 ( 901) -2872 -1271 

(Hom.) (599) (574) Log-OLS (Het.) 

17491 (1159) 1601 

(239)( ) 

Gamma (log link) 

15890 (1104) 

EEE 14722 (1135)

Slide 72
 Multi-part models

• Two-part: 1st part model any costs (logit/p

2nd part models costs given any

• Multi-part: 2-part models for each type of e• CDE: 1st part for multinomial categori

d 2 f n parts or levels within each c

 Mixture models

Slide 73
 Comp g arisons throu

gh simulations or 

empirical applications:

 Duan et al., JBES 1983; Mullahy, JHE 1998; Lipscomb et al

1999; g Mannin

g & Mullahy,J ; HE 2001; Deb & Trivedi,J 

HE 2

Burgess, Hunter College Mimeo 2003; Basu et al., HE 20& Zaslavsky, JHE, 2004; Basu & Rathouz, Biostatistics 2005al., HE 2005; Manning et al., JHE 2006; Basu et al., HE 200Miller HE 2009; Miller, HE 2009;

 Consistently found that no one estima

app p ro

priate for all ex

penditure data

 Theory has difficult time predicting si

covariate effects – provides almost no guidance on functional form

Slide 74
 Longitudinal Longitudinal costs – generalized generalized estimat

estimat

equations versus random effects model

 Right censoring Right censoring Survival effects versus intensity effects Causal estimation Causal estimation – instrumental variabl

instrumental variabl

 Predictions

Basu and Manning, Medic

Slide 75