version 9 mata function covariances(string scalar betas, string scalar doses, string scalar ses, string scalar tot, string scalar cases, real scalar study) { real colvector b, dose, se, n, A, V, cx , ex, sx, sx_1, fitcontrols, SE real scalar i, m1, N, n0, api, a0i, SINCR, DCA, NUMCOR real matrix H1 , H2, AP, INCR, A1, C, X, PASSA, PASSN, R st_view(b=., ., betas, .) st_view(dose=., ., doses, .) st_view(se=., ., ses, .) st_view(n=., ., tot, .) st_view(A=., ., cases, .) V = se:*se // get some number of cases and total obs m1 = colsum(A) N = colsum(n) n0 = n[1] // Loop until converence (no more changes of fitted cases) SINCR = 1 while (abs(SINCR)>1e-5) { api = colsum(A[|2,1 \ ., .|]) a0i = m1 - api ex = J(rows(A), 1,.) if (study==1) { cx = (1:/A) + 1:/(n:-A) for (i=2; i<=rows(A); i++) { ex[i] = b[i] + log(a0i)+log(n[i]-A[i])-log(A[i])-log(n0-a0i) } } if (study==2|study==3) { cx = (1:/A) for (i=2; i<=rows(A); i++) { ex[i] = b[i] + log(m1-api)+log(n[i])-log(A[i])-log(n0) } } H1 = J(rows(A)-1,rows(A)-1,cx[1]) H2 = diag(cx) H = H1 + H2[|2,2\. , .|] AP = A[|2,1 \ ., .|] + invsym(H)*ex[|2,1 \ ., .|] INCR = invsym(H)*ex[|2,1 \ ., .|] SINCR = colsum(INCR) DCA = m1-sum(AP) // New colvector of fitted cases A = (DCA \ AP) } // Covariances and correlations fitcontrols = n :- A sx = J(rows(A), 1,.) // case-control if (study==1) { sx_1 = sqrt(1:/A :+ 1:/fitcontrols :+ 1/A[1] :+ 1/fitcontrols[1]) } // cohort data with person-time if (study==2) { sx_1 = sqrt(1:/A :+ 1/A[1]) } // cohort data with person if (study==3) { sx_1 = sqrt(1:/A :- 1:/n :+ 1/A[1] :- 1/n0) } sx = sx_1[|2,1 \ ., .|] C = diag(V[|2,1 \ ., .|]) SE = diagonal(sqrt(C)) if (study==1) { NUMCOR = (1/A[1] + 1/fitcontrols[1]) } if (study==2) { NUMCOR = (1/A[1]) } if (study==3) { NUMCOR = (1/A[1] - 1/n0) } // fill in the variance-covariance matrix for (j =1; j<=rows(C); j++) { for (i =1; i<=rows(C); i++) { if (i != j) C[i,j] = ( NUMCOR / ( sx[j]*sx[i]) ) * (SE[j]*SE[i]) } } // fill in the correlation matrix R = J(rows(C), cols(C),1) for (j =1; j<=rows(R); j++) { for (i=1; i<=rows(R); i++) { if (i != j) R[i,j] = ( NUMCOR / ( sx[j]*sx[i]) ) // if (i != j) R[i,j] = ( C[i,j] / (SE[j]*SE[i]) ) } } PASSA = A st_matrix("r(fitA)", PASSA) st_matrix("r(C)", C) st_matrix("r(R)", R) } function glsest(string scalar y, string scalar x, string scalar cov) { real matrix Y, X, C, BC, VB real scalar Q, LL Y = st_matrix(y) X = st_matrix(x) C = st_matrix(cov) // Variance and betas VB = invsym(X'*invsym(C)*X) BC = invsym(X'*invsym(C)*X)*X'*invsym(C)*Y // Q goodness of fit Q = (Y-X*BC)'*invsym(C)*(Y-X*BC) C2 = diag(C) // log-likelihood LL = -.5*rows(Y)*log(2*pi())-.5*log(det(C))-.5*Q /* printf("1 %g\n", -.5*rows(Y)*log(2*pi()) ) printf("2 %g\n",-.5*log(det(C)) ) printf("det(C) %g\n", det(C) ) printf("3 %g\n",-.5*Q ) printf("4 %g\n", LL ) */ /* alternative LL = -rows(Y)/2*log(2*pi())-.5*log(det(C))-.5*Q LL2 = log( (2*pi())^(-.5*rows(Y))*det(C)^(-.5)*exp(-.5*(Q))) printf("LL2 %g\n", LL2 ) */ st_matrix("r(BC)", BC) st_matrix("r(VB)", VB) st_numscalar("r(Q)", Q) st_numscalar("r(LL)", LL) st_matrix("r(COVAR)", C) } function crudest(string scalar logrr, string scalar cases, string scalar tot, real scalar study) { real colvector n, A, EXPB, B, OT, COVBV real scalar M_DIFF_C_A_RR real matrix C, R st_view(n=., ., tot, .) st_view(A=., ., cases, .) st_view(LRR=., ., logrr, .) RR = exp(LRR) OT = n - A // estimate crude beta coefficients if (study==1) { EXPB = (A:*OT[1]):/(A[1]:*OT) } if (study==2|study==3) { EXPB = (A:*n[1]):/(A[1]:*n) } B = log(EXPB) // estimate covariances of the crude beta coefficients if (study==1) { COVB = (1/A[1] + 1/OT[1]) } if (study==2) { COVB = (1/A[1]) } if (study==3) { COVB = (1/A[1] - 1/n[1]) } COVBV = J(rows(A)-1,rows(A)-1,COVB) COVBV = COVBV - diag(COVBV) // estimate variances of the crude beta coefficients sx = J(rows(A), 1,.) if (study==1) { sx_1 = (1:/A :+ 1:/OT :+ 1/A[1] :+ 1/OT[1] ) } if (study==2) { sx_1 = (1:/A :+ 1/A[1]) } if (study==3) { sx_1 = (1:/A :- 1:/n :+ 1/A[1] :- 1/n[1]) } sx = sx_1[|2,1 \ ., .|] C = diag(sx[|1,1 \ ., .|]) C = C + COVBV // estimate correlation matrix of the crude beta coefficients R = J(rows(C), cols(C),1) for (j =1; j<=rows(R); j++) { for (i=1; i<=rows(R); i++) { // if (i != j) R[i,j] = COVB / sqrt(C[i,j]*C[j,i]) if (i != j) R[i,j] = COVB / sqrt(sx[j]*sx[i]) } } // relative comparison between crude and adjusted relative risks (check assumption 1) // (crude RR - adjusted RR)/ crude RR * 100 DIFF_C_A_RR = (EXPB :- RR):/(EXPB) * 100 // this provides the percentage of the crude relative risks to be added or subtracted, depending on its // positive or negative sign to the crude relative risks to get the adjusted relative risks // calculate the average relative change (from crude to adjusted) // of the relative risks at different exposure level M_DIFF_C_A_RR = mean( abs(DIFF_C_A_RR[|2,1 \ ., .|] ) ) // I'll also get the maximum and minimum relative change (excluding the referent relative risk) st_matrix("r(RD_C_A)", DIFF_C_A_RR) st_numscalar("r(M_RD_C_A)", M_DIFF_C_A_RR) st_matrix("r(CR_C)", C) st_matrix("r(CR_R)", R) st_matrix("r(CR_B)", B) st_matrix("r(CR_EXPB)", EXPB) st_matrix("r(ADJ_EXPB)", RR) st_numscalar("r(RD_max)", max(DIFF_C_A_RR[|2,1\.,1|]) ) st_numscalar("r(RD_min)", min(DIFF_C_A_RR[|2,1\.,1|]) ) } function boundscov(string scalar a, string scalar u) { real matrix AD, UN, DIF, ADL, ADU real colvector VDIF AD = st_matrix(a) UN = st_matrix(u) DIF = AD - UN VDIF = diagonal(DIF) // create initial values or lower and upper var-cov matrix ADL = AD ADU = AD // replace the off-diagonal elements with a lower and upper limits // as described by Berrington and Cox (Biostatistics, 2003) for (j=1; j<=rows(AD); j++) { for (k=1; k<=rows(AD); k++) { // In the unlikely case (mistake) that the adjusted variance is lower than the un-adjusted variance // in such case we force to zero the covariance if (j!=k) { ADL[j, k] = AD[j,k]-sqrt(VDIF[j]*VDIF[k]) if (ADL[j, k] == .) ADL[j, k] = AD[j,k] } if (j!=k) { ADU[j, k] = AD[j,k]+sqrt(VDIF[j]*VDIF[k]) if (ADU[j, k] == .) ADU[j, k] = AD[j,k] } } } st_matrix("r(UB_C)", ADU) st_matrix("r(LB_C)", ADL) } function blockaccum(string scalar a, string scalar b) { real matrix A, B, GA A = st_matrix(a) B = st_matrix(b) GA = blockdiag(A,B) st_matrix("r(CACC)", GA) } function iterativegls(string scalar y, string scalar b, string scalar x, string scalar cov, real scalar q, real scalar n, real scalar p) { real matrix Y, B, X, C, W, YH, D real scalar tau2, diff, lastau2, stopiter Y = st_matrix(y) B = st_matrix(b) X = st_matrix(x) C = st_matrix(cov) // just get the dose variable, the first in the linear predictor XL = X[.,1] // moment-based estimator of tau2 (between-slope variance) // (q, n, p) // ( q - (n-p) ) // ( trace(invsym(C)) - trace( invsym(C)*X*invsym(X'*invsym(C)*X)*X'*invsym(C) ) ) tau2 = ( q - (n-p) ) / ( trace(invsym(C)) - trace( invsym(C)*X*invsym(X'*invsym(C)*X)*X'*invsym(C) ) ) tau2 = max((tau2, 0)) // printf("Starting tau2=%g\n", tau2) // loop until convergence for tau2 diff = tau2 // if there is not heterogenity (tau2=0) there is no need to loop if (diff==0) { // Variance and betas VB = invsym(X'*invsym(C)*X) BC = invsym(X'*invsym(C)*X)*X'*invsym(C)*Y // Q goodness of fit q = (Y-X*BC)'*invsym(C)*(Y-X*BC) } // printf("Initial diff=%g\n", diff) // printf("Initial tau2=%g\n", tau2) i = 1 D = C while (diff>1e-5) { // printf("i=%g\n", i) lastau2 = tau2 // this should be done for each study and then re-create the block diagonal C = D + tau2*XL*XL' // replace to 0 out of the diagonal block - I'm sure this can be done in a better efficient way for (j =1; j<=rows(C); j++) { for (k =1; k<=rows(C); k++) { if (D[k,j]==0) { C[k,j]=0 } } } // Variance and betas VB = invsym(X'*invsym(C)*X) BC = invsym(X'*invsym(C)*X)*X'*invsym(C)*Y // Q goodness of fit q = (Y-X*BC)'*invsym(C)*(Y-X*BC) // tau2 tau2 = ( q - (n-p) ) / ( trace(invsym(C)) - trace( invsym(C)*X*invsym(X'*invsym(C)*X)*X'*invsym(C) ) ) tau2 = max((tau2, 0)) if (tau2 != 0) { diff = abs(tau2 - lastau2) stopiter = 0 } else { diff = 0 stopiter = 1 } // printf("stop=%g\n", stopiter) // printf("diff=%g\n", diff) // printf("tau2=%g\n", tau2) i++ } // if the loop has been stopped because tau2 < 0 then the results will be equal to fixed-effects if (stopiter==1) { // I have to go back to original sigma C = D // Variance and betas VB = invsym(X'*invsym(C)*X) BC = invsym(X'*invsym(C)*X)*X'*invsym(C)*Y // Q goodness of fit q = (Y-X*BC)'*invsym(C)*(Y-X*BC) } //printf("convergence done\n") // (BC, sqrt(VB)) // (tau2, q) // get the log-likelihood LL = -.5*rows(Y)*log(2*pi())-.5*log(det(C))-.5*q // LL = log( (2*pi())^(-.5*rows(Y))*det(C)^(-.5)*exp(-.5*(q))) st_matrix("r(C)", C) st_matrix("r(BC)", BC) st_matrix("r(VB)", VB) st_numscalar("r(Q)", q) st_numscalar("r(tau2)", tau2) st_numscalar("r(LL)", LL) } function checklowhigh(string scalar sl, string scalar su) { real matrix cov_cru, cov_adj, zero, diffca, low, high real scalar nc, meandiffca cov_cru = st_matrix(sl) cov_adj = st_matrix(su) // number of covariances nc = (rows(cov_cru)*rows(cov_cru)-rows(cov_cru)) /2 lw = 0 uw = 0 zero = J(rows(cov_cru), cols(cov_cru),0) low = J(rows(cov_cru), cols(cov_cru),0) high = J(rows(cov_cru), cols(cov_cru),0) // relative differences between crude and adjusted elements diffca = abs((cov_cru:-cov_adj):/cov_cru)*100 // initialize a vector to containing off-diagonal elements corr = J(nc,1,0) s = 1 // check if cov_cru <=> cov_adj for (j=1; j<=rows(cov_cru); j++) { for (k=1; k<=rows(cov_cru); k++) { zero[j,k] = 1 if (j!=k) { if (zero[j,k]!=1 | zero[k,j]!=1 ) { corr[s++,1] = diffca[j,k] if (cov_cru[k,j]=cov_adj[k,j]) { // printf("Crude>=Adjusted") // (j, k, cov_cru[j,k], cov_adj[j,k]) high[k,j]=1 ++uw } } } } } meandiffca = mean(corr[|1,1 \ ., .|] ) st_matrix("r(diffca)", diffca) st_numscalar("r(meandiffca)", meandiffca) st_numscalar("r(nc)", nc) st_numscalar("r(lw)", lw) st_numscalar("r(uw)", uw) st_matrix("r(low)", low) st_matrix("r(high)", high) } end **************************************************************************************************************************** *! N.Orsini v.1.0.0 18may2005 - Version 8 *! N.Orsini v.2.0.0 16aug2005 - Version 9 *! N.Orsini v.3.0.0 28oct2005 - Version 9 *! N.Orsini v.4.0.0 16nov2005 - Version 9 *! N.Orsini v.5.0.0 20dec2005 - Version 9 *! N.Orsini v.6.0.0 22aug2006 - Version 9.2 *! N.Orsini v.7.0.0 9may2009 - Version 9.2 capture program drop glst program glst, eclass version 9.2 if replay() { if "`e(cmd)'"!="glst" { `e(cmd)' `0' error 301 } Display `0' exit } syntax varlist (min=2) [if] [in] , Se(varname numeric) Cov(varlist numeric min=2) /// [ cc ir ci PFirst(varlist numeric min=2) Level(integer $S_level) eform LBCov(string) UBCov(string) CRudes Random vwls TStage(string) SSest ] // check estimation method - random-effects is allowed only with multiple studies if "`pfirst'" == "" & "`random'" != "" { di as err "random-effects model are allowed only with multiple studies, specify the option pfirst()" exit 198 } if "`pfirst'" != "" & "`crudes'" != "" { di as err "crude estimates are provided only for a single study, do not specify the option pfirst()" exit 198 } // check that the two-stage method is specified with multiple studies if "`pfirst'" == "" & "`tstage'" != "" { di as err "specifies the option pfirst()" exit 198 } // check the two-stage method (fixed or random). if "`tstage'" != "" { * di "`tstage'" if inlist("`tstage'", "f", "r") != 1 { di in red "choose either fixed (f) or random (r) effects two-stage meta-analysis with the option tstage() " exit 198 } } if "`random'" != "" & "`tstage'" != "" { di as err "take out the option random when specifying the option tstage()" exit 198 } // check methods for the confidence bounds of the covariances if "`lbcov'"!="" & "`ubcov'"!="" { di in red "choose lbcov() or ubcov() otpion" exit 198 } if "`lbcov'"!="" { confirm integer number `lbcov' if inrange(`lbcov', 1, 3) != 1 { di in red "choose method either 1,2, or 3 for the option lbcov() " exit 198 } else { global lbcov = `lbcov' global lbounds = "yes" global ubcov = 0 } } if "`ubcov'"!="" { capture confirm integer number `ubcov' if inrange(`ubcov', 1, 3) != 1 { di in red "choose method either 1,2, or 3 for the option ubcov() " exit 198 } else { global ubcov = `ubcov' global ubounds = "yes" global lbcov = 0 } } // check study design (single study) if "`pfirst'" == "" { if "`cc'"=="" & "`ir'"=="" & "`ci'"=="" { di in red "choose cc or ir or ci " exit 198 } if "`cc'"!="" & "`ir'"!="" & "`ci'"!="" { di in red "choose cc or ir or ci " exit 198 } if "`ir'"!="" & "`ci'"!="" { di in red "choose ir or ci " exit 198 } if "`cc'"!="" & "`ir'"!="" { di in red "choose cc or ir " exit 198 } if "`cc'"!="" & "`ci'"!="" { di in red "choose cc or ci " exit 198 } if "`ir'"!="" { local ptcohort = "ptcohort" } if "`cc'"!="" { local casecontrol = "casecontrol" } if "`ci'"!="" { local pecohort = "pecohort" } } // check level if `level' < 10 | `level' > 99 { di in red "level() must be between 10 and 99 inclusive" exit 198 } tempname levelci global levelp = `level'/100 scalar `levelci' = `level' * 0.005 + 0.50 global levelci = `level' * 0.005 + 0.50 // to use the option if/in marksample touse, strok novarlist global touse = `touse' // check observation qui count if `touse' local nobs = r(N) if `nobs'== 0 { di in red "no observations" exit 2000 } if `nobs'== 1 { di in red "insufficient number of observations" exit 2001 } // Parsing input variables tempvar b dose n case control v quietly { tokenize `varlist' local depname = "`1'" mac shift local xname `*' local ncov : word count `xname' gen double `v' = `se'^2 if `touse' parse "`varlist'" , parse(" ") gen double `b' = `1' if `touse' gen double `dose' = `2' if `touse' // check number of covariates and number of observations qui count if `touse' local nobstouse = r(N)-1 if `ncov' > `nobstouse' { di in r "number of covariates greater than the number of observations" exit 198 } // get the variables required to if "`cov'"!="" { parse "`cov'" , parse(" ") gen double `n' = `1' if `touse' gen double `case' = `2' if `touse' capture confirm numeric var `n' capture confirm numeric var `case' } if "`pfirst'"!="" { local multiple = "multiple" tempvar id typestudy parse "`pfirst'" , parse(" ") clonevar `id' = `1' if `touse' gen double `typestudy' = `2' if `touse' capture confirm numeric var `id' capture confirm numeric var `typestudy' capture assert (`typestudy' == 1) | (`typestudy' == 2) | (`typestudy' == 3) if `touse' if _rc != 0 { di _n as err "variable `2' cant take only value 1 (cc), 2 (ir), or 3 (ci)" } } } // end quietly // Multiple studies if "`pfirst'" != "" { global loglik 0 local i = 1 if "`tstage'"!= "" { tempvar ltb ltvar qui gen `ltb' = . qui gen `ltvar' = . } * // Warning message when assuming zero-correlations * if "`vwls'" != "" di _n in g "Assuming zero-correlation within each set of log relative risks" qui levelsof `id' , local(studies) qui tab `pfirst' if `touse' local nstudies = r(r) local countstudy = 1 foreach study of local studies { varcov `varlist' if `id' == `study' & `touse', s(`se') tot(`n') c(`case') /// ts(`typestudy') `crudes' `vwls' // required step for the two-stage meta-analysis tempname M_bc M_vc * get beta and variance of the linear trend mat `M_bc' = r(BC) mat `M_vc' = r(VB) scalar sltb`i' = `M_bc'[1,1] scalar sltvar`i' = `M_vc'[1,1] // display study-specific linear trends, std error and 95% CI if "`ssest'" != "" { local fmt = "%9.0g" local labstudy : label (`id') `study' if `countstudy' == 1 { di _n in g "Study-specific linear trend estimates" di "{dup 78:-}" di in g _col(8) "Study" _col(21) "Coef." _col(29) "Std. Err." _col(44) "z" /// _col(49) "P>|z|" _col(59) "[95% Conf. Interval]" di "{dup 13:-}{c TT}{dup 64:-}" di in g _col(5) %8s abbrev("`labstudy'",8) _col(14) "{c |}" in y _col(17) `fmt' sltb`i' _col(28) `fmt' sqrt(sltvar`i') _col(41) %3.2f sltb`i'/sqrt(sltvar`i') /// _col(49) %4.3f [1- norm(abs(sltb`i'/sqrt(sltvar`i') ) )]*2 _col(59) `fmt' sltb`i'-invnorm(`levelci')*sqrt(sltvar`i') `fmt' _col(69) sltb`i'+invnorm(`levelci')*sqrt(sltvar`i') } else { di in g _col(5) %8s abbrev("`labstudy'",8) _col(14) "{c |}" in y _col(17) `fmt' sltb`i' _col(28) `fmt' sqrt(sltvar`i') _col(41) %3.2f sltb`i'/sqrt(sltvar`i') /// _col(49) %4.3f [1- norm(abs(sltb`i'/sqrt(sltvar`i') ) )]*2 _col(59) `fmt' sltb`i'-invnorm(`levelci')*sqrt(sltvar`i') `fmt' _col(69) sltb`i'+invnorm(`levelci')*sqrt(sltvar`i') } if `countstudy' == `nstudies' di in g "{dup 13:-}{c BT}{dup 64:-}" } local countstudy = `countstudy' + 1 mat C = r(C) mat G`i' = C scalar col`i' = colsof(G`i') local i = `i' + 1 } preserve // I could take out this! tempvar ids tempname X L G quietly { sort `pfirst', stable bys `pfirst' : gen `ids' = cond(_n == 1, 0, 1) keep if `ids' == 1 & `touse' mkmat `xname', matrix(`X') mkmat `depname' , matrix(`L') } // stop if random-effects are required with only one study! if "`pfirst'" != "" & "`random'" != "" & `nstudies' == 1 { di as err "random-effects model are allowed with multiple studies" exit 198 } * create a block-diagonal matrix using my blockaccum mata function quietly { local i = 1 forv l = 1(1)`nstudies' { if `i' == 1 { mat PG = G`i' } if `i' != 1 { mata: blockaccum("PG","G`l'") mat PG = r(CACC) } local ++i } mat `G' = PG } // end quietly // Generalized Least Squares Estimation tempvar b tempname V BETAS b COV tau2 Q ll * just change colnames * noisily di _n in w "Call Mata glsest()" mata: glsest("`L'", "`X'", "`G'") mat `G' = r(COVAR) mat `V' = r(VB) mat `b' = r(BC)' scalar `Q' = r(Q) local Q = `Q' local obs = c(N) // Random-effects estimation if "`random'" != "" { * noisily di _n in w "Call Mata iterativegls()" mata: iterativegls("`L'", "`b'","`X'", "`G'", `Q', `obs', `ncov') scalar `tau2' = r(tau2) eret scalar tau2 = r(tau2) } mat `COV' = r(C) mat `b' = r(BC)' mat `V' = r(VB) global loglik = r(LL) eret scalar ll = r(LL) eret scalar chi2_gf = r(Q) // two-stage fixed and random effects meta-analysis if "`tstage'" != "" { * create a variable with beta and its variance forv i = 1(1)`nstudies' { qui replace `ltb' = sltb`i' in `i' qui replace `ltvar' = sltvar`i' in `i' } * two-step method tempvar w hat v resid tobeused constantone qui gen `constantone' = 1 qui gen `tobeused' = `ltb' != . tempname sumw tau2 b V Q df_Q * use method-of-moments qui gen double `w'=1/`ltvar' summ `w' , meanonly scalar `sumw'=r(sum) * di in w "two-stage fixed-effects meta-analysis" * list `idstudy' `ltb' `ltvar' if `tobeused' qui regress `ltb' `constantone' [iw=`w'] if `tobeused' , mse1 noconstant qui predict `hat' if `tobeused', hat qui replace `hat'=`hat'*`w'^2 scalar `df_Q' = `nstudies' - 1 scalar `Q' = e(rss) summ `hat', meanonly scalar `tau2' = max( (`Q' - `df_Q') / ( `sumw' - r(sum) ), 0) if "`tstage'" == "r" { * di in w "two-stage random-effects meta-analysis" qui gen double `v' = `ltvar' + `tau2' qui regress `ltb' `constantone' [iw=1/`v'], mse1 noconstant } matrix `b' = e(b) matrix `V' = e(V) local df_m = e(df_m) local df = `nstudies'-1 } } // end multiple studies // single study if "`pfirst'" == "" { /* check standard error not equal to zero capture assert `se' != 0 if _rc != 0 { di as err "standard errors equal to zero are not allowed" exit 198 } */ varcov `varlist' `if' `in' , s(`se') tot(`n') c(`case') /// `pecohort' `ptcohort' `casecontrol' `crudes' `vwls' global loglik = r(ll) tempname b V COV R local obs = r(N) mat `COV' = r(C) mat `R' = r(R) mat `b' = r(BC)' mat `V' = r(VB) if "`crudes'" != "" { tempname RD_min RD_max CR_C CR_EXPB CR_B RD_C_A M_RD_C_A CR_R ADJ_EXPB MEANDIFF_C_A DIFF_C_A mat `CR_C' = r(CR_C) mat `CR_R' = r(CR_R) mat `CR_EXPB' = r(CR_EXPB) mat `ADJ_EXPB' = r(ADJ_EXPB) mat `CR_B' = r(CR_B) mat `RD_C_A' = r(RD_C_A) scalar `RD_min' = r(RD_min) scalar `RD_max' = r(RD_max) scalar `M_RD_C_A' = r(M_RD_C_A) scalar `MEANDIFF_C_A' = r(MEANDIFF_C_A) mat `DIFF_C_A' = r(DIFF_C_A) } } // Saved results mat rownames `b' = `depname' mat colnames `b' = `xname' mat rownames `V' = `xname' mat colnames `V' = `xname' if "`tstage'" != "" { ereturn post `b' `V' , dep(`depname') obs(`nstudies') esample(`tobeused') ereturn scalar tau2 = `tau2' ereturn scalar df_m = `df_m' ereturn scalar Q = `Q' ereturn scalar df_Q = `df_Q' tempname sp_betas sp_var_betas mkmat `ltb' if `ltb' != . , matrix(`sp_betas') mkmat `ltvar' if `ltb' != . , matrix(`sp_var_betas') mat colnames `sp_betas' = c1 mat colnames `sp_var_betas' = c1 ereturn matrix sp_betas = `sp_betas' ereturn matrix sp_var_betas = `sp_var_betas' } else { local df = `obs' - rowsof(`V') tempvar tousecopy gen `tousecopy' = `touse' ereturn post `b' `V' , dep(`depname') obs(`obs') esample(`touse') } if "`crudes'" != "" { ereturn matrix C_Sigma =`CR_C' ereturn matrix C_R =`CR_R' ereturn matrix C_y =`CR_B' ereturn matrix C_expy =`CR_EXPB' ereturn matrix A_expy =`ADJ_EXPB' ereturn matrix RD_C_A = `RD_C_A' ereturn scalar RD_min = `RD_min' ereturn scalar RD_max = `RD_max' ereturn scalar M_RD_C_A = `M_RD_C_A' ereturn scalar MEANDIFF_C_A = `MEANDIFF_C_A' ereturn matrix DIFF_C_A = `DIFF_C_A' } if "`pfirst'" == "" { ereturn matrix Sigma = `COV' ereturn matrix R = `R' } else { ereturn matrix Sigma = `G' } if "`random'" != "" { eret scalar tau2 = r(tau2) } ereturn scalar df_gf = `df' ereturn local cmd = "glst" ereturn local depvar = "`depname'" ereturn scalar chi2_gf = r(Q) if "`tstage'" != "" ereturn scalar chi2_gf = `Q' // Q-test heterogeneity if `df' > 0 { if "`pfirst'" == "" eret scalar chi2_gf = r(chi2_gf2) } else { eret scalar chi2_gf = . } /* compute test of indepvars = 0 */ qui testparm `xname' eret scalar df_m = r(df) eret scalar chi2 = r(chi2) eret scalar ll = $loglik if "`pfirst'" != "" { eret scalar S = `nstudies' } Display, level(`level') `eform' `random' `multiple' twostage(`tstage') `vwls' // drop global used macro drop lbcov ubcov loglik ubounds lbounds levelci levelp end // Algorithm to fit the variance-covariance matrix capture program drop varcov program varcov, rclass version 8.2 syntax varlist [if] [in] , s(varname numeric) tot(varname numeric) c(varname numeric) /// [CRudes pecohort ptcohort casecontrol ts(varname numeric) vwls ] tempvar b dose se n case control v tokenize `varlist' local depname = "`1'" mac shift local xname `*' // to use the option if/in marksample touse, strok novarlist preserve qui keep if `touse' & `depname' != . // get the arguments parse "`varlist'" , parse(" ") quietly { gen double `b' = `1' gen double `dose' = `2' gen double `se' = `s' gen double `v' = `se'^2 gen double `n' = `tot' gen double `case' = `c' } // end quietly // get the type of study (either single study or multiple study) if "`ts'" != "" { if `ts'[1] == 1 local typestudy = 1 if `ts'[1] == 2 local typestudy = 2 if `ts'[1] == 3 local typestudy = 3 } else { if "`casecontrol'" != "" local typestudy = 1 if "`ptcohort'" != "" local typestudy = 2 if "`pecohort'" != "" local typestudy = 3 } // check exposure referent level if `b'[1] != 0 { di in r "beta coefficient (log relative risk) at the referent exposure level should be zero, log(1)=0" exit 198 } // Call mata function to calculate covariances * noisily di in w _n "Call Mata covariances()" mata: covariances("`b'", "`dose'", "`se'", "`n'", "`case'", `typestudy') tempname C FCASE R mat `C' = r(C) mat `R' = r(R) mat `FCASE' = r(fitA) * return matrix R = `R' tempname VB BC PASS X XP B bc vb pval ll Q fitcases mkmat `xname' if `b' != 0, matrix(`X') mkmat `b' if `b' != 0, matrix(`B') svmat `FCASE', name(`fitcases') // Estimate beta corrected using weighted least squares for correlated outcomes // Generalized Least Squares * noisily di _n in w "Call Mata glsest()" // Zero-covariance options (replaced 0 off-diagonal) if "`vwls'" != "" { * Replace the previous C with a diagonal matrix (0 off-diagonal elements) tempname CZ mat `CZ' = diag(vecdiag(`C')) mat `C' = `CZ' } mata: glsest("`B'", "`X'", "`C'") mat `VB' = r(VB) mat `BC' = r(BC) scalar `Q' = r(Q) scalar `ll' = r(LL) // Get crude (un-adjusted) estimates of the betas and var-cov matrix if "`crudes'" != "" { * noisily di in w _n "Call Mata crudest()" mata: crudest("`b'", "`case'", "`n'", `typestudy') tempname RD_min RD_max CR_C CR_EXPB CR_B RD_C_A M_RD_C_A CR_R ADJ_EXPB MEANDIFF_C_A DIFF_C_A mat `CR_C' = r(CR_C) mat `CR_R' = r(CR_R) mat `CR_EXPB' = r(CR_EXPB) mat `ADJ_EXPB' = r(ADJ_EXPB) mat `CR_B' = r(CR_B) mat `RD_C_A' = r(RD_C_A) scalar `M_RD_C_A' = r(M_RD_C_A) scalar `RD_min' = r(RD_min) scalar `RD_max' = r(RD_max) // Compare crude and adjusted correlation matrix mata: checklowhigh("`CR_R'", "`R'") scalar `MEANDIFF_C_A' = r(meandiffca) mat `DIFF_C_A' = r(diffca) /* mat list `CR_R' mat list `R' di _n "N. Corr = " r(nc) di "N. Lower corr = " r(lw) di "N. Higher corr = "r(uw) mat list r(low) mat list r(high) */ } // Create Confidence Bounds for the var-cov matrix if ( "$lbounds" != "" | "$ubounds" != "") { if ($lbcov == 1) | ($ubcov == 1) { noi di as res "Method 1" // generate CI of RR tempname lb ub gen double `lb' = `b' - `se'*invnorm($levelci) gen double `ub' = `b' + `se'*invnorm($levelci) // Call mata function to calculate covariances if ($lbcov == 1) { mata: covariances("`lb'", "`dose'", "`se'", "`n'", "`case'", `typestudy') mat `C' = r(C) mata: glsest("`B'", "`X'", "`C'") mat `VB' = r(VB) mat `BC' = r(BC) scalar `Q' = r(Q) scalar `ll' = r(LL) return scalar ll = `ll' return matrix BC = `BC' return matrix VB = `VB' return matrix C = `C' ret scalar chi2_gf2 = `Q' } if ($ubcov == 1) { mata: covariances("`ub'", "`dose'", "`se'", "`n'", "`case'", `typestudy') mat `C' = r(C) mata: glsest("`B'", "`X'", "`C'") mat `VB' = r(VB) mat `BC' = r(BC) scalar `Q' = r(Q) scalar `ll' = r(LL) return scalar ll = `ll' return matrix BC = `BC' return matrix VB = `VB' return matrix C = `C' ret scalar chi2_gf2 = `Q' } } if ($lbcov == 2) | ($ubcov == 2) { noi di as res "Method 2" // find CI for the cases and N using exact poisson // I assume that it's not possible in this setting to have zero cases tempvar lbfitc ubfitc lbn ubn gen `lbfitc' = round(invgammap( `fitcases' , (1-$levelp)/2 ) ) gen `ubfitc' = round(invgammap((`fitcases'+ 1) , (1+$levelp)/2 ) ) gen `lbn' = round(invgammap( `n' , (1-$levelp)/2 )) gen `ubn' = round(invgammap((`n'+ 1) , (1+$levelp)/2 )) // list `lbfitc' `ubfitc' `lbn' `ubn' `case' `n' // call crudest function to get the variance/covariance at lower/upper bounds if ($lbcov == 2) { mata: crudest("`b'", "`ubfitc'", "`ubn'", `typestudy') mat `C' = r(CR_C) mat list `C' mata: glsest("`B'", "`X'", "`C'") mat `VB' = r(VB) mat `BC' = r(BC) scalar `Q' = r(Q) scalar `ll' = r(LL) return scalar ll = `ll' return matrix BC = `BC' return matrix VB = `VB' return matrix C = `C' ret scalar chi2_gf2 = `Q' } if ($ubcov == 2) { mata: crudest("`b'", "`lbfitc'", "`lbn'", `typestudy') mat `C' = r(CR_C) mata: glsest("`B'", "`X'", "`C'") mat `VB' = r(VB) mat `BC' = r(BC) scalar `Q' = r(Q) scalar `ll' = r(LL) return scalar ll = `ll' return matrix BC = `BC' return matrix VB = `VB' return matrix C = `C' ret scalar chi2_gf2 = `Q' } } if ($lbcov == 3) | ($ubcov == 3) { noi di as res "Method 3" // Berrington and Cox method // Difference Adjusted and Un-adjusted (crude) var-cov matrix if $lbcov == 3 | $ubcov == 3 { * noisily di in w _n "Call Mata crudest()" mata: crudest("`b'", "`case'", "`n'", `typestudy') tempname CR_C CR_EXPB CR_B mat `CR_C' = r(CR_C) mat `CR_EXPB' = r(CR_EXPB) mat `CR_B' = r(CR_B) } * noisily di in w _n "Call Mata boundscov()" mata: boundscov("`C'", "`CR_C'") tempname SIGMA_U SIGMA_L mat `SIGMA_U' = r(UB_C) mat `SIGMA_L' = r(LB_C) * noisily di _n in w "Estimates at lower confidence limit of the covariances" if ($lbcov == 3) { mata: glsest("`B'", "`X'", "`SIGMA_L'") tempname VB_L BC_L Q_L ll_L mat `VB_L' = r(VB) mat `BC_L' = r(BC) scalar `Q_L' = r(Q) scalar `ll_L' = r(LL) return scalar ll = `ll_L' return matrix BC = `BC_L' return matrix VB = `VB_L' return matrix C = `SIGMA_L' ret scalar chi2_gf2 = `Q_L' } if ($ubcov == 3) { * noisily di _n in w "Estimates at upper confidence limit of the covariances" mata: glsest("`B'", "`X'", "`SIGMA_U'") tempname VB_U BC_U Q_U ll_U mat `VB_U' = r(VB) mat `BC_U' = r(BC) scalar `Q_U' = r(Q) scalar `ll_U' = r(LL) return scalar ll = `ll_U' return matrix BC = `BC_U' return matrix VB = `VB_U' return matrix C = `SIGMA_U' ret scalar chi2_gf2 = `Q_U' } } } matrix colnames `X' = doses matrix colnames `B' = betas return matrix x = `X' return matrix B = `B' return matrix R = `R' return scalar N = c(N)- 1 if "`crudes'" != "" { return mat CR_C = `CR_C' return mat CR_EXPB = `CR_EXPB' return mat ADJ_EXPB = `ADJ_EXPB' return mat CR_B = `CR_B' return matrix RD_C_A = `RD_C_A' return scalar M_RD_C_A = `M_RD_C_A' return scalar RD_min = `RD_min' return scalar RD_max = `RD_max' return mat CR_R = `CR_R' return scalar MEANDIFF_C_A = `MEANDIFF_C_A' return matrix DIFF_C_A = `DIFF_C_A' } if ("$lbounds" != "yes") & ("$ubounds" != "yes") { * noi di in w "nobounds" // saved results return scalar ll = `ll' return matrix BC = `BC' return matrix VB = `VB' return matrix C = `C' ret scalar chi2_gf2 = `Q' } end capture program drop Display program define Display version 8 syntax [, Level(integer $S_level) eform Random multiple twostage(string) vwls ] if `level' < 10 | `level' > 99 { di in red "level() must be between 10 and 99 inclusive" exit 198 } tempname pgf pm if e(df_gf) > 0 { scalar `pgf' = chiprob(e(df_gf), e(chi2_gf)) } else scalar `pgf' = . scalar `pm' = chiprob(e(df_m), e(chi2)) if "`multiple'" != "" & "`random'" == "" & "`twostage'" == "" { di _n in gr "Fixed-effects dose-response model" _col(50) "Number of studies" _col(70) "= " as res %7.0g e(S) } if "`multiple'" != "" & "`random'" == "" & "`twostage'" == "f" { di _n in gr "Two-stage fixed-effects dose-response model" _col(50) "Number of studies" _col(70) "= " as res %7.0g e(S) } if "`multiple'" != "" & "`random'" != "" & "`twostage'" == "" { di _n in gr "Random-effects dose-response model" _col(50) "Number of studies" _col(70) "= " as res %7.0g e(S) } if "`multiple'" != "" & "`random'" == "" & "`twostage'" == "r" { di _n in gr "Two-stage random-effects dose-response model" _col(50) "Number of studies" _col(70) "= " as res %7.0g e(S) } * if "`multiple'" == "" di _n "" if "`random'" == "" { if "`vwls'" != "" { #delimit ; di _n in gr "Variance-weighted least-squares regression" _col(54) "Number of obs = " in ye %7.0g e(N) _n in gr "Goodness-of-fit chi2(" in ye e(df_gf) in gr ")" _col(28) "= " in ye %7.2f e(chi2_gf) _col(54) in gr "Model chi2(" in ye e(df_m) in gr ")" _col(70) "= " in ye %7.2f e(chi2) _n in gr "Prob > chi2" _col(28) "= " in ye %7.4f `pgf' _col(54) in gr "Prob > chi2 = " in ye %7.4f `pm' ; #delimit cr } else { #delimit ; di _n in gr "Generalized least-squares regression" _col(54) "Number of obs = " in ye %7.0g e(N) _n in gr "Goodness-of-fit chi2(" in ye e(df_gf) in gr ")" _col(28) "= " in ye %7.2f e(chi2_gf) _col(54) in gr "Model chi2(" in ye e(df_m) in gr ")" _col(70) "= " in ye %7.2f e(chi2) _n in gr "Prob > chi2" _col(28) "= " in ye %7.4f `pgf' _col(54) in gr "Prob > chi2 = " in ye %7.4f `pm' ; #delimit cr } } else { if "`vwls'" != "" { #delimit ; di _n in gr "Variance-weighted least-squares regression" _col(54) "Number of obs = " in ye %7.0g e(N) _n in gr "Goodness-of-fit chi2(" in ye e(df_gf) in gr ")" _col(28) "= " in ye %7.2f e(chi2_gf) _col(54) in gr "Model chi2(" in ye e(df_m) in gr ")" _col(70) "= " in ye %7.2f e(chi2) _n in gr "Prob > chi2" _col(28) "= " in ye %7.4f `pgf' _col(54) in gr "Prob > chi2 = " in ye %7.4f `pm' ; #delimit cr } else { #delimit ; di _n in gr "Iterative Generalized least-squares regression" _col(54) "Number of obs = " in ye %7.0g e(N) _n in gr "Goodness-of-fit chi2(" in ye e(df_gf) in gr ")" _col(28) "= " in ye %7.2f e(chi2_gf) _col(54) in gr "Model chi2(" in ye e(df_m) in gr ")" _col(70) "= " in ye %7.2f e(chi2) _n in gr "Prob > chi2" _col(28) "= " in ye %7.4f `pgf' _col(54) in gr "Prob > chi2 = " in ye %7.4f `pm' ; #delimit cr } } if "`eform'" == "" ereturn display, level(`level') else ereturn display, level(`level') eform(exb(b)) if "`random'" != "" { di as res "Moment-based" /* */ as txt " estimate of between-study variance of the slope:" /* */ _col(63) "tau2" _col(66) " = " as res %8.7f e(tau2) } if "`twostage'" == "r" { di as res "Moment-based" /* */ as txt " estimate of between-study variance of the slope:" /* */ _col(63) "tau2" _col(66) " = " as res %8.7f e(tau2) } end