*! version 1.01.2 21august2008 *! Time-stamp: <2009-01-27 14:51:25 prathouz> program define xtlogitods , eclass version 10 , missing // Some issues: // - might not want to use double precision with interactions // Preliminaries syntax varlist(min=1 num ts) [if] [in] [iweight fweight pweight] , /// CASE(varname num) /// [ SAMPLERATio(real 1) SAMPLERATio1(varname num) /// CASEEQ1(varlist num ts) /// CASEEQ2(varlist num ts) CASEOFFset(varname num) /// Corr(string) OFFset(varname num) /// PRINTINT /// *] // Reconcile sample ratio as real or as varname and generate variable // `samprat' for sample ratio no matter what was passed to program // - Note: See stata _help undocumented , parse_ for an easier way if `"`sampleratio1'"'!="" { confirm variable `sampleratio1' local sampleratio `sampleratio1' } tempvar samprat qui gen double `samprat'=`sampleratio' // Marksample and check marksample touse markout `touse' `case' `samprat' `caseeq1' `caseeq2' /// `offset' `caseoffset' // Trick to unmark observations that will be eliminated // by -xtgee- procedure // - Note: Could add message here about dropped observations, // as in -xtgee- capture : xtlogit `varlist' if `touse' [`weight'`exp'] , /// pa c("`corr'") `options' iter(1) tempvar xtgeesample qui gen double `xtgeesample'=e(sample) qui replace `xtgeesample'=. if `xtgeesample'==0 markout `touse' `xtgeesample' // test with -tab `touse' , miss- // Parse to separate response from predictors // - Note: could use -gettoken- to produce two macros // along these lines tokenize `varlist' local resp `1' macro shift local predictors `*' // Check macros formed correctly qui di " `resp' `predictors' " qui di " `case' `caseeq1' " qui di " `caseeq2' " // Set printing macro local print qui if("`printint'" != "") { local print } // Generate predictors (including interections w response) // for case variable capture : drop _Iy* if (_rc==0) { di " " di as text "Note: Dropping all variables with prefix _Iy" } qui gen double _Iy=`resp' local caseeq2y if ("`caseeq2'" != "") { foreach X of varlist `caseeq2' { qui gen double _IyX`X'=`X'*_Iy local caseeq2y `caseeq2y' _IyX`X' } } // Generate offset term for case-model tempvar casebiascorr qui gen double `casebiascorr'=log(`samprat') if ("`caseoffset'" != "") { replace `casebiascorr' = `casebiascorr' + `caseoffset' } // Fit model for case status tempvar lambdaS lambdaP1 lambdaP0 F1 F0 //di //di as text "GEE model for case status:" //xtlogit `case' `caseeq1' _Iy `caseeq2y' if `touse' [`weight'`exp'] , /// // pa corr(ind) offset(`casebiascorr') `options' `print' di `print' di as text "GEE model for case status:" `print' logit `case' `caseeq1' _Iy `caseeq2y' /// if `touse' [`weight'`exp'] , /// offset(`casebiascorr') asis // Print eigen values for Avar matrix if (0) { mata : eigensystem(st_matrix("e(V)"), evecs=., evals=.) mata : "Evals and cond number for Avar for GEE case model" mata : evals mata : evals[1,1]/evals[1,cols(evals)] } tempname gamma // save coefficients for case status model matrix `gamma'=e(b) matrix colnames `gamma'= case: // add equation prefix to colnames qui predict double `lambdaS' if `touse' // need for constructing score // Get predicted values for response=0 or=1 and then reset predicted values qui replace _Iy=1 if ("`caseeq2'" != "") { foreach X of varlist `caseeq2' { qui replace _IyX`X'=`X'*_Iy } } qui predict double `lambdaP1' if `touse' , xb nooff if ("`caseoffset'" != "") { qui replace `lambdaP1'=`lambdaP1' + `caseoffset' } qui replace `lambdaP1'=1/(1+exp(-`lambdaP1')) qui replace _Iy=0 if ("`caseeq2'" != "") { foreach X of varlist `caseeq2' { qui replace _IyX`X'=`X'*_Iy } } qui predict double `lambdaP0' if `touse' , xb nooff if ("`caseoffset'" != "") { qui replace `lambdaP0'=`lambdaP0' + `caseoffset' } qui replace `lambdaP0'=1/(1+exp(-`lambdaP0')) qui replace _Iy=`resp' if ("`caseeq2'" != "") { foreach X of varlist `caseeq2' { qui replace _IyX`X'=`X'*_Iy } } // Compute F1 and F0, objects needed in computation of // (partal Bi/partial gamma) qui gen double `F1'=`lambdaP1'*(1-`lambdaP1')*(1-`samprat')/ /// (1-(1-`samprat')*`lambdaP1') qui gen double `F0'=`lambdaP0'*(1-`lambdaP0')*(1-`samprat')/ /// (1-(1-`samprat')*`lambdaP0') // Compute bias-correcting offset and add to user-specified offset tempvar biascorr qui gen double `biascorr'= /// log( (1 - `lambdaP1' + `samprat'*`lambdaP1') / /// (1 - `lambdaP0' + `samprat'*`lambdaP0') ) if ("`offset'" != "") { qui replace `biascorr' = `biascorr' + `offset' } // Fit main model tempvar muS `print' di `print' di as text "GEE model for response of interest:" `print' xtlogit `resp' `predictors' if `touse' [`weight'`exp'] , /// pa offset(`biascorr') c("`corr'") `options' // estat wcorrelation , compact // Print eigen values for Avar matrix if (0) { mata : eigensystem(st_matrix("e(V)"), evecs=., evals=.) mata : "Evals and cond number for Avar for main GEE model" mata : evals mata : evals[1,1]/evals[1,cols(evals)] } tempname beta matrix `beta'=e(b) // save coefficients for main model matrix colnames `beta'=`resp': // add equation prefix to colnames // Note: Will not work to leave eqs blank for beta when do // ereturn post below qui predict double `muS' if `touse' // need for constructing score // Compute robust standard errors using Mata tempname Avar mata: robust_se() // Combine coeffiencients into a long vector and set names for // beta and Avar matrix `beta'=(`gamma' , `beta') mata: set_names() // Save important results in e() from xtlogit run, repost, // and then display results local scalist N N_g g_max g_min g_avg local maclist ivar tvar corr depvar foreach X in `scalist' { // scalars local `X' = e(`X') } foreach X in `maclist' { // macros local `X' = e(`X') } tempname R matrix `R'=e(R) ereturn clear ereturn post `beta' `Avar' , depname("`resp'") esample("`touse'") foreach X in `scalist' { // scalars ereturn scalar `X'=``X'' } foreach X in `maclist' { // macros ereturn local `X'="``X''" } ereturn matrix R=`R' ereturn local sampleratio="`sampleratio'" __xtlogitods_display end // Mata program to compute robust standard errors version 10 mata: void robust_se() { // Get panel information st_view(id, . , st_global("_dta[iis]") , st_local("touse")) info=panelsetup(id, 1) /// check w -panelstats(info)- // Get response and predictor matrices (see paper for notation) st_view(Z, . , st_local("case"), st_local("touse")) st_view(W1, . , tokens(st_local("caseeq1")) , st_local("touse")) st_view(W2, . , tokens(st_local("caseeq2")) , st_local("touse")) st_view(Y, . , st_local("resp"), st_local("touse")) st_view(X, . , tokens(st_local("predictors")), st_local("touse")) st_view(lambdaS, . , st_local("lambdaS"), st_local("touse")) st_view(F1, . , st_local("F1"), st_local("touse")) st_view(F0, . , st_local("F0"), st_local("touse")) st_view(muS, . , st_local("muS"), st_local("touse")) // Get fitted correlation matrix C = st_matrix("e(R)") // create panel-level score matrices info=panelsetup(id, 1) Meat=J(cols(W1)+cols(W2)+2+cols(X)+1, /// cols(W1)+cols(W2)+2+cols(X)+1,0) Bread=Meat for (i=1; i<=rows(info); i++) { // Objects in this section correspond to notation in Section 3 // of the paper Zi=panelsubmatrix(Z,i,info) ni=rows(Zi) Yi=panelsubmatrix(Y,i,info) Wi=(panelsubmatrix(W1,i,info) , Yi, Yi:*panelsubmatrix(W2,i,info) , J(ni,1,1) ) /* last element is _cons */ // Following two objects, W1istar and W0istar, are two Wi-matrices // in (partial Bi/partial gamma) W1istar=(panelsubmatrix(W1,i,info) , J(ni,1,1), panelsubmatrix(W2,i,info) , J(ni,1,1) ) W0istar=(panelsubmatrix(W1,i,info) , J(ni,1,0), J(ni,cols(W2),0) , J(ni,1,1) ) // Following two objects, Fi1 and Fi0, also appear // in (partial Bi/partial gamma) F1i=diag(panelsubmatrix(F1,i,info)) F0i=diag(panelsubmatrix(F0,i,info)) Xi=(panelsubmatrix(X,i,info),J(ni,1,1)) lambdaSi=panelsubmatrix(lambdaS,i,info) muSi=panelsubmatrix(muS,i,info) Ci=C[|1,1 \ ni,ni|] // upper left ni x ni correlation matrix // Print objects for first few panels if ((i<=1) & (0)) { "i is:" i "Zi, Yi, Wi, Xi, lambdaSi; muSi" Zi; Yi; Wi; Xi; lambdaSi; muSi "Ci" Ci "Wi W1istar W0istar" Wi; W1istar; W0istar "F1i, F0i" F1i; F0i } // Compute scores // gamma score Ti=Wi'*(Zi-lambdaSi) // beta score Aisqrt=(diag(sqrt(muSi:*(J(ni,1,1)-muSi)))) Ui=Xi'*Aisqrt*cholinv(Ci)*cholinv(Aisqrt)*(Yi-muSi) // Note: -cholinv()- is inverse for positive semi-definite // symmetric matrices. See manual for other inverse // functions // complete score vector Scorei=(Ti \ Ui) // Information // gamma-gamma information Itti=Wi'*(diag(lambdaSi:*(J(ni,1,1)-lambdaSi)))*Wi // gamma-beta information Itui=J(cols(Wi),cols(Xi),0) // beta-gamma information // Iuti=Itui' /* Temporary 0 for testing */ partialBipartialgammaT=(-W1istar'*F1i + W0istar'*F0i)' Iuti=Xi'*Aisqrt*cholinv(Ci)*Aisqrt*partialBipartialgammaT // beta-beta information Iuui=Xi'*Aisqrt*cholinv(Ci)*Aisqrt*Xi // complete information matrix Ii = ( (Itti , Itui) \ (Iuti , Iuui) ) // Print objects for first few panels if ((i<=2) & (0)) { "i Meat Bread" i Meat Bread } // Add quantities to total Meat = Meat + (Scorei * Scorei' ) Bread = Bread + Ii } // Compute variance matrix Avar = luinv(Bread)*Meat*(luinv(Bread))' // Check eigenvalues if (0) { "Bread evals and cond number" eigensystem(Bread,Brevec,Brevals) Brevals Brevals[1,1]/Brevals[1,cols(Brevals)] "Meat evals and cond number" eigensystem(Meat,Meevec,Meevals) Meevals Meevals[1,1]/Meevals[1,cols(Meevals)] "Avar evals and cond number" eigensystem(Avar,Avevec,Avevals,rcond=1) Avevals Avevals[1,1]/Avevals[1,cols(Avevals)] rcond } // Compute standard errors as would have been done by GEE // to test rescaling by stata if (0) { sevec = sqrt(diagonal(Avar)) sevec sqrt(1897/(1897-13))*sevec[1::cols(Wi),.] sqrt(1897/(1897-6))*sevec[(cols(Wi)+1)::(cols(Wi)+cols(Xi)),.] } // Return variance estimate to stata st_matrix(st_local("Avar"),Avar) } void set_names() { // Get colstripe from beta colstripe = st_matrixcolstripe(st_local("beta")) // Assign new stripe values st_matrixcolstripe(st_local("Avar"),colstripe) st_matrixrowstripe(st_local("Avar"),colstripe) } end // Program to display results program define __xtlogitods_display version 10 // code swiped and adapted from geeprt in xtgee.ado (geerpt) // for posting results noi di "" local q "_col(68)" local l = 71 + 8 - length("`e(ivar)'") noi di in green "Logistic ODS population-averaged model" /// in gr _col(49) "Number of obs" `q' "=" in ye _col(70) %9.0g e(N) local ivar = abbrev("`e(ivar)'", 12) if "`e(tvar)'"!="" { local tvar = abbrev("`e(tvar)'", 12) local skip = max(43-length("`ivar' `tvar'"),1) di in gr "Group and time vars:" in ye _col(`skip') /// "`ivar' `tvar'" /// in gr _col(49) "Number of groups" `q' "=" /// in ye _col(70) %9.0g e(N_g) } else { local skip = max(43-length("`ivar'"),1) di in gr "Group variable:" in ye _col(`skip') /// "`ivar'" /// in gr _col(49) "Number of groups" `q' "=" /// in ye _col(70) %9.0g e(N_g) } local skip = max(43 - length("`e(sampleratio)'"),1) noi di in green "Sample ratio:" in ye _col(`skip') /// "`e(sampleratio)'" /// in gr _col(49) "Obs per group: min" `q' "=" in ye _col(70) /// %9.0g e(g_min) local skip = max(43 - length("indep"),1) noi di in green "Correlation (case model):" in ye _col(`skip') /// "indep" in gr _col(64) "avg" `q' "=" in ye _col(70) /// %9.1f e(g_avg) local skip = max(43 - length("`e(corr)'"),1) noi di in green "Correlation (`e(depvar)' model):" in ye _col(`skip') /// "`e(corr)'" in gr _col(64) "max" `q' "=" /// in ye _col(70) %9.0g e(g_max) noi di "" noi di in gr _col(26) /// "(Semi-robust Std. Err. adjusted for clustering on `ivar')" ereturn display end