Bioequivalence and Bioavailability Forum

Hi wligtenberg,

❝ I think you have a typo in the definition of your mr. You have a value: 78.45, where in the article I think it is: 78.54.

Eagle eye! I never noticed it before. My copy of Tsong's article is old and is of image format so I copied the data from DDSolver (mentioned in my original post) instead to avoid manually entering the data. I double checked it after you'd pointed out the typo and it seems it is the only individual data point that's "wrong".

❝ And then I started to see if some calculations might have been off, and I checked the mean values calculated with the mean values in the Tsong article.

❝ ...

To summarise, you are right. The individual data and the mean values in the article do not agree (4 time points in reference batch and 2 time points in test batch).

❝ Similar thing with the covariance matrices, if I use your data (with typo) I get results that are more similar (excluding some rounding errors) to the Tsong article, than without the typo.

❝ Which is uncanny, so maybe I missed a memo where this is a known typo, and that you actually used the correct data.

It seems that I missed it too. I guess the only way to be sure is to ask Tsong to provide the original data. But I had a feeling that the authors of DDSolver might have done just that otherwise it's difficult to explain why the sample data generated by DDSolver give much better result with this "typo". So an alternative is to contact them.

❝ (which however does not explain why the mean values are off, in the article)

My guess is that the individual data is the correct one (maybe except the unit 4 of reference at 30 min) otherwise we'd see much bigger difference in the variance-covariance table. The article was printed rather sloppily in my opinion. There are other mistakes too, such as reference to tables are all wrong, one time point is missing for T-R mean difference, data in the last table was never referenced... So I guess we should not be shocked about the discrepancy between individual and mean values.

Actually, there some mistakes in the formula I mentioned in the original post (in the article in Encyclopedia of Biopharmaceutical Statistics by same authors) so i don't know if it's the authors or editors who were sloppy.

❝ I just wanted to verify my implementation before I am going to put it into practice.

I believe that the method is correct but that's just my opinion.

❝ I might even create a package out of this in the end, because I think it would be nice to have an R package with a couple of methods combined.

Good! Please let us know when it's available.

Just to be thorough, I paste the roughly coded function (differ very slightly from code above) that I sometimes use.

after read data from external file as data frame...

msd <- function(test, ref, dg=10){
  mt <- as.matrix(test)
  mr <- as.matrix(ref)
  nt <- nrow(mt)
  nr <- nrow(mr)
  p  <- ncol(mt)
  dg <- dg

  mean_t    <- cbind(colMeans(mt))
  mean_r    <- cbind(colMeans(mr))
  mean_diff <- mean_t - mean_r

  st <- cov(mt)
  sr <- cov(mr)
  sp <- ((nt - 1) * st + (nr - 1) * sr) / (nt  + nr -2)
  sp_inv  <- solve(sp)
  m_dg    <- cbind(rep(dg, p))
  dm_max  <- sqrt(t(m_dg) %*% sp_inv %*% m_dg)

  dm <- sqrt(t(mean_diff) %*% sp_inv %*% mean_diff)

  df    <- nt + nr - p - 1
  k     <- nt * nr / (nt + nr) * df / ((nt + nr - 2) * p)
  fcrit <- qf(0.9, p, df)
  h2    <- nt * nr / (nt + nr) * dm^2

  y1 <- mean_diff %*% (1 + sqrt(fcrit / (k * (t(mean_diff) %*% sp_inv %*% mean_diff))))
  y2 <- mean_diff %*% (1 - sqrt(fcrit / (k * (t(mean_diff) %*% sp_inv %*% mean_diff))))

  ci1 <- sqrt(t(y1) %*% sp_inv %*% y1)
  ci2 <- sqrt(t(y2) %*% sp_inv %*% y2)

  dm_lower <- min(ci1, ci2)
  dm_upper <- max(ci1, ci2)

  result <- data.frame(cbind(dm, dm_lower, dm_upper, dg, dm_max, h2))
  names(result) <- c("Dm", "Dm_lower", "Dm_upper", "Dg (%)", "Dm_max", "Hotelling T Square")
  return(result)
}