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

posted by mittyri – Russia, 2018-03-21 10:35 (875 d 22:26 ago) – Posting: # 18576
Views: 1,739

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,009 posts in 4,379 threads, 1,460 registered users;
online 17 (0 registered, 17 guests [including 12 identified bots]).
Forum time: Thursday 09:01 UTC (Europe/Vienna)

If you obey all the rules,
you will miss all the fun.    Katharine Hepburn

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