Meta-analysis of prevalence rates in R

Question

Currently I am working on a meta-analysis of prevalence rates. There are a lot of packages for meta-analysis in R (i.e. meta, metafor, MAd, ...). Nevertheless, I can't find an appropriate command for prevalence rates. Do you have any experiences with meta-analysis and prevalence rates in R? Which command would you use?

library(metafor)
res <- rma(p(logit), SE(logit), data=data, method="DL")
forest(res)

Any ideas?

Best regards Totti

Answer 1

Here my code I use the control arm of the Olkin95 data which are included in the meta package (only the firs 20 studies). Subgroups were defined cutting the years in three groups of equal size. I decided to put just the results of random-effects model in the plot (comb.fixed = FALSE)

data(Olkin95)
meta <- metaprop(event = event.c,
             n = n.c, studlab = paste(author, year),
             byvar = cut(year, 3),
             data = Olkin95[1:20,])
forest(meta, comb.fixed = FALSE, 
   bylab = "Years subgroup", 
   hetlab = "", print.tau2 = FALSE,
   layout = "RevMan",
   col.square = "black",
   col.square.lines = "black")

In the meta package forest plots could be rendered in different layouts, included the one used by the Cochrane Collaboration (layout = "RevMan"). The forest function use the grid graphics system, which is great but has also its cons. I found impossible to put different forest plots in one panel with par(mfrow()), so I needed to use the grid.grab function to capture a forest plot and grid.arrange to put the forest plots in one panel. A skeleton of the code is given in what follows:

f1 <- metaprop()
forest(f1)
p1 <- grid.grab()

f2 <- metaprop()
forest(f2)
p2 <- grid.grab()

grid.newpage()
grid.arrange(p1, p2, ncol = 1)

Answer 2

The command metaprop of the R package meta can perform meta-analyses of proportions. It has been implemented also in Stata.
A good tutorial on its use is given here. This work shows how to use the Stata version of metaprop, but many things are also true for the R version.
Below I will show the metanalysis of one of the datasets considered in the paper: prevalence of HPV-infection in women with a Pap smear showing ASC-US.

Here is the dataset:

