*! Version 1.0.4 2-24-93 J. Hilbe STB-12 sg16.1 * Generalized Linear Models : Power links program define glmp version 3.0 local options "LEvel(integer $S_level) EForm" local varlist "req ex" local options "`options' LTolerance(real -1) ITerate(integer 0) Fam(string) Pow(real 1.0) SCale(string) Stats Resid Offset(string) EXposure(string)" local in "opt" local if "opt" local weight "fweight aweight" parse "`*'" parse "`varlist'",parse(" ") qui { tempvar mu eta w z dev odev wt touse pz chi2 local y "`1'" mac shift * OFFSET/EXPOSURE OPTION tempvar LNEX if "`exposur'"~="" { if "`offset'"~= "" { error 198 } _crcunab `exposur' local expostr "$S_1" gen double `LNEX' = ln(`expostr') local offset `LNEX' local offmiss "| `offset'==." } else if "`offset'"~="" { _crcunab `offset' local offset "$S_1" local offmiss "| `offset'==." local expostr "exp($S_1)" } gen byte `touse'=1 `if' `in' replace `touse'=0 if `y'==. | `touse'==. `offmiss' if "`exp'"~="" { gen float `wt' `exp' if `touse' replace `wt'=. if `wt'<=0 replace `touse'=0 if `wt'==. } else gen float `wt' = `touse' if `touse' if "`offset'"~="" { replace `touse'=0 if `offset'==. replace `wt'=. if `offset'==. local poffset "+`offset'" local moffset "-`offset'" local offopt "offset(`offset')" } reg `y' `*' if `touse' predict `pz' replace `touse'=0 if `pz'==. drop `pz' local nobs = _result(1) if `ltolera'<0 { local ltolera 0.0001 } if `iterate'<=0 { local iterate 50 } * INITIALIZATION local i = 1 gen `w' = 0 if `touse' gen `z' = 0 if `touse' gen `dev' = 0 if `touse' gen `odev' = 1 if `touse' local ddev = 1 egen `mu' = mean(`y') if `touse' replace `mu' = (`y'+`mu')/2 if `pow'== 0 { gen `eta' = log(`mu') } else { gen `eta' = `mu'^`pow'} * LOCAL SCORING ALGORITHM (IRLS) while ((abs(`ddev') >`ltolera' & `i'<`iterate')| (`i'== 0)) { if `pow'== 0 { replace `w' = `mu' if `touse' } else { replace `w' = 1/(`pow'*`mu'^(`pow'-1)) if `touse'} replace `z' = `eta'+(`y'-`mu')/`w' `moffset' if `touse' if ("`fam'"=="gau" | "`fam'"=="") { replace `w'=`w'^2 } if "`fam'"=="poi" { replace `w' = `w'^2/`mu' } if "`fam'"=="gam" { replace `w' = `w'^2/`mu'^2 } if "`fam'"=="invg" { replace `w' = `w'^2/`mu'^3 } regress `z' `*' [iw=`w'*`wt'], mse1 dof(100000) dep(`y') cap drop `eta' predict `eta' if `touse' replace `eta' = `eta' `poffset' if `pow'== 0 { replace `mu' = exp(`eta') } else { replace `mu' = `eta'^(1/`pow') } replace `odev' = `dev' if ("`fam'" == "gau" | "`fam'" == "") { replace `dev' = (`y'-`mu')^2 } if "`fam'" == "poi" { replace `dev' = (`y'*log(`y'/`mu') - (`y'-`mu')) } if "`fam'" == "gam" { replace `dev' = -log(`y'/`mu') + (`y'-`mu')/`mu' } if "`fam'" == "invg" { replace `dev' = (`y'-`mu')^2 / (`mu'^2*`y') } replace `dev' = sum(`dev'*`wt') if ("`fam'"=="gau" | "`fam'"=="invg") { replace `dev' = `dev'[_N] } if ("`fam'"=="poi" | "`fam'"=="gam") { replace `dev' = 2*`dev'[_N] } local ddev = `dev'-`odev' noi di in gr "Iter `i' : Dev = " in ye %9.4f `dev' local i = `i' + 1 } * CHI SQUARE CALCULATION tempvar chi2 if "`fam'"=="gau" { gen `chi2' = (`y'-`mu')^2 if `touse' } if "`fam'"=="poi" { gen `chi2' = (`y'-`mu')^2/`mu' if `touse' } if "`fam'"=="gam" { gen `chi2' = (`y'-`mu')^2/`mu'^2 if `touse'} if "`fam'"=="invg" { gen `chi2' = (`y'-`mu')^2/`mu'^3 if `touse'} replace `chi2' = sum(`wt'*`chi2') replace `chi2' = `chi2'[_N] local df = `nobs'-_result(3)-1 local dispc = `chi2'/`df' local dispd = `dev'/`df' } di in gr _col(57) "No obs. = " in ye %9.0g `nobs' di in gr "Chi2 = " in ye %9.0g `chi2' in gr _col(57) "Deviance = " in ye %9.0g `dev' di in gr "Prob>chi2 = " in ye %9.0g chiprob(`df',`chi2') in gr _col(57) "Prob>chi2 = " in ye %9.0g chiprob(`df',`dev') di in gr "Dispersion = " in ye %9.0g `dispc' in gr _col(57) "Dispersion = " in ye %9.0g `dispd' _n if ("`fam'"=="gau" | "`fam'"=="") {di in gr "Gaussian: Power = `pow' "} if "`fam'"=="poi" {di in gr "Poisson: Power = `pow' "} if "`fam'"=="gam" {di in gr "Gamma: Power = `pow' "} if "`fam'"=="invg" {di in gr "Inverse Gaussian: Power = `pow' "} if `level'<10 | `level'>99 { local level 95 } if "`eform'"~="" { if ("`link'"=="l" & "$S_E_nc"=="") { local eform "eform(Odds Ratio)" } else if (("`link'"=="p" | "`link'"=="c") & "$S_E_nc"=="") { local eform "eform(ExpB)" } else if ("`fam'"=="poi" & "$S_E_nc"=="") { local eform "eform(IRR)" } else { local eform } } * SCALING AND DISPLAY OF RESULTS if "`scale'"~="n" { di in gr _dup(78) "-" di in gr _col(1) "`y'" _col(10) "|" _col(17) "Coef." _col(25) "Std.Err." _col(41) "t" _col(47) "P>|t|" _col(59) "[`level'% Conf. Interval]" di in gr _dup(78) "-" while "`1'"~="" { local j=1 if ("`scale'"=="d" | "`scale'"=="") { local se`j'=sqrt(`dispd')*_se[`1'] } if "`scale'"=="c" { local se`j'=sqrt(`dispc')*_se[`1'] } local t`j' = _b[`1']/`se`j'' local tpb`j' = tprob(`df',`t`j'') local lci`j' = _b[`1']-(invnorm(`level'/100+(1-`level'/100)/2))*(`se`j'') local uci`j' = _b[`1']+(invnorm(`level'/100+(1-`level'/100)/2))*(`se`j'') di in gr _col(1) "`1'" _col(10) "|" in ye _col(13) %9.0g _b[`1'] _col(24) %9.0g `se`j'' _col(35) %9.3f `t`j'' _col(47) %4.3f `tpb`j'' _col(58) %9.0g `lci`j'' _col(70) %9.0g `uci`j'' local j=`j'+1 mac shift } if ("`scale'"=="d" | "`scale'"=="") { local secon=sqrt(`dispd')*_se[_cons] } if "`scale'"=="c" { local secon=sqrt(`dispc')*_se[_cons] } local tcon = _b[_cons]/`secon' local tpbcon = tprob(`df',`tcon') local lcicon = _b[_cons]-(invnorm(`level'/100+(1-`level'/100)/2))*(`secon') local ucicon = _b[_cons]+(invnorm(`level'/100+(1-`level'/100)/2))*(`secon') di in gr _col(1) "_cons" _col(10) "|" in ye _col(13) %9.0g _b[_cons] _col(24) %9.0g `secon' _col(35) %9.3f `tcon' _col(47) %4.3f `tpbcon' _col(58) %9.0g `lcicon' _col(70) %9.0g `ucicon' di in gr _dup(78) "-" di in gr "Standard errors adjusted by sqrt(dispersion)" } else { noi regress, noheader level(`level') `eform' } * CALCULATION OF LINEAR PREDICTOR (ETA) AND PREDICTED (MU) if "`stats'"~="" { qui { gen _eta = `eta' if `touse' gen _mu = `mu' if `touse' } di in gr "Variables created: " in ye " _eta, _mu" lab var _eta "Linear predictor" lab var _mu "Predicted" } if "`resid'"~="" { qui { if ("`fam'"=="bin" & "`group'"=="") { gen ___c = 1 if `touse'} else if ("`fam'"=="poi") { gen ___c = 3 if `touse' } else if ("`fam'"=="gam") { gen ___c = 4 if `touse'} else if ("`fam'"=="invg") { gen ___c = 5 if `touse' } else if ("`fam'"=="gau") { gen ___c = 6 if `touse' } else if ("`fam'"=="bin" & "`group'"~= "") { gen ___m = `m' if `touse' gen ___c = 2 if `touse' } gen ___y = `y' if `touse' gen ___mu = `mu' if `touse' gen ___w = `wt' if `touse' } _glmresd } end