Leny Mathew wrote:
I'm trying to model data on change in
certain hormones over a period of time. The data on some of these have
a number of values that are set (currently) at a non-zero minimum
detectable value (MDL). Currently, I'm using a GLM with a gamma family
along with robust estimation of SE to model this data. The model works
reasonably fine on the data but I would like to improve this by
properly accounting for the information on the non-detectables.
I'm interested in the change over time of these hormones and so change
from the MDL to a higher value is useful information. Apart from the
GLM, I've looked at the tobit models and zero-inflated models. From
what I understand, the tobit model assumes that the data that is not
censored follows a normal distribution and this not the case in my
data. Zero-inflated models seem to be used for count rather than
continuous data.
Does anybody in the list have experience on modeling continuous data
with zero-inflation or data with a MDL issue?
--------------------------------------------------------------------------------
Take a look at the user-written module -intcens- by Jamie Griffin. It is
available from SSC (-ssc describe intcens-).
-intcens- will handle left-censored assay data (below limit of detection or
limit of quantitation) with a variety of distributions, including gamma.
It doesn't specifically handle panel (longitudinal) data, but it has
a -cluster()- option, which might be adequate for your needs, especially if
you're using -glm- for your analyses, anyway.
Joseph Coveney
clear *
set more off
set seed0 `=date("2008-05-11", "YMD")'
rndgam 1000 2 6
summarize xg, meanonly
display in smcl as text r(mean)
glm xg , family(gamma) nolog
display in smcl as text 1 / _b[_cons]
display in smcl as text 1 / e(phi) // Approx. estimate
display in smcl as text e(phi) / _b[_cons]
gammafit xg, nolog
clonevar a = xg
clonevar b = xg
intcens a b, distribution(gam) nolog
matrix define B = e(b)
display in smcl as text exp(B[1,2])
display in smcl as text exp(B[1,1] - B[1,2])
quietly centile xg, centile(5)
quietly replace a = .b if a < r(c_1)
label define LOQ .b "< LOQ"
label values a LOQ
intcens a b, distribution(gam) nolog
exit
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
