Bioequivalence and Bioavailability Forum

Hi all,

I want to write a Python code for calculating 90% CIs for point estimates in a 2x2 BE setup. As a reference set I chose the NCAResult4BE dataset from Kyun-SeopBae's BE library. When I run the code:

be2x2(NCAResult4BE, c("Cmax"))

I get the following results:

$Cmax$`90% Confidence Interval of Geometric Mean Ratio (T/R)`
                 Lower Limit Point Estimate Upper Limit
90% CI for Ratio   0.9013625      0.9798396    1.065149

I can easily replicate them using nlme library:

cmax_ = lme(log(Cmax) ~ GRP + PRD + TRT, random = ~1|SUBJ,
            data = data,
            control = lmeControl(opt = "optim"))
cmax = exp(intervals(cmax_, 0.9)$fixed["TRTT",])

The problem starts in Python (3.11) and statsmodels (0.14.4).

model = smf.mixedlm("np.log(Cmax) ~ C(GRP) + C(PRD) + C(TRT, Treatment('R'))",
                    data=NCAdata,
                    groups=NCAdata["SUBJ"],
                    re_formula="~1")
result = model.fit()
conf_ints = result.conf_int(alpha=0.1)
treatment_t_index = conf_ints.index[conf_ints.index.str.contains("T")] conf_int_treatment_t = conf_ints.loc[treatment_t_index] exp_coeff = np.exp(conf_int_treatment_t)
print("Exponentiated 90% Confidence Interval for 'TreatmentT':")
print(exp_coeff)

The results are similar, but not the same:

Exponentiated 90% Confidence Interval for 'TreatmentT':
                                    0         1
C(TRT, Treatment('R'))[T.T]  0.903614  1.062496

Honestly, I have no idea what am I missing. Perhaps someone has stumbled upon a similar issue? Could the rounding be a problem or a model? Is there a way I could fix it?

D.