Skip to content

[R-meta] A possible bug in clubSandwich::Wald_test()

4 messages · James Pustejovsky, Stefanou Revesz, Farzad Keyhan

Stefanou Revesz
# 
Dear Dr. Pustejovsky (and other Meta Experts),

I've run into a potential bug in clubSandwich::Wald_test(). The error says:

Error: the leading minor of order 6 is not positive definite

Would you please confirm whether the error in the following
reproducible example is a bug or not? Any possible fix?

Best,
Stefanou

dd <- read.csv("https://raw.githubusercontent.com/fpqq/w/main/1.csv")
dd$obs <- 1:nrow(dd)

res <- rma.mv(gi ~ 0 + teaching_level*time, vi, random = ~1|study/obs,
              data = dd)

clubSandwich::Wald_test(res, constrain_zero(res$btt), "CR2")
1 day later
James Pustejovsky
# 
Stefanou,

I don't think this is a bug per se. The error occurs because your
variance-covariance matrix is not positive-definite (and therefore cannot
be inverted for a Wald test). This happens because two of the levels you're
trying to test occur in just a single study each:
+   group_by(teaching_level, time) %>%
+   summarise(
+     effects = n(),
+     studies = n_distinct(study)
+   )
# A tibble: 9 x 4
# Groups:   teaching_level [3]
  teaching_level time        effects studies
  <chr>          <chr>         <int>   <int>
1 elementary     Baseline         11       8
2 elementary     Post-test 1      18      11
3 mixed          Baseline          5       5
4 mixed          Post-test 1       9       8
5 mixed          Post-test 2       1       1
6 mixed          Post-test 3       1       1
7 secondary      Baseline          8       8
8 secondary      Post-test 1      12       9
9 secondary      Post-test 2       3       3

Sandwich estimators will simply not work for such tests.

Apart from this, it's not clear to me why you would want to simultaneously
test all nine coefficients in this model. The null being tested is that the
average effect size is equal to zero in each of the nine groups.

James

On Sun, Mar 13, 2022 at 1:01 PM Stefanou Revesz <stefanourevesz at gmail.com>
wrote:

  
  


  


      

    

    

      
      

        


  

        

      Stefanou Revesz
    

    

      
      #
    

  


  

      
Dear James,

Thank you for your response. To answer your question, my goal was to
get a robust QM test.

Another doubt that I have is that in the model below, `anova.rma(res1,
X = c(-1,1,0,1,-1,0,0,0,0))` gives me the model-based contrast that I
want.

But to get the same contrast in Wald_test(), shall I use
`Wald_test(res, constrain=rbind(c(-1,1,0,1,-1,0,0,0,0)), "CR2")`?

res1 <- rma.mv(gi ~ 0 + teaching_level:time, vi, random = ~1|study/obs,
              data = dd)

anova.rma(res1, X = c(-1,1,0,1,-1,0,0,0,0))
Wald_test(res1, constrain=rbind(c(-1,1,0,1,-1,0,0,0,0)), "CR2")

Many thanks,
Stefanou

        
On Mon, Mar 14, 2022 at 2:09 PM James Pustejovsky <jepusto at gmail.com> wrote:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Farzad Keyhan
    

    

      
      #
    

  


  

      
Dear Stefanou and James,

I also noticed your robust p-values are often much smaller than your
model-based results. I always expected the opposite to happen. But
maybe that is another indication that robust tests won't work well in
such cases.

Fred

On Mon, Mar 14, 2022 at 2:35 PM Stefanou Revesz

        
<stefanourevesz at gmail.com> wrote: