Bioequivalence and Bioavailability Forum

Thank you very much for your reply.
It would indeed seem that the Tsong paper was not thoroughly checked.

I also started verifying my code with DDsolver, since that is used quite often.
And for most of my tests it went fine, except for one generated data set.
All values are the same apart from the lower part of the confidence region.
In DDsolver it is: -0.361857342251203
In my code (and yours) it is: 0.361857342251203

The rest is all the same, also the conclusion would be the same, but I really dislike the idea that part of the calculation might be wrong. Because (as you know) quite a large bit is the same for the upper and lower part of the confidence region.

Here is my data set (made so you should be able to drop it in your code).
reference <- structure(c(28.7911578671327, 27.7063076528163, 28.2263135554565,
28.3774876815685, 28.2493706626467, 27.9633179734855, 28.8699375412747,
27.9697438683459, 29.1540336099923, 27.9876475987802, 28.7541180117783,
29.3043600382281, 45.5068689363729, 46.5460525217093, 46.682230512727,
47.4023771396442, 46.5409373038651, 44.3202218507459, 44.5224174823224,
48.0428499621152, 46.5199811632234, 45.1394826170731, 46.6460992057861,
47.9441445472879, 61.6499246599356, 58.8871076421716, 59.0928637200424,
57.3039074196922, 59.9450737683508, 58.6381973772281, 59.9395529973835,
60.2358174860198, 60.6281637525813, 58.6751058188512, 59.9983634443677,
57.358737462136, 66.7744608781215, 66.6843050370333, 64.5632571665535,
67.887806253632, 68.118289910638, 67.3489221222566, 68.8674535916327,
66.8528203605585, 65.9823451984023, 68.4260330454663, 66.7379172478314,
69.798117787546, 72.8623052002997, 74.4602574302653, 73.9697085359241,
72.3443125360204, 75.8127101072032, 74.4415179132332, 73.6284427938045,
73.9030111699696, 74.4950079468374, 73.6285489815188, 69.0968679635935,
73.9152589300409), .Dim = c(12L, 5L), .Dimnames = list(NULL,
    c("10", "20", "30", "40", "50")))

test <- structure(c(29.9437005373757, 27.2315750015848, 28.5315897581853,
28.9095250734653, 28.5892325261607, 27.8741008032577, 30.1406497227308,
27.8901655404087, 30.8508898945248, 27.9349248664945, 29.8511008989897,
31.2267059651144, 43.5566152658779, 46.154574229219, 46.4950192067632,
48.2953857740562, 46.1417861846084, 40.5899975518104, 41.0954866307517,
49.8965678302338, 46.0893958330041, 42.6381494676286, 46.4046909394111,
49.6498042931656, 65.0270736686448, 58.1200311242348, 58.6344213189118,
54.1620305680362, 60.7649464396827, 57.4977554618761, 60.7511445122644,
61.4918057338553, 62.4726714002591, 57.5900265659338, 60.898170629725,
54.2991056741457, 65.1779582396222, 64.9525686369017, 59.6499489607023,
67.9613216783985, 68.5375308209135, 66.61411134996, 70.4104400234003,
65.3738569457147, 63.1976690403243, 69.3068886579843, 65.086599163897,
72.7371005131835, 71.9110116134174, 75.9058921883315, 74.6795199524783,
70.616029952719, 79.287023880676, 75.8590433957511, 73.8263555971792,
74.5127765375922, 75.9927684797615, 73.8266210664652, 62.4974185216518,
74.5433959377705), .Dim = c(12L, 5L), .Dimnames = list(NULL,
    c("10", "20", "30", "40", "50")))


If you have any idea why these results would differ, you would really help me out.
I will see if I can figure out what happens in DDsolver.

OK, I now know why these values are different.
DDsolver basically calculates the lower and upper bounds as follows:

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

dm_lower <- dm - sqrt(fcrit/k)
dm_upper <- dm + sqrt(fcrit/k)


That is quite a bit simpler, is it also correct?