/*****************************************************************************************************************************************************
*
* TITLE: Meta-analysis of MDQ for bipolar disorder
* AUTHOR: Yemisi Takwoingi
* DATE CREATED: 22/06/2020
* DATE MODIFIED: 13/02/2023
* PURPOSE: Fit the bivariate model to estimate a summary point using only studies that used cut-off 7 and use all data to estimate a summary curve
*
*
* NOTE: Comments begin with * or a block of text begins with /* and closes with */
* Stata is case sensitive so if you are not familiar with it beware when you create variables and type commands and program names.
* Commands are expected to be lowercase. Also be careful with "=" and "==". For example, after the if command, Stata expects "=="
* for a test of equality; "=" produces an error in this case.
*
*****************************************************************************************************************************************************/
*** You only need to run the following 2 lines if gllamm and metandi are not installed on your computer or not up to date ***
ssc install gllamm, replace
ssc install metandi, replace
*** Set your working directory to the appropriate drive where you saved the file MDQ all.csv ***
*** Replace "C:\Users\takwoiny\Desktop\Chapter 16 datasets and code files\Datasets" with the path where the dataset is saved ***
cd "C:\Users\takwoiny\Desktop\Chapter 16 datasets and code files\Datasets"
*** Read in the data from the .csv file.
insheet using "MDQ all.csv", comma clear
*** Produce a summary of the dataset to check data import was ok ***
describe
********************** 2. FITTING A SUMMARY POINT USING METANDI *************************
*** Use metandi for meta-analysis of MDQ to obtain a summary point of sensitivity and specificity using only the data for cutoff of 7
metandi tp fp fn tn if cutoff==7
*** To obtain a SROC plot as well as parameter estimates, add the plot option to the metandi statement as follows ***
metandi tp fp fn tn if cutoff==7, plot
*** There is no option with metandi to modify the plot but this can be done using metandiplot ***
metandiplot
*** If the optional variables tp fp fn tn are included in the command line, the plot also includes estimates of sensitivity and specificity from each study. ***
metandiplot tp fp fn tn if cutoff==7
/* The default is to scale the plot symbol by the sample size of each study.
To make the symbols all the same size, specify constant weights, e.g. [aw=1]. Try other options too, see help file for metandiplot. */
metandiplot tp fp fn tn [aw=1] if cutoff==7, curve(off) pred(off)
*** Here’s another example including some twoway graph options ***
metandiplot tp fp fn tn [aw=1] if cutoff==7, curve(off) legend(off) title(SROC plot of MDQ for BD) scheme(s1mono)
********************** 3. FITTING A SUMMARY CURVE USING METANDI *************************
*** Using the relationship between the bivariate and HSROC models when there is no covariate in the model, use metandi for meta-analysis of MDQ to obtain a summary curve using all available data
metandi tp fp fn tn, plot
*** A summary point is not clinically meaningful given the mix of thresholds (cutoffs) so show only the summary curve ***
*** Remove the summary point and confidence and prediction regions, and colour the curve black
metandiplot tp fp fn tn, summopts(off) pred(off) conf(off) curveopts(lcolor(black))
*** Here's another example including twoway graph options ***
metandiplot tp fp fn tn, summopts(off) pred(off) conf(off) curveopts(lcolor(black)) title(SROC plot of MDQ for BD) scheme(s1mono)
**************** 4. FITTING THE BIVARIATE MODEL WITH MEQRLOGIT (OR XTMELOGIT) ******************
/*** SET UP THE DATA BEFORE THE META-ANALYSIS
Generate 5 new variables of type long. We need these before we can reshape the data.
• n1 is number diseased
• n0 is number without disease
• true1 is number of true positives
• true0 is the number of true negatives
• study is the unique identifier for each study. _n will generate a sequence of numbers.
IMPORTANT NOTE: after modifying the data, if you need to run metandi again, remember to use the original data. ***/
gen long n1=tp+fn
gen long n0=fp+tn
gen long true1=tp
gen long true0=tn
gen long study= _n
*** Convert data from wide to long form ***
reshape long n true, i(study) j(sens)
*** Generate a new binary variable spec of type byte that takes the value 0 when sens=1 and vice versa ***
gen byte spec=1-sens
*** Sort data to ensure studies are clustered together first by study ***
sort studyid
/*** PERFORM META-ANALYSIS: replace meqrlogit with xtmelogit if your version of Stata is earlier than Stata 13
Note that xtmelogit was introduced in Stata 10 and meqrlogit in Stata 13 so you cannot run meqrlogit or xtmelogit if your version is earlier than Stata 10
Explanation of model specification:
• The variable true specifies the response while sens and spec are the fixed portions of the model similar to if we were using regress or some other Stata estimation command. Our fixed
effects are coefficients on sens and spec without a constant term (nocons)
• With || studyid: the random effects were specified at the level identified by the group variable studyid.
• intpoints(5): the number of integration points for adaptive Gaussian quadrature
• cov(): covariance specifies the structure of the covariance matrix for the random effects. cov(un)specifies unstructured covariance allowing all variances and covariances to be distinct.
• nocons: suppresses the constant (intercept) term and is specified here for both the fixed effects and random-effects equations.
• refineopts(iterate(3)): controls the maximization process during the refinement of starting values. Two iterations is the default. Should the maximization fail because of instability in
the Hessian calculations, one possible solution may be to increase the number of iterations here.
• binomial(n): specifies the data are in binomial form and n as the binomial variable.
• variance: displays the random-effects parameter estimates as variances and covariances. To display standard deviations and correlation replace option variance with stddeviations. ***/
meqrlogit true sens spec if cutoff==7, nocons|| studyid: sens spec, ///
nocons cov(un) binomial(n) refineopts(iterate(5)) intpoints(5) variance
/* To find the covariance between the expected (mean) logit sensitivity and expected logit specificity, display contents of the variance-covariance matrix. */
matrix list e(V)
*** DISPLAY SUMMARY ESTIMATES ***
*** Drop the program in case it is already in Stata’s memory. ***
capture program drop renamematrix
/* Write a little program renaming elements of the coefficient and variance matrices. */
program define renamematrix, eclass
matrix mb = e(b)
matrix mv = e(V)
matrix colnames mb = logitse:_cons logitsp:_cons
matrix colnames mv = logitse:_cons logitsp:_cons
matrix rownames mv = logitse:_cons logitsp:_cons
ereturn post mb mv
end
*** Run the program ***
renamematrix
*** Display summary points ***
_diparm logitse, label(Sensitivity) invlogit
_diparm logitsp, label(Specificity) invlogit
_diparm logitse logitsp, label(LR+) ci(log) f(invlogit(@1)/(1-invlogit(@2))) d( exp(@2-1)*invlogit(@1)^2/invlogit(@2) exp(@2)*invlogit(@1))
_diparm logitse logitsp, label(LR-) ci(log) f((1-invlogit(@1))/invlogit(@2)) d( exp(-@1)*invlogit(@1)^2/invlogit(@2) exp(-@1-@2)*invlogit(@1))
_diparm logitse logitsp, label(DOR) ci(log) f(exp(@1+@2)) d(exp(@1+@2) exp(@1+@2))