qrule mathematic interpolation in quantile estimation in R survey package

Question

I've been having difficulty understanding the code for qrule_math in the survey package. This document here says this about it:

However, this is not what the code seems to do. Given two vectors of values (x) and weights (w) and a quantile (p) of 0.5 to find the median value:

x <- structure(c(577450.791535772, 919590.385520716, 755016.909933021, 
1515808.3277814), .Dim = c(1L, 4L), .Dimnames = list("ing_t_p_R", 
    c("11110", "31250", "113113", "113118")))
w <- c(`11110` = 311.909858216766, `31250` = 59.2529837088652, `113113` = 34.1247483726328, 
`113118` = 69.4310625894021)
p <- 0.5

#First define helper function last():
last <- function (a) 
   {
     if (any(a)) 
       max(which(a))
     else 1
   }

Function qs finds an upper and a lower quantile:

    qs <- function (x, w, p) 
   {
     n <- length(x)
     ii <- order(x)
     x <- x[ii]
     cumw <- cumsum(w[ii])
     pos <- last(cumw <= p * sum(w))
     posnext <- if (pos == length(x)) 
       pos
     else pos + 1
     list(qlow = x[pos], qup = x[posnext], ilow = pos, iup = posnext, 
          wlow = p - cumw[pos]/sum(w), wup = cumw[posnext]/sum(w) - 
            p)
   }

Then, qrule_math has a condition where it outputs the lower quantile only if wlow is equal to zero:

qrule_math <- function (x, w, p) {
      
      qdata <- qs(x, w, p)
            
      if (qdata$wlow == 0) 
        qdata$qlow
      else qdata$qup
    }   

So applying qrule_math to given x, w and p you get:

qrule_math(x,w,p)
[1] 755016.9

It gets you the upper quantile. This is different from output given by Stata and calculation "by hand" (577451).

So, my question is:

1. Is this output correct?
2. What does it mean the condition given inside qrule_math function, where it gives the lower quantile only if wlow is equal to zero? Shouldn't this function give always the quantile that satisfies to be <= p?

Answer 1

Good catch! This looks to me like a bug. Inside qrule_math(), the condition should be if (qdata$wlow <= 0) rather than ==.

The == makes sense in the "typical" situation where p >= min(w)/sum(w), so that wlow is non-negative. Then qrule_math() should always return the upper quantile qup unless you're exactly on the lower quantile.

But it misses the possibility that p might be smaller than the smallest w's proportion of sum(w).

In your case, the weight w corresponding to the smallest x is over half of sum(w), so you get a negative wlow. In that situation, qrule_math() should be returning min(x)... but because it (mistakenly) assumed wlow will always be non-negative, it only checks for wlow==0 and incorrectly returns the 2nd-smallest x.

---

Simplified example, based on your example data:

x <- c(5,9,7,15)
w <- c(300, 60, 30, 70)
p <- 0.5

# Reorder to make it easier to think about
w <- w[order(x)]
x <- x[order(x)]

x
## [1]  5  7  9 15
cumsum(w)
## [1] 300 330 390 460
cumsum(w)/sum(w)
## [1] 0.6521739 0.7173913 0.8478261 1.0000000

So for instance you'd expect:

• p=0.5 should return x[1] which is 5
• p=0.7 should return x[2] which is 7
• p=sum(w[1:2])/sum(w) should also return x[2] which is 7
• p=0.72 should return x[3] which is 9

BUT the current qrule_math() wrongly returns 7 for the first of these.

qrule_math(x,w,.5)
## [1] 7
qrule_math(x,w,.7)
## [1] 7
qrule_math(x,w,sum(w[1:2]/sum(w)))
## [1] 7
qrule_math(x,w,.72)
## [1] 9

However, if we fix == to <=, it works:

qrule_math_fixed <- function (x, w, p) {
  qdata <- qs(x, w, p)
  if (qdata$wlow <= 0) 
    qdata$qlow
  else qdata$qup
} 

qrule_math_fixed(x,w,.5)
## [1] 5
qrule_math_fixed(x,w,.7)
## [1] 7
qrule_math_fixed(x,w,sum(w[1:2]/sum(w)))
## [1] 7
qrule_math_fixed(x,w,.72)
## [1] 9

Besides qrule_math() aka qrule_hf1(), there seems to be a related issue in qrule_school() aka qrule_hf2(), and in qrule_hf3(). I'll submit a bug report.