df <- structure(list(study = c("Manos, 1999", "Bergeron, 2000", "Lytwyn, 2000", 
"Shlay, 2000", "Morin, 2001", "Rebello, 2001", "Solomon, 2001", 
"Zielinski, 2001", "Kulasingam, 2002", "Pretorius, 2002", "Cuzick, 2003", 
"Guyot, 2003", "Lonky, 2003", "Wensveen, 2003", "Bruner, 2004", 
"Rowe, 2004", "Andersson, 2005", "Dalla Palma, 2005", "Giovannelli, 2005", 
"Kendall, 2005", "Nieh, 2005", "Bergeron, 2006", "Kelly, 2006", 
"Kiatpongsan, 2006", "Ko, 2006", "Monsonego, 2006", "Moss, 2006", 
"Selvaggi, 2006", "Wright, 2006", "Cuschieri, 2007", "Ronco, 2007", 
"You, 2007"), author = c("Manos", "Bergeron", "Lytwyn", "Shlay", 
"Morin", "Rebello", "Solomon", "Zielinski", "Kulasingam", "Pretorius", 
"Cuzick", "Guyot", "Lonky", "Wensveen", "Bruner", "Rowe", "Andersson", 
"Palma", "Giovannelli", "Kendall", "Nieh", "Bergeron", "Kelly", 
"Kiatpongsan", "Ko", "Monsonego", "Moss", "Selvaggi", "Wright", 
"Cuschieri", "Ronco", "You"), year = c(1999L, 2000L, 2000L, 2000L, 
2001L, 2001L, 2001L, 2001L, 2002L, 2002L, 2003L, 2003L, 2003L, 
2003L, 2004L, 2004L, 2005L, 2005L, 2005L, 2005L, 2005L, 2006L, 
2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2006L, 2007L, 2007L, 
2007L), tgroup = structure(c(1L, 1L, 1L, 1L, 1L, 9L, 1L, 9L, 
1L, 1L, 9L, 9L, 1L, 1L, 6L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 6L, 1L, 
6L, 1L, 9L, 6L, 6L, 9L, 1L, 6L), .Label = c("ASCUS", "LSIL", 
"ASCUS or SIL", "ASC-R", "ASC-R/US", "ASC-US", "ASC-H", "AGC/AGUS", 
"BORDERLINE DYSKARYOSIS"), class = "factor"), num = c(384, 48, 
23, 61, 105, 31, 1306, 74, 138, 306, 32, 12, 128, 67, 25, 121, 
23, 109, 21, 2501, 49, 835, 37, 35, 930, 34, 1680, 266, 438, 
115, 238, 542), denom = c(973, 111, 57, 195, 360, 75, 2301, 213, 
270, 949, 123, 23, 278, 148, 93, 275, 52, 156, 92, 7334, 66, 
1880, 51, 90, 2319, 71, 3681, 672, 1285, 190, 757, 1171), frac = c(0.394699990749359, 
0.432399988174438, 0.403499990701675, 0.312799990177155, 0.291700005531311, 
0.413300007581711, 0.567600011825562, 0.347400009632111, 0.511099994182587, 
0.322400003671646, 0.260199993848801, 0.521700024604797, 0.460399985313416, 
0.452699989080429, 0.268799990415573, 0.439999997615814, 0.442299991846085, 
0.69870001077652, 0.228300005197525, 0.340999990701675, 0.742399990558624, 
0.444099992513657, 0.725499987602234, 0.388900011777878, 0.400999993085861, 
0.478899985551834, 0.456400007009506, 0.395799994468689, 0.340856045484543, 
0.605300009250641, 0.314399987459183, 0.462900012731552), se = c(0.0156999994069338, 
0.0469999983906746, 0.0649999976158142, 0.0331999994814396, 0.0240000002086163, 
0.0568999983370304, 0.0103000001981854, 0.032600000500679, 0.0304000005125999, 
0.0152000002563, 0.0395999997854233, 0.10419999808073, 0.029899999499321, 
0.0408999994397163, 0.046000000089407, 0.029899999499321, 0.0688999965786934, 
0.0366999991238117, 0.0438000001013279, 0.00549999997019768, 
0.0538000017404556, 0.0115000000223517, 0.0625, 0.0513999983668327, 
0.0102000003680587, 0.0593000017106533, 0.00820000004023314, 
0.0188999995589256, 0.0132228191941977, 0.0355000011622906, 0.0168999992311001, 
0.0146000003442168), up = c(0.42616218328476, 0.529843688011169, 
0.54178661108017, 0.382952570915222, 0.341589659452438, 0.532972872257233, 
0.587941765785217, 0.415492951869965, 0.572178602218628, 0.353226482868195, 
0.346979528665543, 0.731803834438324, 0.520975351333618, 0.536525011062622, 
0.370758354663849, 0.500865280628204, 0.586724460124969, 0.769499897956848, 
0.327510416507721, 0.351993471384048, 0.842233598232269, 0.46694752573967, 
0.841072738170624, 0.497442901134491, 0.421312838792801, 0.600784838199615, 
0.47265362739563, 0.433942914009094, 0.367502719163895, 0.675257563591003, 
0.348808646202087, 0.491900950670242), lo = c(0.363788783550262, 
0.338717311620712, 0.275612711906433, 0.248484954237938, 0.245208755135536, 
0.300753623247147, 0.547044515609741, 0.283657312393188, 0.449797093868256, 
0.292769432067871, 0.185249075293541, 0.305878013372421, 0.40074035525322, 
0.370811879634857, 0.182117596268654, 0.380441725254059, 0.304695636034012, 
0.62021142244339, 0.147192806005478, 0.330161929130554, 0.619938731193542, 
0.421524852514267, 0.582552492618561, 0.287862300872803, 0.381007432937622, 
0.358779191970825, 0.440211117267609, 0.358642756938934, 0.314939588308334, 
0.531930029392242, 0.281442880630493, 0.433990597724915)), datalabel = "", time.stamp = "15 Jan 2014 14:08", .Names = c("study", 
"author", "year", "tgroup", "num", "denom", "frac", "se", "up", 
"lo"), formats = c("%20s", "%14s", "%8.0g", "%22.0g", "%9.0g", 
"%9.0g", "%9.0g", "%9.0g", "%9.0g", "%9.0g"), types = c(20L, 
14L, 252L, 251L, 254L, 254L, 254L, 254L, 254L, 254L), val.labels = c("", 
"", "", "tgroup", "", "", "", "", "", ""), var.labels = c("", 
"", "", "", "", "", "", "", "Upper limit", "Lower limit"), row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", 
"25", "26", "27", "28", "29", "30", "31", "32"), version = 12L, label.table = structure(list(
    tgroup = structure(c(1L, 2L, 3L, 11L, 12L, 13L, 14L, 15L, 
    16L), .Names = c("ASCUS", "LSIL", "ASCUS or SIL", "ASC-R", 
    "ASC-R/US", "ASC-US", "ASC-H", "AGC/AGUS", "BORDERLINE DYSKARYOSIS"
    ))), .Names = "tgroup"), class = "data.frame")

and here the code for the meta-analysis:

library(meta)
(mtprop <- metaprop(event=num, n=denom, studlab=study, byvar=droplevels(tgroup), data=df))

The command metaprop yields the following results:

                  proportion           95%-CI %W(fixed) %W(random) byvar
