Study with negative treatment arm variance

[Vanessa-WV]

Hi Guido,

Thanks for your great contribution to the R Community!

I am trying to run a network meta-analysis using SMD as effect size, since I'm new to nma, I selected part of my three six-arm studies dataset as first attempt, but I am struggling with the negative variance error. My data and code are as follow:

Study	treat1	treat2	TE	seTE
Batch 7	a	Control	-1.3178118	0.2898016
Batch 9	b	Control	-0.8340301	0.2207003
Batch 10	b	c	1.3208092	0.2110693
Batch 10	b	Control	-0.780556	0.1811051
Batch 10	c	Control	-2.0122344	0.2879232

try(netmeta(TE = TE,
            seTE = seTE,
            treat1 = treat1,
            treat2 = treat2,
            studlab = Study,
            data = es_table_nma,
            sm = "SMD",
            method.tau = "REML",
            fixed = FALSE,
            random = TRUE,
            reference.group = "Control",
            details.chkmultiarm = TRUE,
            sep.trts = " vs "
            ))

And received the following error:

Warning: Matrix inversion resulted in negative variances for multi-arm study 'Batch 10'. 
Consider changing argument 'func.inverse', e.g., ginv() from MASS package.
Multi-arm studies with inconsistent treatment effects:

  studlab      treat1    treat2        TE          resid
 Batch 10      c         Control      -2.012234    0.02971026
 Batch 10      c         b            -1.320809    -0.02971026
 Batch 10      Control   b            0.780556     0.02971026

Multi-arm studies with negative treatment arm variance:

  studlab       treat    var.treat
 Batch 10        c            0.047325504
 Batch 10        Control      0.035574289
 Batch 10        b            -0.002775241

Legend:
 resid - residual deviation (observed minus expected)
 TE - treatment estimate
 var.treat - treatment arm variance

Error : Problems in multi-arm studies!
  - Study with inconsistent treatment estimates: 'Batch 10'
  - Study with negative treatment arm variance: 'Batch 10'
  - Please check original data used as input to netmeta().
  - Argument 'tol.multiarm' in netmeta() can be used to relax consistency
    assumption for multi-arm studies (if appropriate).

I tried to add tol.multiarm, tol.multiarm.se and func.inverse = MASS::ginv after looking at the fix for #6 and #9 , the inconsistent treatment estimates error was fixed while the negative variance error is still there.

try(netmeta(TE = TE,
            seTE = seTE,
            treat1 = treat1,
            treat2 = treat2,
            studlab = Study,
            data = es_table_nma,
            sm = "SMD",
            method.tau = "REML",
            fixed = FALSE,
            random = TRUE,
            reference.group = "Control",
            details.chkmultiarm = TRUE,
            sep.trts = " vs ",
            tol.multiarm = 0.5,
            tol.multiarm.se = 0.8,
            func.inverse = MASS::ginv
            ))

Warning: Matrix inversion resulted in negative variances for multi-arm study 'Batch 10'. 
Consider changing argument 'func.inverse', e.g., ginv() from MASS package.
Multi-arm studies with negative treatment arm variance:

  studlab       treat    var.treat
 Batch 10        c            0.047325504
 Batch 10        Control      0.035574289
 Batch 10        b            -0.002775241

Legend:
 var.treat - treatment arm variance

Error: Problem in multi-arm studies!
- Study with negative treatment arm variance: 'Batch 10'
- Please check original data used as input to netmeta().

