Statistical Methods for Analysis of H

Costs

Xiao-Hua Andrew Zhou

azhou@u.washington.edu

HSR&D Center of Excellence, Seattle VA Medical C

Department of Biostatistics, University of Washing
Outline

• Challenges in the analysis of health care costs.• Misuse of statistical methods in analysis of costs in medical lite• Some new statistical methods for the analysis of health care co
Some challenges in analysis of costs

• A proportion of the population can usually be expected to incu

typical study period.

• Non-zero cost observations are typically skewed to the right; th

be approximated by a log-normal distribution (Diehr et al., 200

• Health care costs typically exhibit heteroscedasticity.
Focus on population means

• Although it is a common practice to use the median instead the 

measure of central location in skewed data, many applications of means.

• This is particularly common in the analyses of medical cost da

mean can be used to recover the total cost, which reflects the on health care in a given patient population.
Misuse of methods for analysis of costs in medical literatur

e

• Zhou et al. (1997a, Annals of Internal Medicine) reviewed stati

studies of medical costs published in medical journals between January, 1996 and found that at least 26% of the studies might conclusions.

• Barber and Thompson (1998, British Medical Journal) critically 

statistical methods used in analysis of health care costs in randthey found that in at least two thirds of the published papers, thregarding costs were not justified.
Some newly developed statistical methods

• One population• Two populations• Three or more populations• Regression models.
Statistical problem in one population

• Let

W1, . . . , Wn be a random sample from a skewed distributio

variance

τ

2. We are interested in point and interval estimators 
Point estimators of θ

• Let W1, . . . , Wn be a random sample from a skewed distributio

variance τ 2. We are interested in point and interval estimators 
Point estimators of θ, cont

Assuming that Y = log W ∼ N(μ, σ2) (Zhou, 1998, Stat Med), we h

estimators for θ.

• The sample mean, W ̄ .• the ML estimator:

exp(Y ̄ + m 2(m + 1)S2),

where m = n − 1, Y ̄ and S2 are the sample mean and varianc

• A uniformly minimum variance unbiased (UMVU) estimator:

exp(Y ̄ )gm(S2/2),

where

gm(t) = X

∞

r=1

1

r!

m + 2r

m

( m

m + 1

t)r Y

r

i=1

m

m + 2

• A conditionally minimal MSE estimator:

exp(Y ̄ )gm(m − 3 2m S2).
Results

• After deriving the mean square errors of the four estimators, w

following results:

• The standard estimator, the sample mean, has the largest MS• We recommend the use of the ML estimator when skewness is • Otherwise we recommend the use of the conditionally minimal • Future research problem: searching for better estimators, part

log-normal distribution does not hold.
Interval estimation of θ

Assume that Yi = log Wi ∼ N(μ, σ2). A commonly used (naive) inte

• Construct a confidence interval based on log-transformed outc

Y ̄ ± Z1−α/2 √Sn .

• Transform the interval back to the original outcome,

exp(Y ̄ ± Z1−α/2 S√n).

• Since this interval is for exp(μ), not for θ, the naive interval is b
Interval estimation of θ, based on the original data

• Large-sample central limit theory:

W ̄ ± Z1−α/2√τb2n,

where W ̄ and τb2 are the sample mean and variance of the orig

Wi’s.

• Student’s interval of θ, based on the original data:

W ̄ ± t1−α/2,n−1√τb2n,
New Interval One

• Note that log(θ) = μ + σ2/2.• The biased corrected ML estimator of log θ: Y ̄ + S2/2 with the 

S2/n + S4/(2m).

• a 1 − α level confidence interval for θ:

exp(Y ̄ + S2/2 ± Z1−α/2sS2n + S4 2(n − 1) 

• A simulation study suggests this interval has good coverage pr

Gao (1997, Stat in Med)).

• Or, the modified version (Olsson, 2005, Journal of Statistical E

exp(Y ̄ + S2/2 ± t1−α/2,n−1sS2n + S4 2(n − 1

• The modified version outperforms the original one when the sa
Generalized confidence intervals

• Krishnamoorthy and Mathew (2003) proposed an interval for

θ 

generalized confidence intervals ( Journal of statistical plannin115, 103-121)
Simulation study

• Olsson (2005) compared the performance of the naive approa

approach, the modified method with t instead of z as multiplierconfidence intervals, and the large-sample central limit method

• The large-sample method gives a consistently lower coverage 

level.

• The easy to compute new and modified confidence intervals p

modified method being better in small sizes.

• The generalized confidence interval approach also works well; 

disadvantage is that it requires a computer to simulate the sam
A real example

• We illustrate the method using data from a study of the effect o

charges following knee replacement (KR) procedures.

• This dataset consists of hospital charges for 355 obese patien

operations.

• The distribution of costs was skewed significantly toward highe

the log transformed data approximate a normal distribution.

• A formal Shapiro-Wilk test for the normality on the log-transfor

p-value of 0.25.
Results

The 95% confidence intervals for the mean of hospital charges:

• The naive method,

[$8, 839.6, $9, 363.7]

• The new method,

[$9, 326.0, $9, 893.2]
A general skewed distribution

• If the distribution of Wi is unknown, some modified t-intervals t

limited success have been proposed using an Edgeworth expastandard t-statistic.

◦ Here are some references: Hall (1992, Biometrika) Sutton 

(1995, JASA), Zhou and Gao (2000, Amer. Statist.).

• Future research: finding better non-parametric intervals
Statistical problem in one population with additional zero

s

• Let

W1, . . . , Wn be a random sample from a lognormal distrib

additional zero values. That is, if

Wi >

0, logWi has

N(μ, σ

2

)

• Construction of an interval for

θ

=

E(Wi) =

P(Wi > 0)

E(Wi 
Statistical problem in one population with additional zero

s

• Let

W1, . . . , Wn be a random sample from a lognormal distrib

additional zero values. That is, if

Wi >

0, logWi has

N(μ, σ

2

)

• Construction of an interval for

θ

=

E(Wi) =

P(Wi > 0)

E(Wi 
Interval estimation

• Owen and DeRouen (1980, Biometrics) derived a minimum va

estimator (MVUE) confidence interval for θ.

• Zhou and Tu (2000, Biometrics) have proposed a percentile-t b

based on the sufficient statistics, a biased-corrected maximum estimation using normal approximation, and an interval based log-likelihood ratio test statistic. The bootstrap and likelihood raintervals were recommended for means of lognormal data with methods perform not so well in very small sample situations.

• Tian (2005, Stat in Med) proposed an alternative generalized i

for confidence interval estimation and hypothesis testing.
Two populations without any zeros

•

Wi1: the outcome variable of the ith patient in the first sample,

the outcome variable of the jth patient in the second sample,

i 

Mj

=

E(W ij), and

σ

2

j

= V ar(W ij).

• Assume that

Xij

= logWij

∼

N(μ

j , σ

2

j

).

• The null hypothesis of interest is

H

0 : M

1

=

M2.
Some existing methods

• Student’s test on Wij ’s; it is valid when both n1 and n2 are lar• Wilcoxon non-parametric test, based on Wij ’s• Student’s test on logWij ’s• Two-sample bootstrap test.
A simple valid test

• Zhou et al (1997, Biometrics) proposed the following test statis

Z = μ

b

2 − μ

b

1 + (1/2)(S22 − S21 )

rS21 n1 + S22 n2 + (1/2)( S41 n1−1 + S42 n2−1 ) 

where μbj and S2j are the jth sample mean and variance.

• The two-sided p-value = 2Φ(− | Z |).
Simulation results

• The Z-score method has the observed type I error rate that is t

pre-set nominal level even with small sample sizes (n1 = n2 = 

• The type I error rates of the other existing four tests are all larg

level and become larger as the difference in two skewness inc

• With unequal sample sizes, the Z-score method has a greater 

other tests.
More complicated tests

• Modified signed Likelihood ratio tests for the ratio of means of 

normal distributions (Wu et al, 2002).
Two extensions

• Zhou and Tu (2000, Computational Statistics & Data Analysis) 

confidence intervals for the difference in and the ratio of the mewith additional zeros.

• Zhou et al. (2000, Stat in Med) proposed new tests for compar

care costs in a paired design study.
Three or more populations with additional zeros

• Let

W1

j , . . . , Wn

j

j be a random sample from the jth populatio

additional zeros and

Mj

=

E(Wij ), j = 1, . . . , K.

• Assume that for

Wij

>

0, Yij

= logWij

∼

N(μ

j , σ

2

j

). The nu

interest is

H

0 : M

1

=

. . .

=

MK.
New Tests

• Zhou and Tu (1999, Biometrics) proposed a computationally m

likelihood ratio test.

• Tu and Zhou (1999, Stat in Med) proposed a computationally s

test.

• Simulation results suggest (1) that the likelihood ratio test has 

rate (closest to the nominal level), and is closely followed by thfor unequal sample sizes, the likelihood ratio test has better cothan the Wald test, and (3) that when the sample sizes are largrates of the two tests are quite close to the nominal level.
Regression models for costs

• We are interested in the effect of patient-level factors (such as 

compliance and patients¡
R satisfaction with their health care p

care costs of patients.
Existing regression models for skewed data

In statistical literature, there are four common ways of modeling skew

1. The standard linear regression model with Ordinary least squa

any transformation.

2. The Cox proportional hazards model (Dudley et al, 1993).3. A parametric skewed distribution family for

ǫ, leading to genera

with an exponential family (Blough et al (1999)).

4. Transformation model so that transformed costs have a particu

distribution, e.g. normal, homoscedastic, symmetric distributioextreme skewness with more efficient estimation (Ruppert (200
Comments on OLS linear models

Advantages

1. Easy2. No retransformation problem3. Easy to compute marginal and Incremental effects

Disadvantages

1. Clear violation on normality and homoscedasticity.2. Not robust in small to medium sized data set or in large datase

observations

3. Can obtain predictions with negative costs
Comments on Cox regression models

Advantages

1. Semi=parametric model without assuming the normality assum

Disadvantages

1. The regression coefficients in the Cox proportional hazards mo

hazard ratio, it is difficult to interpret them in the context of hea
Comments on GLM

Advantages

1. No retransformation problems2. Gains in precision from estimator of the assumed model holds

3. Consistent even if the variance function is misspecified.

Disadvantages

1. Can suffer substantial precision losses if heavy-tailed (log) erro

log-scale residuals have high kurtosis (> 3) or if variance funct
Extension of GLMs

• A non-parametric GLM by Chiou and Muller (1998, JASA): g(E

V ar(W) = σ2(E(W)), where both g(.) and σ(.) are unknown 

• A semi-parametric GLM by Basu and Rathouz (2005, Biostatis

gλ(E(W)) = X′β, and V ar(W) = h(E(W); γ1, γ2), a param

mean with two unknown parameters. Here gλ(y) is a Box-Cox 

function, and h(E(W); γ1, γ2) = γ1(E(W))γ2 , a power family 

h(E(W); γ1, γ2) = γ1E(W) + γ2(E(W))2, a quadratic varian
Comparison of these two approaches

• Usually, generalized linear models and transformation models 

non-linear regression relationships between

E(W) and

X.

• Which model is correct will depend on a particular application 

have additive or multiplicative errors.

• Unlike the transformation model, with the GLM we do not have 

re-transformation bias.

• The GLM addresses skewness by the choice of a distribution f

used one being a Gamma, Poisson, or negative binomial distrithe non-linearity by the choice of its link function, commonly usor square root link.

• Estimation of GLM still only use the first two moments of cost d

information from the third and higher moments.
Advantages of transformation models over GLM

• When the expected value of

W is related to a vector of covaria

way, often a transformation of

W will simplify this relationship b

linearities or removing interactions (Ruppert, 2000). Suppose t

W

=

β

0

X

β1

1 . . . X

β

J

J

+ ǫ,

a complicated non-linear model, where

ǫ is a small random err

log(μ

+

ǫ) ≈ log(μ) + ǫ/μ for small values of

ǫ, we obtain that 

approximately the linear model

log(W) =

β

∗

0

+

β

1

X

∗

1

+ . . .

+

β

J

X

∗

J + (

β

0

X

β1

1 . . . 

where

β

∗

0

= log(β

0), and

X

∗

j

= log(Xj).
Advantages of transformation models over GLM, cont

• The transformation method can induce a particular type of dist

extreme skewness so that more efficient estimators and more can be obtained. For example, for highly skewed health care cthe data are crowded into the lower left-hand corner of the plotwhat type of regression model is appropriate between the outc

• A transformation that can induce the symmetric distribution of 

would make it much easier to see a relationship between the trand a predictor.

• The transformation method for estimating

E(W) uses informat

moments.
Re-transformation bias

• One complication in the use of a transformation model is poss

re-transformation bias.

• As we are interested in dollar scale not in transformed-dollar s

transform regression results on the transformed scale back to dollar-scale.

• If this re-transformation is not properly done, we may introduce 

called "re-transformation bias".
Existing semi-parametric transformation model

• Welsh and Zhou (2006) and Zhou and Cheng (2008) proposes 

semi-parametric transformation model to analyze skewed and variance data:

H(Y ) =

X′β

+

σ(X′

γ)ε, 

where

Y is a scalar dependent variable,

H(·) is a known incre

function,

σ(·) is the known variance function,

X is a

q

×

1 vect

explanatory variables with the first element being

1, β and

γ a

unknown parameters, and

ε is an error term that has a unknow

mean 0 and variance 1. Here the error term

ε is independent o

• In the model (0.4), we allow the effect of

X on the mean and v

different.
Software

• We developed a computer program in R.• It compiles and runs on a wide variety of UNIX platforms, Wind• Our program can run on R 2.0.0 and later releases.• The program is available from http://faculty.washington.edu/∼a
Software, cont

• Our program reads in patient data files in plain text format• Estimates regression parameters αb, β

b

, and θ

b

• Then computes both the externally weighted estimator ub∗ and 

weighted estimators ub of mean on the original scale.

• Welsh and Zhou (2005) showed that both the internally weight

externally weighted estimator have very similar bias and MSE.

• While the internally weighted estimator has slightly smaller bia

weighted estimator, the externally weighted estimator has slighthan the internally weighted estimator.

• Therefore, we give both estimators so that users can choose w
Software, cont

• In addition to these two mean estimators, our program also ou

including standard deviation, asymptotic confidence interval, abootstrap confidence interval.

• Our program can run in both interactive mode and batch mode
An example

• We illustrate the use of our computer program in a data set on 

patients from a prospective drug utilization review (DUR) study

• A goal of our analysis is to estimate the average of in-patient c

given his/her age, gender, race, and general health status as m

• Since the in-patient charges are zero for some patients, we ap

heteroscedastic regression model to our data set.
An example, cont

• Let Yi be the in-patient charge of the ith patient, and correspo

defined as follows.

• Xi1 is the age of the patient; Xi2 represents the patient’s race 

Caucasians and Xi2 = 0 for African Americans); Xi3 represen

patient (Xi3 = 1 for males and Xi3 = 0 for females); Xi4 is th

100 representing ith patient’s general health status.
An example, cont

• For i = 1, . . . , n, we model the probability of non-zero in-patien

logistic regression model,

log

P(Yi = 0 | Xi1, . . . , Xi4)P(Yi > 0 | Xi1, . . . , Xi4)

= α0 + α1Xi1 + . . 

• and we model the conditional magnitude of the positive charge

by the log-transformed, heteroscedastic linear regression mod

log Yi = β0 + β1Xi1 + . . . + β4Xi4 + exp{(θ0 + θ1Xi1 + . . . 
Input data

0 56 0 0 320 64 1 1 251952.05 68 1 1 42

0 54 1 1 50......
Input data

Please specify the name of input data file: examplechoose method for estimating parameters beta and th

1: MLE estimator, maximize log-likelihood functio2: MLE estimator, solve estimating equation

1

Do you want to give initial guess of parameters (d

zeros)? (Y/N)n alpha.0 = 0 0 0 0 0 beta.0 = 0 0 0 0 0 0

Please specify (1-r) confidence interval (default 

r =r = 0.05

Do you want to calculate bootstrap confidence interPlease specify bootstrap sample size (default 100): 

Do you want to assign seed for random number gener

Seed not assigned.
Commands

Please specify the name of input data file: example

Please choose method for estimating parameters bet

1: MLE estimator, maximize log-likelihood functio2: MLE estimator, solve estimating equation

1

Do you want to give initial guess of parameters (d

zeros)? (Y/N)n alpha.0 = 0 0 0 0 0 beta.0 = 0 0 0 0 0 0

Please specify (1-r) confidence interval (default 

r =r = 0.05

Do you want to calculate bootstrap confidence inte

Please specify bootstrap sample size (default 100)Do you want to assign seed for random number generSeed not assigned.
Commands, cont

Please input covariates values: 1: 55 1 0 50 CovarDo you want to see estimation results for another Please input covariates values: 1: 65 1 0 50 CovarDo you want to see estimation results for another 
Results

estimator of parameter alpha = 1.513531 -0.00802566

0.3740364 0.007434743std of alpha estimator = 0.7334805 0.01321834 0.310.006093177

estimator of parameter beta = 9.538691 -0.004213263 

0.02954837 0.003566381

std of beta estimator = 0.701519 0.0126878

0.3566817 0.2814012 0.005726485

estimator of parameter theta =

-0.9736211 0.05375912 -1.048779 -0.5785688 -0.01572

std of theta

estimator = 0.9676426 0.01870368 0.4402539 0.376311
Results, cont

Source data file: example.dat Number of obseNumber of covariates: 4 ---------------------------------------------------Covariate = 2 2 2 2

Externally weighted estimator:

mean = 111.6968

standard deviation = 145.8094

95% confidence interval = [0, 397.478]

95% bootstrap confidence interval = [0.7717479, 352Internally weighted estimator:

mean = 111.3817

standard deviation = 144.0542

95% confidence interval = [0, 393.7227]

95% bootstrap confidence interval = [0, 3539.867] ---------------------------------------------------
A new non-parametric transformation model

• We propose the following nonparametric transformation model 

and heteroscedastic variance data:

H(Y ) =

X′β

+

σ(X′

γ)ε, 

where

Y is a scalar dependent variable,

H(·) is a unknown inc

transformation function,

σ(·) is the known variance function,

X 

observed explanatory variables with the first element being

1, β 

of unknown parameters, and

ε is an error term that has a unkn

with mean 0 and variance 1. Here the error term

ε is independ

• In the model (0.4), we allow the effect of

X on the mean and v

different.
Identifiability Assumptions

• To make the model (0.4) identifiable, we need 

following assumptions.

• Like Horwitz (1996), we assume that there ex

that H(y0) = 0.

• Re-arrange X so that its first component X1 h

absolute continuous density conditional on X2

β1 be the corresponding coefficient. | β1 |= 1.
Some special cases

• Model (0.4) includes as special cases a large number of widely 

extensively investigated models that make stronger assumptiomade in the paper about H and F.

• Linear regression models, log-linear regression models, the Co

hazard model, and accelerated failure time models

• Transformation models in which H is specified up to a vector o

parameters (e.g.,Box and Cox, 1964; Bickel and Doksum, 198

• Transformation models in which F is specified up to a vector o

parameters, and H is nonparametric (Cheng et al.1995, Dabro1988).
Some notation

• Let {Yi, Xi, i = 1, · · · , n} be a random sample of (Y, X) that 

(0.4).

• Denote Z1 = X′β , Z2 = X′γ, Z1i = X′iβ, and Z2i = X′iγ.• Let G(·|z1, z2) be the cumulative distribution function (CDF) of 

Z1 = z1 and Z2 = z2, and p(·, ·) be the probability density fun

• Assume that H, F, and G are differentiable with all their argum• Define h(y) = dH(y)/dy, f(y) = dF(y)/dy, p(y|z1, z2) = ∂G

gj (y|z1, z2) = ∂G(y|z1, z2)/∂zj , j = 1, 2.
Estimation method

• Since under the model (0.4), Y depends on X only through th

the model (0.4) implies that

G(y|z1, z2) = F( H(y) − z1

σ(z2) ),

and we can show that p(y|z1, z2) = −g1(y | z1, z2)h(y).

• Denote g1(y, z1, z2) = g1(y|z1, z2)p(z1, z2) and

p(y, z1, z2) = p(y|z1, z2)p(z1, z2), we get

g1(y, z1, z2)h(y) = −p(y, z1, z2). 

• Hence we obtain that h(y) = − Pni=1 P p(y,Z1i,Z2i) ni=1 g1(y,Z1i,Z2i) , and

H(y) = − Z yy0 Pni=1 p(u, Z1i, Z2i) Pni=1 g1(u, Z1i, Z2i) du. 

• The expression (0.6) forms the basis for the estimator of H pro
Estimation of H(.)

• From (0.6), we see that to derive an estimator of H(.), we nee

p(z1, z2), G(y|z1, z2) and derivatives of G(y|z1, z2) when the 

are given.
Estimation of H(.), cont

• We estimate G(y|z1, z2) by the following kernel estimator:

Gn(y|z1, z2) = 1 nh1h2pn(z1, z2) Xni=1 I(Yi ≤ y)K1 „Z1i − z1 h1 where K1 and K2 be bounded and symmetric kernel functions [−1, 1] with h1 and h2 being bandwidths.

• Here pn(z1, z2) is the kernel density estimate of p(z1, z2), and 

pn(z1, z2) = 1 nh1h2 Xni=1 K1 „Z1i − z1 h1 « K2 „Z2i − 

h2 
Estimation of H(.), cont

• Since g1(y|z1, z2) = dG(y|z1, z2)/dz1, we obtain an estimato

differentiating Gn(y|z1, z2) with respect to z1,

g1n(y|z1, z2) = ∂Gn(y|z1, z2)/∂z1. 

• Although p(y|z1, z2) is the probability density function of Y con

Z1 = z1, Z2 = z2, it can not be estimated by ∂Gn(y|z1, z2)/∂Gn(y|z1, z2) is a step function of y.

• Instead, we use the following kernel density estimator for p(y|z

pn(y|z1, z2) = 1 nh1h2h0pn(z1, z2) Xni=1 K0(Yi − y h0 )K1 „Z1i − 

h

where K0 be bounded and symmetric kernel functions with the 

with h0 being a bandwidth.
Estimation of H(.), cont

By substituting (0.8), (0.9) and (0.10) into (0.6), we obtain the estima

Hn(y) = − Z yy0 Pni=1 pn(u|Z1i, Z2i)pn(Z1i, Z2i) Pni=1 g1n(u|Z1i, Z2i)pn(Z1i, Z2i) du
Estimation of β and γ

• Since E((H(Y ) − Z1)2 | X) = σ2(X′γ), we can use the follo

equations to simultaneously estimate β and γ:

X

n

i=1

(H(Yi) − X′iβ)Xi

σ2(X′iγ) = 0, 

and

X

n

i=1

 ̆(H(Yi) − X′iβ)2 − σ2(X′iγ) ̄ Xi = 0, 

when given H.

• From the equation (0.12), we obtain a closed-form estimator of 

βn = Xni=1 XiX′i σ2(X′iγ) !−1Xni=1 XiH(Yi

σ2(X′iγ
Estimation algorithm of β, γ and H(·).

1. Selection of initial values.

(a) Initial values for β and H. We can still obtain consistent e

H even we misspecify the variance function, and hence, w

reasonable starting values for β and H with estimates obthomoscedasticity model,

H(Y ) = X′β + σε. 

Under this homoscedastic model, we can estimate β by thcorrelation (MRC) method proposed by Han(1987); that is

β

e

= argmaxβWn(β), where Wn(β) = Pi=6 j{Yi > Yj}{

then we can estimate H using the proposed method with that K2( Z2i−z2 h2 ) = 1 for any z2 and i = 1, · · · , n.

(b) Initial values for γ. Given β and H, we estimate γ by the e

2. Estimation of H(.). Given β and γ, we estimate H by (0.11).3. Estimation of β and γ. Given H, we estimate β and γ by (0.144. Iteration. Repeat Steps 2 and 3 until two successive values of 

significantly.
Prediction of μ(x) = E(Y | X = x)

• For given covariates x of a patient, we are interested in predic

model (0.4), we can write

μ(x) = Z H−1(xT β + σ(xT γ)u)dF(u). 

•• We propose to estimate F by the empirical distribution F

b

of the 

residuals, ebi = H

b

(Yi)−XTi βb

σ(XTi γb) , where H,

b

β

b

and γb are the estima
Non-parametric estimator

• Therefore, replacing

H, β, γ, and

F by

Hb, βb, γb, and

Fb in (0.16

following estimator of

μ(x):

μb(x) =

1

n

Xn

i=1

Hb

−

1

 

x

′βb

+

σ(x

′

γb)

Hb(Yi) −

X′

iβb

σ(X′

i

γb)

• This estimator can be considered as an extension of Duan’s sm

(Duan, 1982) to the heteroscedastic transformation model with transformation and error distribution functions.
Simulation studies

• We conduct a simulation study to assess the finite-sample perf

proposed method.

• Unlike the existing parametric or the single-semiparametric mo

the transformation and error distribution functions is specified, method does not rely on parametric specifications for both the error distribution functions.

• Hence we expect our estimators of the untransformed scale ex

regression parameters are more robust than the ones derived parametric and single-semiparametric methods.

• We also want to know whether the added robustness is gained 

reduced efficiency.
Simulation studies, cont

• To investigate this, we compare the performance of the propos

following models:

1. the CTCD model, where the transformation and error distr

correctly specified by a parametric model, the case that sestandard,

2. the CTMD model, where the transformation is correctly sp

distribution is misspecified, and

3. the MTCD model, where the error distribution is correctly 

transformation function is misspecified. The CTCD model investigate the efficiency of the proposed method, and the models are used to investigate the robustness of the prop

4. Two existing methods, Basu’s method and Chiou’s method
First simulation study

• In our first simulation study, we consider a true transformation 

with one binary covariate, one continuous covariate, and a nontransformation function. For n = 2000 subjects, we generate c

X2 from the binomial distribution with p = 0.5 and the uniform 

respectively, the random error ε from the standard normal distroutcome follow the following transformation model:

H(Y ) = β0 + X1β1 + X2β2 + p0.4 + γX1

where H(y) = Φ−1{exp(y − 10)}, β0 = −1.8, β1 = 1.4 β2 = 

• We assess the performance of the various estimators of the re

parameters and the untransformed scale expectation in terms deviation (SD), bias, and the root of mean squared error (RMS

• In the MTCV model, the transformation function is misspecified 

H(y) = exp(y − 10).
Simulation studies, cont

• We present the results of RMSE and related quantities for the 

expectation at the combination of x1 = 1, 1.5, 2 and x2 = 0.5 i

on 200 simulated data sets.
First simulation results on predicted means

x

1

x

2 Average value Method Bias SD 

0 0.0 6.5023 Proposed 0.1362 0.216

CTCV 0.0015 0.088

CTMV 0.0022 0.087MTCV 2.8103 0.018

CHIOU 0.8027 0.226

BASU Failed to co

0 1.0 8.7948 Proposed -0.0062 0.042

CTCV -0.0030 0.027

CTMV -0.0015 0.025MTCV 0.7702 0.022

CHIOU -0.3023 0.062

BASU Failed to co
First simulation results on predicted means, cont

x

1

x

2 Average value Method Bias SD 

1 0.0 8.9171 Proposed 0.0384 0.033

CTCV 0.0015 0.035

CTMV 0.0019 0.034MTCV 0.6480 0.005

CHIOU 0.1796 0.189

BASU Failed to co

1 1.0 9.8181 Proposed -0.0057 0.008

CTCV -0.0001 0.002

CTMV -0.0001 0.002MTCV -0.0724 0.004

CHIOU -0.0637 0.022

BASU Failed to co
Conclusion

• Among the two existing GLM-based estimators, we found that 

procedure failed to converge for all of 200 simulated samples, Baus & Rathouz’s estimator is not stable.

• Chiou and Muller’s estimator has much larger bias and SD tha

proposed estimator, and is badly biased.

• By comparing results among the parametric CTCV, CTMV, and 

we conclude that misspecification of the transformation functiobias and large RMSE while misspecification of variance functioon bias and RMSE.

• Comparing our new estimator with the gold standard estimator

estimator, derived under correctly specified transformation and the empirical efficiency of our new estimator is around 60% on 
Second simulation study

• In our second simulation study, we considers a true transforma

model with a logarithm transformation function, a constant variasymmetrical distribution error.

• For n = 2000 subjects, we generate covariates X1 and X2 fro

distribution with p = 0.5 and the uniform distribution on [0, 1], rrandom error ε from the Gamma distribution with scale of 1 an

• We let our outcome follow the following transformation model:

H(Y ) = β0 + X1β1 + X2β2 + ε/2,

where H(y) = 2 log(y), β0 = −8, β1 = 4, and β2 = 4.
Second simulation results on predicted means

x

1

x

2 Average cost Method Bias SD 

1 0.0 0.4287 Proposed 0.0276 0.019

CTCV -0.0012 0.0125 

MTCV 0.2068 0.0292 

BASU* -0.0035 0.0176 

CHIOU** 0.1308 0.0639 

1 0.5 1.1653 Proposed -0.0113 0.0496 

CTCV -0.0038 0.0238 

MTCV 0.2009 0.0418 

BASU* -0.0078 0.0265 

CHIOU** 0.1470 0.0466 

* The procedure failed to converge in 37 out of the 200 samples, the based on the 163 remaining samples. ** The procedure failed to conthe 200 samples, the summaries are based on the 177 remaining sa
An example

• The sample used here was from a study on the effectiveness o

Mood-Promoting Access to Collaborative Treatment (IMPACT) management program for late-life depression (Unutzer, et al., 2

• In this talk, we focus on the total outpatient cost in the first yea

standard deviation of Y are $6258.442 and $5065.507, respectcoefficients of skewness and kurtosis of Y are 3.36 and 26.94, 

• We fit the data using our model with the outcome variable bein

the first year, and the two independent variable, X1 and X2.

• Here X1 was the binary treatment indicator, and X2 was the m

depression items from the Symptom Checklist.
An example, cont

• We set the variance function to be a polynomial function σ2(x; 

p = 1, 2, · · · , where p was chosen to minimize the following GF

GF(p) = min

γ X

n

i=1 “H ̃ (Yi) − X′iβ ̃”2 − σ2(X′iγ

where H

e

and β ̃ were the initial values of H and β, respectively

• The results showed that GF(p) did not substantially change w• Hence to assure σ2(x) ≥ 0, here we selected σ2(x; p) = x2.
The estimates of the expected costs of a patient

Proposed CHIOU B

Rand SCL Expection(se) Expection(se) Expe

1 0.04 5008.6(444.4) 5239.9(74.7) 53340 0.04 4424.0(500.6) 4639.4(10.5) 49911 1.50 6916.1(392.9) 6779.1(26.7) 67170 1.50 6172.1(216.4) 6177.7(40.7) 6161 3.20 9177.3(1156.9) 8574.4(31.9) 98020 3.20 8349.5(816.3) 7971.8(97.5) 8622
The estimated transformation and its 95% confidence limits for IMPACT data.

0 20000 40000 60000

Y

-10 -5 0 5 10 15

H(Y)