[David Fisher]

Hi Suhail,

Great, thanks -- it sounds like we've solved this one and I just need to roll out a fix.

In answer to your question, user-specified quality weights are always divided by the maximum value, to rescale to 0-1, prior to processing. It's actually only just occurred to me, but in Appendix A of your 2015 article it states that a correction for negative weights is necessary if there exist any Q_j < 1. But in my metan implementation, the Q_j are always re-scaled, as I say ... and therefore the "correction" to Q_dash_j is always done. Is that how you anticipated it?

Thanks,
David.

[Suhail Doi]

Hi David,

Thanks for the update. I just did the analysis using the "qiplus1" column in the dataset above and there is no error. I checked both ways using the rank (qiplus1 divided by 4) as well as the raw qiplus1 count of safeguards variable. Both give the same result and no errors so your observation seems correct. I also checked other datasets and this only happens if there is a zero qwt for one or more studies.

Just to be clear - are all user entered qwt divided by the maximum value prior to processing in metan or do you check the maximum value and only do this if the qwt maximum value is not 1? Either way would work but perhaps just doing this in all cases would be better.


Regards
Suhail

[David Fisher]

Hi Suhail Doi

Having taken a look at your example, I think this is to do with rounding error. With reference to Appendix A of your 2015 article (dx.doi.org/10.1016/j.cct.2015.05.010), we see that tau_hat_j is defined as a function (let's call it "f") of the inverse-variance weights and quality weights, minus the original tau_j. When applied to the last two observations in your example dataset, that function "f" works out to be equal to tau_j, so that the difference tau_hat_j is equal to zero. But in the presence of rounding error, Stata is calculating this to be a very very small negative quantity -- hence the error message "negative weights encountered".

It can easily be shown (and as you may already be aware) that this happens exactly when the "quality weight" is zero -- as it is for the last two observations of your example dataset. This simplifies the expression for the adjusted/corrected quality weight, leading to other terms cancelling out in the expression for tau_hat_j.

So, I should be able to fix this, either by detecting when the quality weight is zero, and/or by being cleverer with the way I arrange the mathematical expressions so that Stata evaluates tau_hat_j as exactly zero rather than a (very) small negative number. In the meantime: is it possible to run the same analysis but with a different set of quality weights which do not take on zero values?

Thanks,

David.


P.S. I am a little surprised that your code ran using the ranks (score/maxscore). On my machine, both your commands failed, and for the same reason. I suppose it's possible that Stata successfully evaluated tau_hat_j as zero on your machine using ranks but not using raw scores. If anything, I might have expected the opposite to be true, because the raw scores are whole numbers and hence less at risk of rounding error. Could you check this again, and confirm your Stata and metan version numbers, and the error messages and codes that you are seeing? Thanks!

[Suhail Doi]

Hi David,

With the below dataset

* Example generated by -dataex-. To install: ssc install dataex
clear
input str15 study int(n1 c1 nc1 n2 c2 nc2) float rank byte(qi qiplus1)
"Sandberg2000"      40  19   21   41  13   28        1 3 4
"Clapp1989"         56   0   56   59   5   54 .6666667 2 3
"Bussel1990"        61  20   41   65  23   42 .6666667 2 3
"Fanaroff1994"    1204 186 1018 1212 209 1003 .6666667 2 3
"Haque1986"        100   4   96   50   5   45 .3333333 1 2
"Weisman1994a"     372  40  332  381  39  342 .3333333 1 2
"Chirico1987"       43   2   41   43   8   35 .3333333 1 2
"Conway1990"        34   8   26   32  14   18 .3333333 1 2
"Ratrisawadi1991"   68  10   58   34  13   21        0 0 1
"Tanzer1997"        40   3   37   40   8   32        0 0 1
end

And the following code with ranks (score/maxscore) gives the correct result

metan c1 nc1 c2 nc2 , or model(qe,qwt(rank))  study( study )  forestplot(astext(85) textsize(100) boxscale(35) spacing(1) leftjustify range(0.3 2) dp(2))  extraline(yes) hetinfo(q pvalue tausq)

However this code using raw scores gives an error

metan c1 nc1 c2 nc2 , or model(qe,qwt(qi) \ re )  study( study )  forestplot(astext(85) textsize(100) boxscale(35) spacing(1) leftjustify range(0.3 2) dp(2))  extraline(yes) hetinfo(q pvalue tausq)

I assumed that metan automatically converted all qwt to ranks by dividing by maximum in the list as even if this is done with a rank input by a user instead of a score there would be no issues as all ranks always start from 1. If this is not the case can this behavior please be updated in the next metan update as may otherwise cause errors for some users? Am I missing something?

Regards
Suhail

[oy lf]

Hi David, that works great, thanks!

Regards

Feng Li

[David Fisher]

Hi again,

Just to say that it occurred to me that the code I used above (for the early data manipulation) could actually be simplified quite a lot! So, here is some alternative code. The ultimate destination, and use of metan and forestplot commands, are identical in both cases.

Data entry (same as above):

* Example generated by -dataex-. For more info, type help dataex
clear
input str26 Disease str23 StudyName str25 Intervention str11 Control int(IntN ControlN) float(HR_pos HRlci_pos HRuci_pos HR_neg HRlci_neg HRuci_neg)
"Non-small-cell lung cancer" "CheckMate 057, 2017"     "Nivolumab"                 "Docetaxel"   231 224 .43  .3  .63 1.01 .76 1.33
"Non-small-cell lung cancer" "CheckMate 017, 2017"     "Nivolumab"                 "Docetaxel"   117 108 .53 .31  .89   .7 .47 1.02
"Non-small-cell lung cancer" "POPLAR, 2016"            "Atezolizumab"              "Docetaxel"   144 143 .54 .33  .89  .84 .53 1.32
"Non-small-cell lung cancer" "OAK, 2016"               "Atezolizumab"              "Docetaxel"   607 608 .64 .49  .84  .85 .73    1
"Melanoma"                   "CheckMate 067, 2018 (a)" "Nivolumab"                 "Ipilimumab"  288 277 .65 .42  .99  .64  .5  .82
"Melanoma"                   "CheckMate 067, 2018 (b)" "Nivolumab plus ipilimumab" "Ipilimumab"  278 277 .56 .35   .9  .54 .42  .69
"Melanoma"                   "CheckMate 069, 2016"     "Nivolumab plus ipilimumab" "Ipilimumab"   35  83 .78 .26  2.4  .74 .38 1.47
"Melanoma"                   "CheckMate 066, 2018"     "Nivolumab"                 "Dacarbazine" 120 243 .16 .07  .38  .56 .37  .85
"Melanoma"                   "CheckMate 037, 2017"     "Nivolumab"                 "CTx"         201 195 .73 .49 1.09 1.15 .82 1.62
"Head and neck cancer"       "CheckMate 141, 2018"     "Nivolumab"                 "SOC"         172 103  .5  .3  .83  .81 .55 1.21
end

Setup prior to meta-analysis:

gen double lnHR_pos = ln(HR_pos)
gen double lnHRlci_pos = ln(HRlci_pos)
gen double lnHRuci_pos = ln(HRuci_pos)

gen double lnHR_neg = ln(HR_neg)
gen double lnHRlci_neg = ln(HRlci_neg)
gen double lnHRuci_neg = ln(HRuci_neg)

gen IntString = Intervention + " (n=" + strofreal(IntN, "%3.0f") + ")"
label variable IntString `""Intervention " "(No. of patients analysed)""'
format %-1s IntString

gen ControlString = Control + " (n=" + strofreal(ControlN, "%3.0f") + ")"
label variable ControlString `""Control" "(No. of patients analysed)""'
format %-1s ControlString

As before, perform two meta-analyses, for PD-L1 positive and negative separately; then append and re-sort:

preserve
    tempfile neg
    metan lnHR_neg lnHRlci_neg lnHRuci_neg, random study(StudyName) by(Disease) nosubgroup nograph clear
    save `neg'
restore
metan lnHR_pos lnHRlci_pos lnHRuci_pos, random study(StudyName) by(Disease) lcols(IntString ControlString) nosubgroup nograph clear
append using `neg', gen(PDL1)

isid _BY _USE _STUDY PDL1, sort miss

A little more data manipulation prior to running forestplot:

replace PDL1 = PDL1 + 1
label define PDL1_ 1 "PD-L1 Positive" 2 "PD-L1 Negative"
label values PDL1 PDL1_

drop if inlist(_USE, 0, 6) & PDL1==1
replace PDL1 = . if inlist(_USE, 0, 6)

replace _LABELS = "Pooled estimate in PD-L1 positive at 5% cut-off" if _USE==5 & PDL1==1
replace _LABELS = "Pooled estimate in PD-L1 negative at 5% cut-off" if _USE==5 & PDL1==2

replace _LABELS = "" if _USE!=5 & PDL1==2
label variable _LABELS "Study name"
format %-1s _LABELS

And finally, creation of the plot (same as before):

forestplot, hr nowt plotid(PDL1) lcols(IntString ControlString) xlabel(.5 "0.5" 1 2) range(.0625 2.5) cirange(.08 2.5) ///
    astext(70) textsize(90) boxsca(80) favours(Favours intervention # Favours control, fp(3)) ///
    box1opts(mcolor(edkblue)) box2opts(mcolor(orange)) ciopts(lcolor(black) lw(thin)) ///
    diam1opts(lc(edkblue)) diam2opts(lc(orange)) oline1opts(lp(dash) lc(edkblue)) oline2opts(lp(dash) lc(orange)) graphregion(color(white))

Thanks,

David.

[David Fisher]

Hi oy lf

This is a great question, thanks!

Since I don't have your specific underlying dataset, I'm going to instead use the PDL1 data in your attachment, which a quick Google informs me is taken from the Supplementary Information (Figure S1) of Liu X et al. Int. J. Cancer 2020; 147: 116–127. doi:10.1002/ijc.32744

The data can be entered into Stata as follows:

* Example generated by -dataex-. For more info, type help dataex
clear
input str26 Disease str23 StudyName str25 Intervention str11 Control int(IntN ControlN) float(HR_pos HRlci_pos HRuci_pos HR_neg HRlci_neg HRuci_neg)
"Non-small-cell lung cancer" "CheckMate 057, 2017"     "Nivolumab"                 "Docetaxel"   231 224 .43  .3  .63 1.01 .76 1.33
"Non-small-cell lung cancer" "CheckMate 017, 2017"     "Nivolumab"                 "Docetaxel"   117 108 .53 .31  .89   .7 .47 1.02
"Non-small-cell lung cancer" "POPLAR, 2016"            "Atezolizumab"              "Docetaxel"   144 143 .54 .33  .89  .84 .53 1.32
"Non-small-cell lung cancer" "OAK, 2016"               "Atezolizumab"              "Docetaxel"   607 608 .64 .49  .84  .85 .73    1
"Melanoma"                   "CheckMate 067, 2018 (a)" "Nivolumab"                 "Ipilimumab"  288 277 .65 .42  .99  .64  .5  .82
"Melanoma"                   "CheckMate 067, 2018 (b)" "Nivolumab plus ipilimumab" "Ipilimumab"  278 277 .56 .35   .9  .54 .42  .69
"Melanoma"                   "CheckMate 069, 2016"     "Nivolumab plus ipilimumab" "Ipilimumab"   35  83 .78 .26  2.4  .74 .38 1.47
"Melanoma"                   "CheckMate 066, 2018"     "Nivolumab"                 "Dacarbazine" 120 243 .16 .07  .38  .56 .37  .85
"Melanoma"                   "CheckMate 037, 2017"     "Nivolumab"                 "CTx"         201 195 .73 .49 1.09 1.15 .82 1.62
"Head and neck cancer"       "CheckMate 141, 2018"     "Nivolumab"                 "SOC"         172 103  .5  .3  .83  .81 .55 1.21
end

We begin by processing the data into a form suitable for meta-analysis; note that the first encode command is simply so that we maintain the same order of studies as in Figure S1 (otherwise they would be re-sorted alphabetically):

label define StudyName_ 1 "CheckMate 057, 2017"
label define StudyName_ 2 "CheckMate 017, 2017", add
label define StudyName_ 3 "POPLAR, 2016", add
label define StudyName_ 4 "OAK, 2016", add
label define StudyName_ 5 "CheckMate 067, 2018 (a)", add
label define StudyName_ 6 "CheckMate 067, 2018 (b)", add
label define StudyName_ 7 "CheckMate 069, 2016", add
label define StudyName_ 8 "CheckMate 066, 2018", add
label define StudyName_ 9 "CheckMate 037, 2017", add
label define StudyName_ 10 "CheckMate 141, 2018", add
rename StudyName StudyNameStr
encode StudyNameStr, gen(StudyName) label(StudyName_)
order StudyName, after(StudyNameStr)
drop StudyNameStr

rename (HR_pos HRlci_pos HRuci_pos) (HR1 HR1_lci HR1_uci)
rename (HR_neg HRlci_neg HRuci_neg) (HR2 HR2_lci HR2_uci)
reshape long HR@ HR@_lci HR@_uci, i(StudyName) j(PDL1)

gen double lnHR = ln(HR)
gen double lnHR_lci = ln(HR_lci)
gen double lnHR_uci = ln(HR_uci)

Now, Figure S1 effectively shows two separate meta-analyses, but presents them as interlaced. So to recreate the figure, we will do the same thing: fit both analyses separately, take the results, and re-sort them so that the PD-L1 positive and negative results alternate.

preserve
    tempfile neg
    metan lnHR lnHR_lci lnHR_uci if PDL1==2, random study(StudyName) by(Disease) lcols(Intervention Control IntN ControlN) nosubgroup nograph clear
    save `neg'
restore
metan lnHR lnHR_lci lnHR_uci if PDL1==1, random study(StudyName) by(Disease) lcols(Intervention Control IntN ControlN) nosubgroup nograph clear
append using `neg', gen(PDL1)

isid _BY _USE _STUDY PDL1, sort miss

Finally, we need to do a bit of tidying up before we create the final plot:

replace PDL1 = PDL1 + 1
label define PDL1_ 1 "PD-L1 Positive" 2 "PD-L1 Negative"
label values PDL1 PDL1_

drop if inlist(_USE, 0, 6) & PDL1==1
replace PDL1 = . if inlist(_USE, 0, 6)

replace _LABELS = "Pooled estimate in PD-L1 positive at 5% cut-off" if _USE==5 & PDL1==1
replace _LABELS = "Pooled estimate in PD-L1 negative at 5% cut-off" if _USE==5 & PDL1==2

replace _LABELS = "" if _USE!=5 & PDL1==2
label variable _LABELS "Study name"
format %-1s _LABELS

gen IntString = Intervention + " (n=" + strofreal(IntN, "%3.0f") + ")" if !missing(IntN) & PDL1==1
label variable IntString `""Intervention " "(No. of patients analysed)""'
format %-1s IntString

gen ControlString = Control + " (n=" + strofreal(ControlN, "%3.0f") + ")" if !missing(ControlN) & PDL1==1
label variable ControlString `""Control" "(No. of patients analysed)""'
format %-1s ControlString

The following line of code will now almost recreate Figure S1:

forestplot, hr nowt plotid(PDL1) lcols(IntString ControlString) xlabel(.5 "0.5" 1 2) range(.0625 2.5) cirange(.08 2.5) ///
    astext(70) textsize(90) boxsca(80) favours(Favours intervention # Favours control, fp(3)) ///
    box1opts(mcolor(edkblue)) box2opts(mcolor(orange)) ciopts(lcolor(black) lw(thin)) ///
    diam1opts(lc(edkblue)) diam2opts(lc(orange)) oline1opts(lp(dash) lc(edkblue)) oline2opts(lp(dash) lc(orange)) graphregion(color(white))

The only thing missing is the second vertical line (the orange line corresponding to PD-L1 negative). Unfortunately, there is currently a bug in the forestplot command which prevents this from displaying. We would like to add the option dataid(PDL1), but as you would see if you tried it, it gives an error message. I have corrected the error in a beta version, and can confirm that the line does appear once the bug is fixed (see attachment). I shall release a new version of the metan package as soon as possible -- sorry about that!


I hope that this gives you some useful pointers.

Best wishes,

David.

[oy lf]

Hi David,

Nice package the updated metan!
I have come across a problem. Hope you can help

How can I get this picture with two pooled results with different colors



I only get the picture below with the code"metan lnhr lnll lnul,eform effect(HR) label(namevar=trial) plotid(line) forestplot(favours(Favours treatment # Favours control) boxopts(msize(*.12)) box1opts(mcolor(orange)) ciopts(lcolor(black) lw(thin)) box2opts(mcolor(edkblue)) diam1opts(lc(orange)) diam2opts(lc(edkblue)) oline1opts(lp(dash) lc(orange)) oline2opts(lp(dash) lc(edkblue)) xlabel(0.2 2.5)) scheme(s1color) "


Thank you very much

[David Fisher]

I am in contact with Mark Ebell via email.

For the benefit other readers with similar issues: when Stata's user-created packages start misbehaving, there are various options. Sometimes, the approach suggested by Suhail Doi in the post immediately above works; that is, to force an install of the latest package, over-writing any conflicting files. The downside is that you cannot be sure exactly which files were conflicting, which were over-written, and which are now in use by Stata.

Therefore, a more comprehensive approach would be to identify all the relevant (user-created!!) packages which are causing issues, uninstall them all (creating a backup if necessary, e.g. if the packages are no longer available for download!), and re-install one by one. Tools for helping here include:

. ado dir                     // list all installed packages, with internal Stata ID numbers in square brackets; optionally restrict to a particular pkgname
. ado dir pkgname             // same but restricted to a particular pkgname
. which pkgname, all          // list all file locations and versions for a particular pkgname
. ado uninstall [x]           // uninstall the package with internal Stata ID number [x]  (including the square brackets, and where x is an integer)

Finally: this shouldn't really matter, but it might be advisable to install metan before installing ipdmetan, rather than the other way around.

I hope that helps!

David.

[Suhail Doi]

try

ssc install metan, replace

[Mark Ebell]

As shown below, I have the old version of metan installed. I have tried both ado update, update and ssc install metan to try and update it with no luck. Any help would be much appreciated.

Thanks, Mark

. which metan
/Users/markebell/Library/Application Support/Stata/ado/plus/m/metan.ado
*! 3.03 19May2009

. ssc uninstall metan
package not found
r(111);

. ssc install metan
checking metan consistency and verifying not already installed...

the following files already exist and are different:
/Users/markebell/Library/Application Support/Stata/ado/plus/m/metan.ado
/Users/markebell/Library/Application Support/Stata/ado/plus/f/forestplot.ado
/Users/markebell/Library/Application Support/Stata/ado/plus/f/forestplot.sthlp

no files installed or copied
(no action taken)
r(602);

[Skye Newton]

Nevermind, I managed to get it to work by using either

“metan loghr logll logul, eform effect("HR") re study(study)” OR
“admetan loghr logll logul, hr re study(study)“

Many thanks anyway!
Kind regards,
Skye

Hi David,

I'm still a bit clueless here, so apologies if I'm getting errors due to my ignorance!

I'm using Stata 15.1, and the answer to 'which metan' and 'which forestplo't are both version 4.03 28apr2021.

I'm trying to do a meta-analysis of hazard ratios from 3 studies using summary measures, so have been using the command:
metan hr ll ul, effect(Hazard Ratio) null(1) xlabel(0 .5 1 1.5) /// lcols(study) texts(200)
This produces the overall HR, but is saying 'Error in forestplot' r(198)

Any suggestions on what I should try to get it produce a forestplot?
Many thanks,
Skye

[Skye Newton]

Hi David,

I'm still a bit clueless here, so apologies if I'm getting errors due to my ignorance!

I'm using Stata 15.1, and the answer to 'which metan' and 'which forestplo't are both version 4.03 28apr2021.

I'm trying to do a meta-analysis of hazard ratios from 3 studies using summary measures, so have been using the command:
metan hr ll ul, effect(Hazard Ratio) null(1) xlabel(0 .5 1 1.5) /// lcols(study) texts(200)
This produces the overall HR, but is saying 'Error in forestplot' r(198)

Any suggestions on what I should try to get it produce a forestplot?
Many thanks,
Skye

[Suhail Doi]

Hi David,

Seems like another bug report:

For this data:

* Example generated by -dataex-. To install: ssc install dataex
clear
input str3 studyname long n int cases byte qi
"S1"  217154 422 1
"S13"    676   1 1
"S18"     44   1 1
"S26"     29   1 1
end

and this code:

metan cases n  , pr model(ivhet \ qe,qwt(qi) \ re \ iv ) transform(ftukey, iv) study( study )  denom(1000)  nograph nohet

the metan results are:

Overall,	IVhet	0.831	0.000	7.773
Overall,	QE	0.831	0.000	7.773
Overall,	DL	2.931	0.034	8.898
Overall,	IV	1.949	1.768	2.138

The results should be (changes in bold font):

Overall,	IVhet	1.951	0.000	7.773
Overall,	QE	1.951	0.000	7.773
Overall,	DL	2.931	0.000	8.898
Overall,	IV	1.949	1.768	2.138

Regards
Suhail

[Suhail Doi]

Hi David, that works great, thanks - seems I missed this in the help file!

Regards
Suhail