95% CI of CV [General Sta­tis­tics]

posted by mittyri – Russia, 2018-03-21 11:35 (1312 d 16:46 ago) – Posting: # 18576
Views: 1,969

Hi Seppe,

cannot get exactly the same (assuming normal distribution)
confint.cv <- function(x,alpha=.05, method="modmckay"){
  # Calculate the confidence interval of the cv of the vector x
  # Author: Kevin Wright
  # See: Vangel, Mark.  Confidence Intervals for a Normal Coefficient
  # of Variation. American Statistician, Vol 15, No1, p. 21--26.
  # x <- c(326,302,307,299,329)
  # confint.cv(x,.05,"modmckay")
 
  x <- na.omit(x)
  v <- length(x)-1
  mu <- mean(x)
  sigma <- sqrt(var(x))
  k <- sigma/mu
  # CV > .33 may give poor results, so warn the user
  if(k>.33) warning("Confidence interval may be very approximate.")
 
  method <- casefold(method) # In case we see "McKay"
 
  if(method=="mckay"){
    # McKay method.  See equation 15.
    t1 <- qchisq(1-alpha/2,v)/v
    t2 <- qchisq(alpha/2,v)/v
    u1 <- v*t1
    u2 <- v*t2
    lower <- k/sqrt(( u1/(v+1) -1)*k*k + u1/v)
    upper <- k/sqrt(( u2/(v+1) -1)*k*k + u2/v)
  } else if (method=="naive"){
    # Naive method.  See equation 17.
    t1 <- qchisq(1-alpha/2,v)/v
    t2 <- qchisq(alpha/2,v)/v
    lower <- k/sqrt(t1)
    upper <- k/sqrt(t2)
  } else {
    # Modified McKay method. See equation 16.
    method="modmckay"
    u1 <- qchisq(1-alpha/2,v)
    u2 <- qchisq(alpha/2,v)
    lower <- k/sqrt(( (u1+2)/(v+1) -1)*k*k + u1/v)
    upper <- k/sqrt(( (u2+2)/(v+1) -1)*k*k + u2/v)
  }
  ci <- cbind.data.frame(method = method, lower=lower,upper=upper, CV= k, alpha=alpha)
  return(ci)
}
x=c(3.8,3.9,3.9,3.8,4.1,4.0,3.4,3.6,3.8,3.7)
results <- rbind(confint.cv(x), confint.cv(x, method = "mckay"), confint.cv(x, method = "naive"))
print(results)


gives
    method      lower      upper         CV alpha
1 modmckay 0.03617573 0.09632067 0.05263158  0.05
2    mckay 0.03618047 0.09641016 0.05263158  0.05
3    naive 0.03620185 0.09608475 0.05263158  0.05

Kind regards,
Mittyri

Complete thread:

Activity
 Admin contact
21,753 posts in 4,548 threads, 1,544 registered users;
online 5 (0 registered, 5 guests [including 3 identified bots]).
Forum time: Sunday 05:21 CEST (Europe/Vienna)

They were “so intent of making everything numerical”
that they frequently missed seeing
what was there to be seen.    Barbara McClintock

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5