Roger Harbord <[email protected]> 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
[email protected]
------------------------------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/
