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

st: -swaic- for Stata V7+ -- Attn: Lucas Habib


From   "Steichen, Thomas J." <SteichT@rjrt.com>
To   <statalist@hsphsun2.harvard.edu>
Subject   st: -swaic- for Stata V7+ -- Attn: Lucas Habib
Date   Mon, 14 Mar 2005 08:33:06 -0500

My apologies to the list.
  
I received a private request for an older program that I had modified to run 
on Version 7 and later.  Three attempts to reply directly have all ended in 
failure. As a last and alternative attempt, I'm posting the program the old-
fashioned way, embedded in Statalist, in hopes Lucas will receive it this way.

To Lucas:

The program and help file are embedded below. Note that the name has been
changed to -swaic7-.  You should cut the program from this message and save
it as swaic7.ado in your personal ado directory.  Do similarly with swaic7.hlp.
You may need to rejoin some long lines if they are split during mailing.
 
Tom

> -----Original Message-----
> From: Lucas Habib [mailto:lucas.habib@ualberta.ca]
> Sent: Monday, March 07, 2005 10:26 PM
> To: Steichen, Thomas J.
> Subject: swaic and xi3
> 
> Hi Tom,
> I came across your posting to the STATA list from September 2004 where
> you mentioned that you had a version of swaic that would work 
> with STATA 8.  Would you be able to send it to me?

------- swaic7.ado ----cut here---------------------------------------------
*! version 1.3 TJS    24sep2004 xi: works under both pre and post version 7
*  version 1.2 TJS    23sep2004 added model list, min aic output
*  version 1.1 Z.WANG 22NOV1999 (STB-54: sg134)
*  version 1.0 Z.WANG 07Nov1999 
program define swaic7, rclass     /* changed to rclass */
  version 6.0
  syntax [, FP(str) FChi(str) Back Model ]
  local yvar "`e(depvar)'"

* v7 fix
  global caller = _caller()      /* added */