Manos, 1999           0.3947 [0.3638; 0.4262]       3.8        3.7     1
Bergeron, 2000        0.4324 [0.3387; 0.5298]       0.4        2.9     1
Lytwyn, 2000          0.4035 [0.2756; 0.5418]       0.2        2.4     1
Shlay, 2000           0.3128 [0.2485; 0.3830]       0.7        3.2     1
Morin, 2001           0.2917 [0.2452; 0.3416]       1.2        3.4     1
Rebello, 2001         0.4133 [0.3008; 0.5330]       0.3        2.6     3
Solomon, 2001         0.5676 [0.5470; 0.5879]       9.1        3.7     1
Zielinski, 2001       0.3474 [0.2837; 0.4155]       0.8        3.3     3
Kulasingam, 2002      0.5111 [0.4498; 0.5722]       1.1        3.4     1
Pretorius, 2002       0.3224 [0.2928; 0.3532]       3.4        3.7     1
Cuzick, 2003          0.2602 [0.1852; 0.3470]       0.4        2.8     3
Guyot, 2003           0.5217 [0.3059; 0.7318]       0.1        1.6     3
Lonky, 2003           0.4604 [0.4007; 0.5210]       1.1        3.4     1
Wensveen, 2003        0.4527 [0.3708; 0.5365]       0.6        3.1     1
Bruner, 2004          0.2688 [0.1821; 0.3708]       0.3        2.6     2
Rowe, 2004            0.4400 [0.3804; 0.5009]       1.1        3.4     1
Andersson, 2005       0.4423 [0.3047; 0.5867]       0.2        2.3     1
Dalla Palma, 2005     0.6987 [0.6202; 0.7695]       0.5        3.1     1
Giovannelli, 2005     0.2283 [0.1472; 0.3275]       0.3        2.5     1
Kendall, 2005         0.3410 [0.3302; 0.3520]      26.7        3.8     1
Nieh, 2005            0.7424 [0.6199; 0.8422]       0.2        2.3     1
Bergeron, 2006        0.4441 [0.4215; 0.4669]       7.5        3.7     1
Kelly, 2006           0.7255 [0.5826; 0.8411]       0.2        2.1     2
Kiatpongsan, 2006     0.3889 [0.2879; 0.4974]       0.3        2.8     1
Ko, 2006              0.4010 [0.3810; 0.4213]       9.0        3.7     2
Monsonego, 2006       0.4789 [0.3588; 0.6008]       0.3        2.6     1
Moss, 2006            0.4564 [0.4402; 0.4727]      14.8        3.8     3
Selvaggi, 2006        0.3958 [0.3586; 0.4339]       2.6        3.6     2
Wright, 2006          0.3409 [0.3149; 0.3675]       4.7        3.7     2
Cuschieri, 2007       0.6053 [0.5319; 0.6753]       0.7        3.2     3
Ronco, 2007           0.3144 [0.2814; 0.3488]       2.6        3.6     1
You, 2007             0.4629 [0.4340; 0.4919]       4.7        3.7     2

Number of studies combined: k = 32

                     proportion           95%-CI  z  p-value
Fixed effect model       0.4084 [0.4024; 0.4144] --       --
Random effects model     0.4249 [0.3921; 0.4583] --       --

Quantifying heterogeneity:
 tau^2 = 0.1262; H = 4.79 [4.33; 5.31]; I^2 = 95.7% [94.7%; 96.5%]

Test of heterogeneity:
      Q d.f.  p-value
 712.74   31 < 0.0001

Results for subgroups (fixed effect model):
                                 k proportion           95%-CI      Q    tau^2   I^2
byvar = ASCUS                   20     0.3990 [0.3914; 0.4067] 562.93   0.1879 96.6%
byvar = ASC-US                   6     0.4007 [0.3879; 0.4137]  63.59   0.062  92.1%
byvar = BORDERLINE DYSKARYOSIS   6     0.4524 [0.4375; 0.4674]  45.10   0.1551 88.9%

Test for subgroup differences (fixed effect model):
                    Q d.f.  p-value
Between groups  41.13    2 < 0.0001
Within groups  671.62   29 < 0.0001

Results for subgroups (random effects model):
                                 k proportion           95%-CI      Q    tau^2   I^2
byvar = ASCUS                   20     0.4281 [0.3793; 0.4783] 562.93   0.1879 96.6%
byvar = ASC-US                   6     0.4106 [0.3574; 0.4660]  63.59   0.062  92.1%
byvar = BORDERLINE DYSKARYOSIS   6     0.4271 [0.3434; 0.5151]  45.10   0.1551 88.9%

Test for subgroup differences (random effects model):
                    Q d.f.  p-value
Between groups   0.24    2   0.8891

Details on meta-analytical method:
- Inverse variance method
- DerSimonian-Laird estimator for tau^2
- Logit transformation
- Clopper-Pearson confidence interval for individual studies

And the forest plot is:

forest(mtprop)