Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: predict returning incorrect values for mu following glm


From   "Paul Dickman" <paul.dickman@ki.se>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: predict returning incorrect values for mu following glm
Date   Tue, 8 Nov 2005 11:48:01 +0100

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 (paul.dickman@mep.ki.se)
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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index