Is there a way to solve this problem? (FYI, I'm using netmeta 2.9-0 and R 4.2.2 on Windows)

Look forward to your reply and your help on this is really appreciated!

Many thanks~

[gertaruecker]

Hi Vanessa, Guido forwarded your e-mail to me. I could reproduce all your error messages and tried the same actions you had tried. I could not solve the problem, because it is in the data. Here a bit of background. In a three-arm study, all treatment effects and also the variances should be consistent. For the treatment effects, this means that b - c = (b - Control) - (c - Control) But: 1.3208092 is not equal to -0.780556 - (-2.0122344) = 1.231678. You solved this by granting a tolerance of 0.5 (0.03 would be enough).

For the variances, we must have (by the rules for the variance of a difference of independent variables): var(b - c) = var(b) + var(c) var(b - Control) = var(b) + var(Control) var(c - Control) = var(c) + var(Control) This is a linear equation which can be solved for the arm-based variances var(b), var(c), var(Control). But it turns out that one of them (var(b)) is negative, which is impossible. This means that the contrast-based variances are inconsistent. Apparently, the extent of inconsistency is too large to be fixed by setting a tolerance. The question is were this inconsistency comes from. I cannot answer this without more knowledge of the data. Could you check how the data were derived and whether they are correct?

Best, Gerta

[Vanessa-WV]

Hi Gerta,

Thank you so much for the reply, it's quite helpful :) At least now I know how the negative variance error came from and why the tolerance adjustment was not working~

To give you more information about the data, I calculated TE and seTE based on the within-group SMD and its corresponding standard error being introduced in section 3.3.1.3 of this guideline Doing Meta-Analysis in R. As all of the treatments within each study were actually applied to a same group of people (i.e. different zones on the arm), which means they're paired samples instead of independent ones.

Therefore, the way I calculated TE is by using the function in the effectsize package.

TE<- hedges_g(Score ~ treatment, 
              data = data %>% filter(treatment %in% comparison)  # where comparison is each contrast i.e. c("a", "b")
              conf.level = 0.95, paired=TRUE)

And calculated seTE based on the formula (3.23) in the guideline as follow

$$SE{SMD{within}} = \sqrt{\frac{2(1-r{t{1}t{2}})}{n} + \frac{SMD^{2}{within}}{2n}}$$

where $r_{t{1}t{2}}$ is the correlation between the assessment points, the related code is as follow.

seTE<- sqrt(((2*(1-cor(data %>% filter(treatment== comparison[1]) %>% .$Score,
                       data %>% filter(treatment== comparison[2]) %>% .$Score, 
                       method = "spearman")))/n1) + (TE$Hedges_g^2/(2*n)))

Then I combined the results with the corresponding treatment pair and obtained the final data as shown in the table above.

Please feel to let me know if there's anything unclear, or you need further information!

I also have one more naive question, how did you get that 0.03 saying it would be enough for the tolerance of inconsistent treatment effects?

Thanks again for your help, I really appreciate it!

Regards, Vanessa

[gertaruecker]

Hi Vanessa, Though I did not go through the code in detail, the reason for inconsistency of both TE and seTE might by the use of Hedges‘ g. It is known that Hedges‘ g does not produce consistent effect estimates in NMA, see the help page to function pairwise() in R package netmeta. Therefore pairwise() uses Cohen’s d instead. Perhaps it might be worth a trial with Cohen’s d to calculate your data for Batch 10 before conducting NMA?

Best, Gerta

[gertaruecker]

Hi Vanessa, I forgot to answer your other question. I used the information of the (absolute) residual for the treatment effects which is given as 0.02971026. Rounding this provides 0.03. Best, Gerta

[Vanessa-WV]

Hi Gerta,

Thanks for your answer! I tried to use Cohen's d for both TE and seTE, while the negative treatment arm variance situation still exists, as the new data didn't differ from the Hedges' g results that much.

Study | treat1 | treat2 | TE_Hedges | seTE_Hedges | TE_Cohen | seTE_Cohen -- | -- | -- | -- | -- | -- | -- Batch 7 | a | Control | -1.3178118 | 0.2898016 | -1.3519425 | 0.2923263 Batch 9 | b | Control | -0.8340301 | 0.2207003 | -0.8581275| 0.2223438 Batch 10 | b | c | 1.3208092 | 0.2110693 | 1.3508658| 0.2137696 Batch 10 | b | Control | -0.780556 | 0.1811051 | -0.7983184 | 0.1822078 Batch 10 | c | Control | -2.0122344 | 0.2879232 | -2.0580251 | 0.2925105

Would calculate TE and seTE based on the paired samples method be a problem when conducting NMA?

Regards, Vanessa

[guido-s]

Vanessa,

I assume that the inconsistency comes from using different standard deviations to calculate SMD_within for the three pairwise comparisons in Batch 10, i.e., equation (3.20) or (3.21) in section 3.3.1.3, Doing Meta-Analysis in R.

I would suggest that you calculate a pooled standard deviation and use this in equation (3.20) or (3.21) for Cohen's d.

[gertaruecker]

I agree with Guido. A pooled standard deviation could be obtained by a simple generalization of eq. (3.16), see text immediately before (3.20).

[Vanessa-WV]

Hi Guido and Gerta,

Thanks for your suggestion! I tried to use both equations (3.20 and 3.21) to calculate TE instead of directly using the build-in function in the package, while the error still exists, the SMD and its standard error didn't change that much just as our previous attempt. Please kindly find the detailed code and the corresponding results below:

    n1<- nrow(data %>% filter(treatment == comparison[1]))
    n2<- nrow(data %>% filter(treatment == comparison[2]))  # n1=n2

    m1<- mean(data %>% filter(treatment == comparison[1]) %>% .$Score)
    m2<- mean(data %>% filter(treatment == comparison[2]) %>% .$Score)

    s1<- sd(data %>% filter(treatment == comparison[1]) %>% .$Score)
    s2<- sd(data %>% filter(treatment == comparison[2]) %>% .$Score)

    s_pooled <- sqrt((((n1-1)*s1^2) + ((n2-1)*s2^2))/((n1-1)+(n2-1)))    

    TE<- (m1-m2)/s_pooled
    # or TE<- (m1-m2)/s2

    seTE<- sqrt(((2*(1-cor(data %>% filter(treatment == comparison[1]) %>% .$Score,
                           data %>% filter(treatment == comparison[2]) %>% .$Score, method = "spearman")))/n1) +
                          (TE^2/(2*n1)))

Study | treat1 | treat2 | TE_spooled | seTE_spooled | TE_s2 | seTE_s2 -- | -- | -- | -- | -- | -- | -- Batch 7 | a | Control | -1.8305545 | 0.3316958 | -1.5100314 | 0.3045524 Batch 9 | b | Control | -0.8826485 | 0.2240514 | -0.9224535 | 0.2268966 Batch 10 | b | c | 1.2433751 | 0.2042395 | 1.2132827 | 0.2016374 Batch 10 | b | Control | -0.6659560 | 0.1744446 | -0.6699755 | 0.1746643 Batch 10 | c | Control | -1.8992382 | 0.2767413 | -1.9611543 | 0.2828433

FYI, the correlation between b and c, b and Control, and c and Control is 0.6565047, 0.5783336, and 0.5615255 respectively.

Thanks, Vanessa

[guido-s]

Vanessa,

You have to calculate one pooled SD for a multi-arm study.

R function pairwise() uses the following function to calculate the pooled SD:

pooled_sd <- function(sd, n) {
  sel <- !is.na(sd) & !is.na(n)
  #
  if (any(sel))
    res <- sqrt(sum((n[sel] - 1) * sd[sel]^2) / sum(n[sel] - 1))
  else
    res <- NA
  #
  res
}

Here, arguments 'sd' and 'n' should be vectors, e.g., of length 3 for a three-arm study.

Concerning the correlations, I could be that you also have to calculate a pooled correlation to get consistent variances.

Best, Guido

[Vanessa-WV]

Hi Guido,

Thanks for pointing that out! After using one pooled SD (based on three treatments) to calculate SMD for each pairwise comparison, the inconsistent problem in treatment estimates was solved even without setting the tol.multiarm argument!!

However, the negative variance error still exists. Therefore, I tried to calculate the pooled correlation to see if this can solve the problem. I searched for the way of doing it and follow the instruction: first, transform the correlation coefficient toward normality to Fisher's $z$ using equation (11.1), then calculate the pooled correlation $z$, finally back transform the result to obtain the pooled correlation using equation (11.4). Then I used the pooled correaltion to calcualte seTE, the negative variance error is still here. The related code is as follow:

# calcualte a pooled sd
study_sd<- study %>% group_by(treatment) %>% dplyr::summarise(sd = sd(Score), n = n())
s_pooled_overall<- sqrt(sum((study_sd$n-1)*study_sd$sd^2)/(sum(study_sd$n)-nrow(study_sd)))

# also need to calcualte a pooled correlation coefficient
cor_pair<- cor(study[,c(1,2,4)] %>% pivot_wider(names_from = treatment, values_from = Score) %>% .[,-1],
               method = "spearman") 

cor_unique<- unique(as.numeric(cor_pair))[-1]  ## extract the correlation result for each pairwise comparison
cor_fisher<- 0.5*(log((1+cor_unique)/(1-cor_unique))) ## Fishers Z transformation
fisher_pooled<- sum(cor_fisher)/length(cor_unique)  ## calculate pooled z using equal weighted average
cor_pooled<- (exp(2*fisher_pooled)-1)/(exp(2*fisher_pooled)+1)  ## back transform

TE<- (m1-m2)/s_pooled_overall
seTE<- sqrt(((2*(1-cor_pooled))/n1) + (TE^2/(2*n1)))

The pooled correlation for Batch 10 is 0.6004729 and the correspond TE and seTE results are as below:

Study | treat1 | treat2 | TE | seTE -- | -- | -- | -- | -- Batch 7 | a | Control | -1.8305545 | 0.3316958 Batch 9 | b | Control | -0.8826485| 0.2240514 Batch 10 | b | c | 1.2587150| 0.2132226 Batch 10 | b | Control | -0.6531304| 0.1700709 Batch 10 | c | Control | -1.9118454| 0.2739463

Regards, Vanessa

[guido-s]

Hi Vanessa,

I talked about this with Gerta. It seems that the linear equation system for the arm-based variances (which Gerta described in her first post) does not hold for paired data. Looks like this is (another) methodological issue with multi-arm studies.

The argument 'tol.multiarm.se' only works for inconsistent variances, however, does not help with a negative arm-based variance.

Instead, you could add or subtract a small increment to / from one of the three standard errors in order to push the negative arm-based variance across zero. E.g., subtracting 0.0013 from the third standard error worked for me.

> library("netmeta")
> load("batch10.rda")
> data.frame(study, TE, seTE, t1, t2)
     study         TE      seTE t1      t2
1 Batch 10  1.2587150 0.2132226  b       c
2 Batch 10 -0.6531304 0.1700709  b Control
3 Batch 10 -1.9118454 0.2739463  c Control
> try(netmeta(TE, seTE, t1, t2, study, sm = "SMD", details.chkmultiarm = TRUE))

Multi-arm studies with negative treatment arm variance:

  studlab   treat     var.treat
 Batch 10       b -0.0003292936
 Batch 10       c  0.0457931707
 Batch 10 Control  0.0292534046

Legend:
 var.treat - treatment arm variance

Error : Problem in multi-arm studies!
  - Study with negative treatment arm variance: 'Batch 10'
  - Please check original data used as input to netmeta().

In addition: Warning message:
Matrix inversion resulted in negative variances for multi-arm study 'Batch 10'. 
Consider changing argument 'func.inverse', e.g., ginv() from MASS package. 
> try(netmeta(TE, c(seTE[1:2], seTE[3] - 0.001), t1, t2, study, sm = "SMD",
+  details.chkmultiarm = TRUE))

Multi-arm studies with negative treatment arm variance:

  studlab   treat     var.treat
 Batch 10       b -5.584725e-05
 Batch 10       c  4.551972e-02
 Batch 10 Control  2.897996e-02

Legend:
 var.treat - treatment arm variance

Error : Problem in multi-arm studies!
  - Study with negative treatment arm variance: 'Batch 10'
  - Please check original data used as input to netmeta().

> try(netmeta(TE, c(seTE[1:2], seTE[3] - 0.0012), t1, t2, study, sm = "SMD",
+   details.chkmultiarm = TRUE))

Multi-arm studies with negative treatment arm variance:

  studlab   treat     var.treat
 Batch 10       b -1.277993e-06
 Batch 10       c  4.546516e-02
 Batch 10 Control  2.892539e-02

Legend:
 var.treat - treatment arm variance

Error : Problem in multi-arm studies!
  - Study with negative treatment arm variance: 'Batch 10'
  - Please check original data used as input to netmeta().

> netmeta(TE, c(seTE[1:2], seTE[3] - 0.0013), t1, t2, study, sm = "SMD",
+   random = FALSE, ref = "cont")
Number of studies: k = 1
Number of pairwise comparisons: m = 3
Number of treatments: n = 3
Number of designs: d = 1

Common effects model

Treatment estimate (sm = 'SMD', comparison: other treatments vs 'Control'):
            SMD             95%-CI     z  p-value
b       -0.6531 [-0.9865; -0.3198] -3.84   0.0001
c       -1.9118 [-2.4462; -1.3775] -7.01 < 0.0001
Control       .                  .     .        .

Quantifying heterogeneity:
tau^2 = NA; tau = NA

Test of heterogeneity:
 Q d.f. p-value
 0    0      --
> netmeta(TE, c(seTE[1:2], seTE[3] - 0.002), t1, t2, study, sm = "SMD",
+   random = FALSE, ref = "cont")
Number of studies: k = 1
Number of pairwise comparisons: m = 3
Number of treatments: n = 3
Number of designs: d = 1

Common effects model

Treatment estimate (sm = 'SMD', comparison: other treatments vs 'Control'):
            SMD             95%-CI     z  p-value
b       -0.6531 [-0.9865; -0.3198] -3.84   0.0001
c       -1.9118 [-2.4449; -1.3788] -7.03 < 0.0001
Control       .                  .     .        .

Quantifying heterogeneity:
tau^2 = NA; tau = NA

Test of heterogeneity:
 Q d.f. p-value
 0    0      --
> try(netmeta(TE, c(seTE[1] + 0.0015, seTE[2:3]), t1, t2, study, sm = "SMD", details.chkmultiarm = TRUE))

Multi-arm studies with negative treatment arm variance:

  studlab   treat     var.treat
 Batch 10       b -8.334653e-06
 Batch 10       c  4.611413e-02
 Batch 10 Control  2.893245e-02

Legend:
 var.treat - treatment arm variance

Error : Problem in multi-arm studies!
  - Study with negative treatment arm variance: 'Batch 10'
  - Please check original data used as input to netmeta().

> netmeta(TE, c(seTE[1] + 0.0016, seTE[2:3]), t1, t2, study, sm = "SMD",
+   random = FALSE, ref = "cont")
Number of studies: k = 1
Number of pairwise comparisons: m = 3
Number of treatments: n = 3
Number of designs: d = 1

Common effects model

Treatment estimate (sm = 'SMD', comparison: other treatments vs 'Control'):
            SMD             95%-CI     z  p-value
b       -0.6531 [-0.9865; -0.3198] -3.84   0.0001
c       -1.9118 [-2.4488; -1.3749] -6.98 < 0.0001
Control       .                  .     .        .

Quantifying heterogeneity:
tau^2 = NA; tau = NA

Test of heterogeneity:
 Q d.f. p-value
 0    0      --

An alternative approach would be to ignore that the three results come from the same observations which would result in narrower confidence intervals (I am not sure whether this is an appropriate approach).

> netmeta(TE, seTE, t1, t2, sm = "SMD", random = FALSE, ref = "cont")
Number of studies: k = 3
Number of pairwise comparisons: m = 3
Number of treatments: n = 3
Number of designs: d = 3

Common effects model

Treatment estimate (sm = 'SMD', comparison: other treatments vs 'Control'):
            SMD             95%-CI     z  p-value
b       -0.6531 [-0.9525; -0.3538] -4.28 < 0.0001
c       -1.9118 [-2.2907; -1.5330] -9.89 < 0.0001
Control       .                  .     .        .

Quantifying heterogeneity / inconsistency:
tau^2 = 0; tau = 0; I^2 = 0%

Tests of heterogeneity (within designs) and inconsistency (between designs):
                Q d.f. p-value
Total           0    1  1.0000
Within designs  0    0      --
Between designs 0    1  1.0000
Warning message:
No information given for argument 'studlab'. Assuming that comparisons are from independent studies. 

Best, Guido

P.S. You could use cor2z() and z2cor() from R package meta to (back)transform between correlations and Fisher's z transformed correlations.

[Vanessa-WV]

Hi Guido,

Thanks for your solution! I tried to subtract 0.0013 to the third seTE and it worked as well!!

Therefore, I'm just wondering how am I supposed to which number to add to or subtract from which standard error among a multi-arm study? Is there a rule or it's just we try different combinations and the negative variance error would disappear by chance?

I also have a quick question, why ref = "cont" could work the same as setting ref = "Control", don't we need specify the base-arm treatment name that match the original input data?

Regards, Vanessa

[guido-s]

Hi Vanessa,

I got to the minimal number to add / subtract by trial and error. ;-)

The input for argument 'ref' (or 'reference.group') is compared with the full intervention names. Accordingly, ref = "cont" is matched against the intervention 'Control'. In your example, ref = "c" or ref = "C" would not have worked due to the interventions 'Control' and 'c'.

Best, Guido