predict is giving me incorrect values following glm. Following is some code
that illustrates the problem. rs.ado (which defines the model) is at the end
of this message.
This is a Poisson regression model for excess mortality.
ln(mu-d*)=ln(y)+xbeta
mu=E(d) is assumed Poisson. If we didn't have the "-d*" component in the
link function (d* is expected deaths) then this would be a standard Poisson
regression model.
I have defined "`mu' = exp(`eta')+$SGLM_p" but predict seems to be ignoring
the $SGLM_p component. Where have I gone wrong?
This mail was motivated by problems Enzo Coviello brought to my attention
with a similar model I had coded in Stata.
clear
input end d d_star y betahat
1 49 21 468 0
2 52 20 421 .2393669
3 47 19 373 .2268899
4 33 17 332 -.2162825
5 39 16 296 .2613985
end
xi: glm d i.end, fam(pois) link(rs d_star) lnoffset(y)
predict xbeta, xb
predict mu, mu
scal constant=-2.816264
gen xbeta2=constant+betahat+ln(y)
gen mu2=exp(xbeta2)+d_star
list end d d_star y xbeta xbeta2 mu mu2
. list end d d_star y xbeta xbeta2 mu mu2
+---------------------------------------------------------------+
| end d d_star y xbeta xbeta2 mu mu2 |
|---------------------------------------------------------------|
1. | 1 49 21 468 3.332205 3.332204 28 49 |
2. | 2 52 20 421 3.465736 3.465736 32 51.99999 |
3. | 3 47 19 373 3.332205 3.332204 28 47 |
4. | 4 33 17 332 2.772589 2.772588 16 33 |
5. | 5 39 16 296 3.135494 3.135494 23 39 |
+---------------------------------------------------------------+
This is a saturated model so the predicted values should equal the observed
values. The values of xbeta given by predict are correct but the values of
mu are not. It is apparent that predict is calculating
mu=exp(xbeta)
when I was hoping that it would calculate
mu=exp(xbeta)+d_star
Here are details of my setup.
. version
version 9.1
. which glm
C:\Program Files\Stata9\ado\updates\g\glm.ado
*! version 5.6.14 04aug2005
. which predict
C:\Program Files\Stata9\ado\base\p\predict.ado
*! version 2.0.4 24sep2004
I have defined "`mu' = exp(`eta')+$SGLM_p" but predict seems to be ignoring
the $SGLM_p component. Following is my code for defining the model. Where
have I gone wrong?
program define rs
version 7
args todo eta mu return
if `todo' == -1 {
global SGLM_lt "Relative survival"
global SGLM_lf "log(u-d*)"
exit
}
if `todo' == 0 { /* eta = g(mu) */
gen double `eta' = ln(`mu'-$SGLM_p)
exit
}
if `todo' == 1 { /* mu = g^-1(eta) */
gen double `mu' = exp(`eta')+$SGLM_p
exit
}
if `todo' == 2 { /* (d mu)/(d eta) */
gen double `return' = exp(`eta')
exit
}
if `todo' == 3 { /* (d^2 mu)(d eta^2) */
gen double `return' = exp(`eta')
exit
}
di as error "Unknown call to glm link function"
exit 198
end
-----
Paul Dickman ([email protected])
Department of Medical Epidemiology and Biostatistics,
Karolinska Institutet, Box 281, 171 77 Stockholm, Sweden
Ph: +46 8 5248 6186 Fax: +46 8 314 975
*
* 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/