* di in red "caller: " $caller
    
  if "`e(cmd)'" == "" {
    di in r "last estimates not found"
    exit 301
  } 
  if index(`"logitlogisticpoissonprobit"', `"`e(cmd)'"') > 0 {
    local cmd "`e(cmd)'"
    local yvar1 yvar(`yvar')
  }    
  if "`e(cmd)'" == "cox" {
    if "`e(cmd2)'" != "" {
      local cmd `e(cmd2)'
      local yvar1
      local yvar
    } 
    else {
      di in gr "Please use stcox instead of cox in estimation command."
      exit 199
    }
  }
  if `"`e(cmd2)'"' == `"streg"' {
    local dist `"dist(`e(cmd)')"' 
    local cmd `e(cmd2)'
    local yvar1
    local yvar
  } 
  if "`cmd'" == "" {
    di in w "`e(cmd)' " in gr "not supported by " in w "swaic"
    di in gr _con `"Current"' in w `" swaic"' in gr `" only supports "' 
    di in w  _con "logit logistic poisson stcox streg" 
    di in gr `" and"' in w " probit" in gr "."
    exit 199
  }
  if "`e(offset)'" != "" {
    tempvar off 
    gen `off'=`e(offset)'
    lab var `off' "`e(offset)'"
    local offset offset(`off')
  }
  if "`e(vcetype)'" ~= "" { 
    local robust r
  }
  if "`e(clustvar)'" ~= "" {
    local cluster cluster(`e(clustvar)')
  }
  local ops `offset' `robust' `cluster' `dist'

  tempvar smpl
  qui gen `smpl' = e(sample)
  mat bs = e(b) 
  local xnames : colnames(bs) 
  _xnam `xnames'
  
  local xnames  `r(alln)'
  local vnum  = `r(vnum)'
  local factn = `r(factn)'  

  di in gr "Stepwise Model Selection by AIC"  
  di in gr "`cmd' regression. `dist'"
  di in gr "number of obs = " in ye e(N)
  di in gr _dup(71) "-" /* formatting changes added for longer var names */
/*           outcome    Df     Chi2   P>Chi2    -2*ll   Df Res.  AIC     */  
  di in gr _col(9) %12s "`yvar'" _col(25) "Df" _col(32) "Chi2" /*    
    */ _col(39) "P>Chi2" _col(49) "-2*ll" _col(57) "Df Res." _col(66) "AIC"
  di in gr _dup(71) "-"
  token `xnames' 
* Display Table
  if "`fp'" == "" {
    local fp "%9.4f" 
  }
  if "`fchi'" == "" {
    local fchi "%9.2f" 
  }

  local decreas = 0
  local endhere = 0
* backward stepping  
  if "`back'" ~= "" {
    qui `cmd' `yvar' `xnames' if `smpl', `ops'
    local aic0 = -2 * e(ll) + 2 * (e(df_m) + `factn')
    local dev0 = -2 * e(ll)
    local rdf0 = e(N) - e(df_m) - `factn'
    local mdf0 = e(df_m)
    di in gr "Full Model" in ye _col(46) `fchi' `dev0' _col(56) /*
      */ %6.0f `rdf0' _col(60) `fchi' `aic0'
    local i = 1
    local xs `xnames'
    while `i' <= `vnum' {
      local i_1 = `i' - 1
      _stepb, `yvar1' xvar(`xs') cmd(`cmd') smpl(`smpl') /*
        */ factn(`factn')  `ops' 
      local pv`i'  `r(pv)'
      local mdf`i' `r(mdf)'
      local aic`i' `r(aic)'
      local dev`i' `r(dev)'
      local chi`i' = `dev`i'' - `dev`i_1''
      local df`i'  = `mdf`i_1'' - `mdf`i''
      local p`i'   = chiprob(`df`i'', `chi`i'')
      local rdf`i' = `rdf`i_1'' + `df`i''
      if `decreas' == 1 {
        if `aic`i'' > `aic`i_1'' { 
          local endhere = 1
        }  
      } 
      if `aic`i'' < `aic`i_1'' {
        if `endhere' == 0 {
          local maic = `aic`i''               /* added */
          local mlist `r(mlist)'
        }
        local decreas = 1
      }
     di in gr " Step `i':" _col(9) %12s "-`pv`i''" in ye %6.0f /*
        */ _col(14) `df`i'' _col(21) `fchi' `chi`i'' _col(33)  /*
        */ `fp' `p`i'' _col(46) `fchi' `dev`i'' _col(56) %6.0f /*
        */ `rdf`i'' _col(60) `fchi' `aic`i''

      _newx, xvar(`xs') pick(`r(pv)')
      local xs `r(newlst)'
      local i = `i' + 1
    }
  }
* forward stepping  
  else {
    qui `cmd' `yvar' `xnames' if `smpl', `ops'
    local aic0 = -2 * e(ll_0) + 2 * `factn'
    local dev0 = -2 * e(ll_0)
    local rdf0 = e(N) - `factn'
    local mdf0 = 0
    local df0  = 0
    di in gr "Null Model" in ye _col(46) `fchi' `dev0' /*
      */ _col(56) %6.0f `rdf0' _col(60) `fchi' `aic0'
    local xs `xnames'
    local i = 1
    while `i' <= `vnum' { 
      local i_1 = `i' - 1
      _stepf, `yvar1' xvar(`"`xs'"') cmd(`cmd') smpl(`smpl') /*
        */ `plist' factn(`factn') `ops' 
      local pv`i'    `r(pv)'
      local mdf`i' = `r(mdf)'
      local aic`i' = `r(aic)'
      local dev`i' = `r(dev)'
      local chi`i' = `dev`i_1'' - `dev`i''
      local df`i'  = `mdf`i'' - `mdf`i_1''
      local p`i'   = chiprob(`df`i'', `chi`i'')
      local rdf`i' = `rdf0' - `mdf`i''
      local picked   `picked' `pv`i''
      if `decreas' == 1 {
        if `aic`i'' > `aic`i_1'' { 
          local endhere = 1
        }  
      } 
      if `aic`i'' < `aic`i_1'' {
        local maic = `aic`i''                 /* added */
        if `endhere' == 0 {
          local mlist `picked'
        }
        local decreas = 1
      }
      di in gr " Step `i':" _col(9) %12s "`pv`i''" in ye %6.0f /*
        */ _col(14) `df`i'' _col(21) `fchi' `chi`i'' _col(33)  /*
        */ `fp' `p`i'' _col(46) `fchi' `dev`i'' _col(56) %6.0f /*
        */ `rdf`i'' _col(60) `fchi' `aic`i''

      local plist plist(`picked')
      _newx, xvar(`xs') pick(`pv`i'')
      local xs `r(newlst)'
      local i = `i' + 1
    } 
  }

  di in gr _dup(71) "-"

  if "`model'" != "" { 
    `cmd' `yvar' `mlist' if `smpl', `options' nolog
  }
  di _n as txt "minimun AIC =" as res %9.3f `maic' /*      added
  */    as txt ";  model:" as res " `mlist'"            /* added */
  qui `cmd' `yvar' `xnames' if `smpl', `options'

  return local model   = "`mlist'"            /* added */
  return local min_aic = `maic'               /* added */

end

program define _stepb, rclass 
  syntax [, Yvar(str) (str Xvar(str) CMD(str) smpl(str) factn(str) *]
  token `xvar'
  local aicmin = 10000
  local mdf    = .
  local dev    = .
  local i = 1
  while "``i''" ~= "" {
    local i_1 = `i' - 1  
    _newx, xvar(`xvar') pick(``i'') 
    local newx `r(newlst)'
    qui `cmd' `yvar' `newx' if `smpl', `options'
    local dev`i' = -2 * e(ll)
    local df`i' = e(df_m)   
    local aic`i' = `dev`i'' + 2 * (`df`i'' + `factn')
    if `aic`i'' < `aicmin' {
      local pv       "``i''"
      local mdf    = `df`i''
      local dev    = `dev`i''
      local mlist    `newx'
      local aicmin = `aic`i''
    }
    local i = `i' + 1
  }
  return local pv    `pv'
  return local aic = `aicmin'
  return local mdf = `mdf'
  return local dev = `dev'
  return local mlist `mlist'
end

program define _stepf, rclass 
  syntax [, Yvar(str) plist(str) Xvar(str) factn(str) CMD(str) smpl(str) *]
  token `xvar'
  local aicmin = 10000
  local mdf    = .
  local dev    = .
  local i = 1
  while "``i''" ~= "" {
    local i_1 = `i' - 1  
    qui `cmd' `yvar' `plist' ``i'' if `smpl', `options'
    local dev`i' = -2 * e(ll)
    local aic`i' = `dev`i'' + 2 * (e(df_m) + `factn')
    if `aic`i'' < `aicmin' {
      local pv       "``i''"
      local mdf    = e(df_m)
      local dev    = `dev`i''
      local aicmin = `aic`i''
    }
    local i = `i' + 1
  }
  return local pv    `pv'
  return local aic = `aicmin'
  return local mdf = `mdf'
  return local dev = `dev'
end

program define _ab, rclass
  token `0', parse(" _")
  if $caller < 7 { unab cat:  `1'* }          /* added */
    else         { unab cat: _`1'* }   
  return local cname "`1'*"
  token `cat'
  local i = 1
  while "``i''" != "" {
    local i = `i' + 1
  }
  return scalar catdf = `i' - 1
end

program define _xnam, rclass
  local factn = 0
  while "`1'" != "" {
    if $caller < 7 {                          /* added  */
      local if1 = substr("`1'", 1, 1)
      local if2 = "_"
      local i1  = 1
      local if3 = "I"
      local u   = ""
    }
    else {
      local if1 = "`1'"
      local if2 = "_cons"
      local i1  = 2
      local if3 = "_I"
      local u   = "_"
    }
    if "`if1'" != "`if2'" {                   /* modified */
      local xlist `xlist' `1'
    }
    else {
      local factn = `factn' + 1
    }
    mac shift
  }
  token `xlist'
  local i = 1 
  while "`1'" != "" {
    if substr("`1'", 1, `i1') == "`if3'" {    /* modified */
      if $caller < 7 { local x = "`1'" }      /* added  */
        else         { local x = substr("`1'", 2, .) }
      _ab `x'                                 /* modified */
      local catdf = r(catdf)
      local xnam`i' `u'`r(cname)'             /* modified */
      mac shift `catdf'
    }
    else {
      local xnam`i' `1'
      mac shift
    }
    local alln `alln' `xnam`i''
    local i = `i' + 1
  }
  return local alln `alln'
  return local vnum = `i' - 1
  return local factn `factn'
end  

program define _newx, rclass
  syntax [, Xvar(str) Pick(str)]  
  token `xvar'
  while "`1'" != "" {
    if  "`1'" != "`pick'" {
      local newlst "`newlst' `1'"
    }
    mac shift
  }
  return local newlst `newlst'
end

---- swaic7.hlp -------cut here--------------------------------------------------
..-
help for ^swaic7^                 modified version of swaic from (STB-54: sg134)
..-

Stepwise model selection using AIC
----------------------------------

	^swaic7^ [^,^ ^fp(^format^)^ ^fc^hi^(^format^)^ ^m^odel ^b^ack ]

Description
-----------

 ^swaic7^ performs automatic model selection using Akaike information 
 criterion (AIC), where 
 
     AIC = -2 * log-likelihood + 2 * (number of parameters)

 ^swaic7^ is identical to ^swaic^ except that it returns the selected model
 and its AIC value in r(model) and r(min_aic), and it works correctly
 with ^xi:^ regardless of Stata version number.

 
Options
-------
^fp(^format^)^ specifies the output format for p values, default is %9.4f. 

^fc(^format^)^ specifies the output format for chi2 values, default is %9.2f. 

^model^ reports the model which reachs the minimum AIC with either forward
    or backward method.

^back^ specifies a backward method, the default is a forward method.


Examples
--------

   . xi: logit outcome age sex i.expose hibp bmi
   . swaic7, m

   . swaic7, b m

   . stset time outcome
   . xi: stcox age sex bmi hibp i.expose, nolog
   . swaic7, b m

   . xi:streg age sex bmi hibp i.expose, dist(weibull)
   . swaic7, b m

   . xi: poisson outcome age sex bmi hibp i.expose, e(time) nolog
   . swaic7, m


Author of swaic
---------------

    Zhiqiang Wang
    Menzies School of Health Research
    Darwin, Australia
    wang@@menzies.edu.au

    
Author of swaic7    
----------------

    T. J. Steichen
    steicht@@rjrt.com    

    
Also see
--------

 Manual:  ^[U] 23 Estimation and post-estimation commands^
	  ^[R] lrtest^
     
On-line:  help for @linktest@, @test@, @testnl@, @lrtest2@ (if installed) @streg@

-----------cut here-----------------------------------------------------------

-----------------------------------------
CONFIDENTIALITY NOTE:  This e-mail message, including any  attachment(s),
contains information that may be confidential,  protected by the attorney-
client or other legal privileges, and/or  proprietary non-public
information.  If you are not an intended  recipient of this message or an
authorized assistant to an intended  recipient, please notify the sender by
replying to this message and  then delete it from your system.  Use,
dissemination, distribution,  or reproduction of this message and/or any of
its attachments (if  any) by unintended recipients is not authorized and
may be unlawful.


*
*   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