Bug reporting

[kai0511]

Dear authors,

I found the variance ratio test has a bug. In the line, (kurtosis(Y1) - 1)/N1 + (kurtosis(Y0) - 1)/N0 may produce a negative value, so if sqrt a negative value, a Nan will be yielded out.

The code to reproduce the error:

n <- 500
p <- 10
X <- matrix(rnorm(n * p), n, p)
W <- rbinom(n, 1, 0.5)
Y <- pmax(X[, 1], 0) * W + X[, 2] + pmin(X[, 3], 0) + rnorm(n)
variance.ratio.test(Y, W)

I am wondering whether we can abs it before sqrt it.

[bfifield]

Thanks for this, @kai0511. @lmiratrix, looking back at equation A.1 (https://arxiv.org/pdf/1412.5000.pdf) I guess this situation could arise when the kurtosis of Y_1 and/or Y_0 is less than 1? Maybe it's best to just throw a warning in this scenario that the variance ratio test will be invalid, but you might have other thoughts.

[bfifield]

Hi @kai0511, hiding my previous comment because it was incorrect given that kurtosis measures can never fall below 1 (kurtosis \geq skew^2 + 1) As a result, the denominator of the variance ratio test statistic can never be negative, so I'm not certain what's causing the error yet.

Quick question - I tried running your code to reproduce the error, and it ran without a problem. Did I have to set a seed to get it to error out? If you're able to get it to reproduce, and can re-send, I think I can debug pretty quickly.

[kai0511]

@bfifield I found the problem is due to different definition of kurtosis. See the page https://rviews.rstudio.com/2018/01/04/introduction-to-kurtosis/. Subtracting 3 is contained in the definition of kurtosis in the page. I used this definition.

Could you please also shed light on the different definition of kurtosis? Thanks!

[bfifield]

My understanding of the implementation in the moments package, which is what we use here, is that the kurtosis() function is just the standard Pearson's kurtosis - the fourth standardized moment - whereas what you've posted in the link is a related measure called "excess kurtosis". Looking at the raw code in the kurtosis() function, it also looks like they've implemented Pearson's kurtosis rather than excess kurtosis: n * sum((x - mean(x))^4)/(sum((x - mean(x))^2)^2)

The second equation in this section https://en.wikipedia.org/wiki/Kurtosis#Pearson_moments shows how the lowest possible bound of any kurtosis measure is 1, which is when the skew of the random variable Y is equal to 0. If the skew of the random variable is not equal to 0, then the lower bound of the kurtosis measure will be greater than 1.

That's all to say, if you're able to send along a fully reproducible example (the version you posted above ran without throwing an error for me), that would be a big help.

[lmiratrix]

Thanks for handling this. I can drag Peng into it if you think it is useful? He was the primary architect for the variance test.

  --> Due to my RSI (wrist trouble), e-mail often abrupt <--

On Thu, Sep 3, 2020 at 1:10 AM Ben Fifield notifications@github.com wrote:

My understanding of the implementation in the moments package, which is what we use here, is that the kurtosis() function is just the standard Pearson's kurtosis - the fourth standardized moment - whereas what you've posted in the link is a related measure called "excess kurtosis". Looking at the raw code in the kurtosis() function, it also looks like they've implemented Pearson's kurtosis rather than excess kurtosis: n * sum((x - mean(x))^4)/(sum((x - mean(x))^2)^2)

The second equation in this section https://en.wikipedia.org/wiki/Kurtosis#Pearson_moments shows how the lowest possible bound of any kurtosis measure is 1, which is when the skew of the random variable Y is equal to 0. If the skew of the random variable is greater than 0, then the lower bound of the kurtosis measure will also be greater than 1.

That's all to say, if you're able to send along a fully reproducible example (the version you posted above ran without throwing an error for me), that would be a big help.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/bfifield/hettx/issues/19#issuecomment-686256539, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAI63ZOKGRC2FOAYUKA2E5DSD4QNLANCNFSM4PWLVMAQ .

[kai0511]

@lmiratrix @bfifield Thanks for looking at this problem! Here is the code for kurtosis I wrote. kurtosis <- function(x){ m <- mean(x) v <- var(x) return((sum((x - m)^4) / length(x)) / v^2 - 3) }

I used the above kurtosis function in the following test. n <- 500 p <- 10 X <- matrix(rnorm(n p), n, p) W <- rbinom(n, 1, 0.5) Y <- pmax(X[, 1], 0) W + X[, 2] + pmin(X[, 3], 0) + rnorm(n) variance.ratio.test(Y, W)

If I remove the -3 in kurtosis, then the test will be okay. To my understanding, I think either implementation of kurtosis should be enable to use here. You may also find another problem. In the example, I incorporated heterogeneity into the data, but the test cannot detect it.

Hope you can help me with this issue! Thanks!

[bfifield]

Hi @kai0511, thanks for the follow-up. However, I think I'm still a bit confused about the issue you're seeing.

First, in the above example, which version of variance.ratio.test() are you using? Is this the version in the hettx package, or did you write your own version that uses the definition of kurtosis you provided above (which, again, I believe is excess kurtosis rather than the standard kurtosis measure)?

When I run the above example, using the variance.ratio.test() in hettx, I get the following:

> library(hettx)
> set.seed(20200907)
> n <- 500
> p <- 10
> X <- matrix(rnorm(n * p), n, p)
> W <- rbinom(n, 1, 0.5)
> Y <- pmax(X[, 1], 0) * W + X[, 2] + pmin(X[, 3], 0) + rnorm(n)
> variance.ratio.test(Y, W)
      pvalue     var1     var0    ratio log.ratio   asy.se        z
1 0.09802064 2.616641 2.228032 1.174418 0.1607724 0.124349 1.292913

Are you getting something different? Also, could you please provide replication code that includes (1) a set.seed for replicability and (2) the output that shows the exact error you're seeing?

[kai0511]

@bfifield I do not use the kurtosis from the package moments. I use kurtosis with excess, which was coded by myself. To my understanding, I think either moment or excess version of kurtosis should work well. It's not the case here.

[bfifield]

Hi @kai0511 - OK, that's helpful. So the variance.ratio.test works as coded for you (with the standard kurtosis measure), but when you use a hand-coded version of variance.ratio.test that uses an excess kurtosis measure instead, you get NAs. Is that correct?

I coded up what I believe is the problem you're seeing, can you confirm that this is an accurate representation of the function you built and the problem you're seeing?

set.seed(807787)

## Create functions
excess.kurtosis <- function(x){
    m <- mean(x)
    v <- var(x)
    return((sum((x - m)^4) / length(x)) / v^2 - 3)
}

variance.ratio.test.cond <- function(Yobs, Z, data= NULL){
    if (!is.null( data ) ) {
        Yobs = eval( substitute( Yobs ), data )
        Z = eval( substitute( Z ), data )
    }
    Y1 = Yobs[Z==1]
    Y0 = Yobs[Z==0]

    N1 = length(Y1)
    N0 = length(Y0)

    log.varR = log(var(Y1)/var(Y0))

    asy.se   = sqrt(  (excess.kurtosis(Y1) - 1)/N1 + (excess.kurtosis(Y0) - 1)/N0   )

    pvalue   = as.numeric((1 - pnorm(abs(log.varR), 0, asy.se)))

    res = data.frame( pvalue = pvalue, var1 = var( Y1 ), var0 = var( Y0 ))
    res$ratio = res$var1 / res$var0
    res$log.ratio = log( res$ratio )
    res$asy.se = asy.se
    res$z = res$log.ratio / res$asy.se
    return( res )
}

## Test
n <- 500
p <- 10
X <- matrix(rnorm(n * p), n, p)
W <- rbinom(n, 1, 0.5)
Y <- pmax(X[, 1], 0) * W + X[, 2] + pmin(X[, 3], 0) + rnorm(n)
variance.ratio.test.cond(Y, W)

If this is the case, I'll follow up with Luke and Peng about whether excess kurtosis can be used in the function. However, it would be helpful to know (a) why you want to use excess kurtosis instead of standard kurtosis - is there a use case you have in mind where it would be useful? and (b) why you think it should also work here.