wt_var <- function(x, w, na.rm = FALSE) {
  if (na.rm) {
    na <- is.na(x) | is.na(w)
    x <- x[!na]
    w <- w[!na]
  }
  wm <- weighted.mean(x, w)
  sqrdev <- (x - wm)^2
  (sum(sqrdev*w))
}