Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Binomial regression


From   rgutierrez@stata.com (Roberto G. Gutierrez, StataCorp LP)
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Binomial regression
Date   Sat, 04 Aug 2007 10:01:43 -0500

Roger Harbord <rogerharbord@bigfoot.com> writes:

> Thanks Bobby, very interesting.  Personally I'm particularly interested in
> your mention that this is the result of similar work for the log link -
> presumably a log/logit link? Is this available somewhere, or will it be?
 
As of five minutes ago, it was only available in my home directory.  It was
just something I kicked around a while back, someday to be written up.  

The -loglogit- link code is appended below.  Since the range of the inverse
log link is (0, infty), there is only one endpoint to worry about.  By
default, the loglogit link is equivalent to the log link in the range (0,0.99)
and logit past 0.99.  As such, the interpretation is that the resulting
coefficients are log(relative risk)'s for predicted probabilities up to 0.99,
log(odds ratio)'s in the extreme upper tail.

As was the case with -identitylogit-, you can override the default behavior
with an optional link argument p, which makes the valid log range (0,p) with
logit past p.

-Bobby
rgutierrez@stata.com

------------------------------begin loglogit.ado-----------------------------
*! version 1.1.0   04aug2007
program define loglogit		
	version 8.2 
	args todo eta mu return

        if `todo' == -1 {                       /* Title */
		if "$SGLM_p" == "" {
			global SGLM_p 0.99
		}
		else {
			confirm number $SGLM_p
			if $SGLM_p <= 0 | $SGLM_p >= 1 {
				di as error ///
			            "link parameter p must be between 0 and 1"
				exit 459
			}
		}
		if "$SGLM_m" == "1" {
			global SGLM_lf "ln(u) -> logit(u)"
		}
		else    {
			global SGLM_lf "ln(u/$SGLM_m) -> logit(u/$SGLM_m)"
		}
                global SGLM_lt "Log to $SGLM_p, then logit"
                exit
        }

	tempname a b c
	scalar `a' = 1/(1-$SGLM_p)
	scalar `b' = logit($SGLM_p) - `a'*ln($SGLM_p)
	scalar `c' = ln($SGLM_p)

	if `todo' == 0 {			// eta = g(mu) 
		gen double `eta' = ln(`mu'/$SGLM_m)
		replace `eta' = (logit(`mu'/$SGLM_m)-`b')/`a'	///
		 	  	if `eta' >= `c'
		exit 
	}
	if `todo' == 1 {			// mu = g^-1(eta) 
		gen double `mu' = $SGLM_m*exp(`eta')
		replace `mu' = $SGLM_m*invlogit(`a'*`eta'+`b') 	///
		 	  	if `eta' >= `c'
		exit 
	}
	if `todo' == 2 {			// (d mu)/(d eta) 
		gen double `return' = `mu'
		replace `return' = (`mu'*(1-`mu'/$SGLM_m)) * `a' ///
		 	  	if `eta' >= `c'
		exit 
	}
	if `todo' == 3 {			// (d^2 mu)(d eta^2) 
		gen double `return' = `mu'
		replace `return' = (`mu'*(1-`mu'/$SGLM_m)* ///
			   (1-2*`mu'/$SGLM_m)) * (`a'^2)   ///
			   if `eta' >= `c'
		exit
	}
	noi di as err "Unknown call to glm link function"
	exit 198
end
--------------------------------end loglogit.ado-----------------------------
*
*   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