Ensayos de Bioequivalencia de Promedios (Borrador)

Modelos de efectos fijos

Introducción

En la guía Guideline on the Investigation of Bioequivalence (EMA 2010) se describe el diseño experimental y, de forma literal, el método de análisis estadístico recomendados para la realización de los ensayos de bioequivalencia: transformación logarítmica de la variable utilizada (ABC, C_MAX); modelo de análisis de la varianza que incluye los factores tratamiento, período, secuencia e individuos dentro de secuencia, todos como efectos fijos; la prueba de bioequivalencia es el intervalo de confianza 90% de la diferencia entre los promedios del medicamento ensayado (T) y el medicamento de referencia (R).

El procedimiento habitual para realizar el análisis de la varianza se basa en la división de la varianza total entre los orígenes de variación del modelo estadístico. El resultado final es una tabla con las fórmulas para calcular las sumas cuadráticas, grados de libertad y las medias cuadráticas para cada origen de variación. Para interpretar esta tabla es necesario además conocer los valores esperados de las medias cuadráticas para cada origen de variación. Para una descripción detallada de este procedimiento aplicado al ensayo cruzado de dos períodos y dos secuencias ver Chow y Liu (2009, página 66). Este procedimiento está implementado en la función test2x2() de la libraría BE.

En el documento Annex I: statistical analysis methods compatible with CHMP guideline (EMA, 2016), y siguiendo la tendencia iniciada por el FDA, el método para el análisis de los datos se especifica mediante un procedimiento SAS. Se describen tres métodos:

• Método A: modelo compatible con las directrices CHMP. Modelo de efectos fijos idéntico al descrito en la guía, salvo que el código SAS es igualmente aplicable a los ensayos cruzados de dos o más secuencias y dos o más períodos. Se utiliza el procedimiento glm que implementa la estimación por el método de mínimos cuadrados ordinario de los parámetros del modelo general lineal.
• Método B: ligera modificación del modelo compatible con las directrices CHMP. Modelo de efectos mixtos que incluye dos términos aleatorios para explicar las variabilidades intra-individual e inter-individual. Se utiliza el procedimiento mixed que implementa la estimación por el método máxima verisimilitud de los parámetros de modelos lineales de efectos mixtos.
• Método C (sín título) Es el método incluido en la guía Statistical approaches to establishing bioequivalence (FDA, 2001). Recientemente se ha publicado Bioequivalence Studies with Pharmacokinetic Endpoints for Drugs Submitted Under an ANDA. Guidance for Industry (FDA 2021).

  Al igual que el método B, está basado en la estimación de los parámetros de un modelo lineal de efectos mixtos por el método de máxima verisimilitud.

Los métodos A y B están implementados en la librería replicateBE (Schütz, Tomashevskiy y Labes(2021)).

Para más información sobre las guias para los ensayos de bioequivalencia publicados por la EMA ver CE

Objetivos

El objetivo de este documento es explicar la implementación del método A en R. Los puntos a desarrollar son:

• análisis del código SAS del método A,
• descripción del modelo estadístico y de la prueba de bioequivalencia de promedios implícitos en el método A,
• descripción de la función aov() y su aplicación a la estimación de la biodisponibilidad relativa y de su intervalo de confianza,
• descripción de los modelos I y III del análisis de la varianza,
• validación de la implementación del método A en R utilizando los datos publicados por la EMA con tal fin.

Ejemplos de trabajo

Ver código.

rm(list=ls(all=TRUE))

## data from Ekbohm and Melander, Biometrics 45 (1989), 1249

subj <- rep(1:8, 4)
per <- rep(1:4, each = 8)
seq <- rep(c(rep("TRTR", 4), rep("RTRT",4)),4)
trt <- rep(c(rep(c("T", "R"), each = 4), rep(c("R", "T"), each = 4)), 2)
auc <- c(164.9, 268.9, 95.0, 62.3, 184.3, 173.5, 194.7, 151.3,
    137.3, 166.7, 107.8, 190.2, 147.2, 175.7, 161.0, 168.8,
    162.4, 234.9, 108.0, 127.7, 91.2, 146.5, 153.7, 189.9,
    110.9, 209.4, 173.4, 189.8, 109.2, 183.7, 207.1, 241.9)

furo.dat <- data.frame(subj = factor(subj), per = factor(per), seq = factor(seq),
    trt = trt, auc = auc)
rm("subj", "per", "seq", "trt")

subj <- c(1,5,6,8,11,15,16,18,19,21,22,
    2,3,4,7,9,10,12,13,14,17,20,23,24)
subj <- rep(subj, 2)
seq <- rep(c(rep('RT', 11), rep('TR',13)), 2)
trt <- c(rep('R', 11), rep('T', 13), rep('T', 11), rep('R', 13))
per <- c(rep(1, 24), rep(2, 24))
auc <- c(110.0, 63.5, 34.6, 163.0, 53.9, 20.3, 26.0, 17.8, 48.1, 23.8, 68.7,
    64.5, 84.2, 63.5, 20.7, 18.7, 13.8, 14.2, 22.7, 91.8, 120.0, 95.8, 31.1, 14.7,
    126.0, 55.3, 49.9, 83.7, 117.0, 25.6, 10.2, 23.4, 27.2, 15.4, 112.0,
    47.4, 79.8, 61.9, 13.8, 24.9, 4.57, 22.2, 8.7, 45.0, 76.3, 190.0, 19.8, 26.69)

gran.dat <- data.frame(seq = factor(seq), subj = factor(subj), per = factor(per),
    trt = factor(trt), auc = auc)
rm('subj', 'seq', 'trt', 'per', 'auc')

eryth.dat <- data.frame(
  subj = factor(rep(1:18, 2)),
  seq = as.factor(rep(c(rep("TR", 9), rep("RT", 9)), 2)),
  per = factor(rep(c(1,2), each = 18)),
  trt = as.factor(c(rep("T", 9), rep("R", 9), rep("R", 9), rep("T", 9))),
  auc = c(2.52, 8.87, 0.79, 1.68, 6.95, 1.05, 0.99, 5.60, 3.16,
          4.98, 7.14, 1.81, 7.34, 4.25, 6.66, 4.76, 7.16, 5.52,
          5.47, 4.84, 2.25, 1.82, 7.87, 3.25, 12.39, 4.77, 1.88,
          3.19, 9.83, 2.91, 4.58, 7.05, 3.41, 2.49, 6.18, 2.85
          )
)

 

Código SAS método A

El código SAS incluido en la descripción del método A del anexo I es el siguiente:

proc glm data=replicate;
class formulation subject period sequence;
model logDATA= sequence subject (sequence) period formulation;
estimate "test-ref" formulation -1+1;
test h=sequence e=subject(sequence);
lsmeans formulation / adjust=t pdiff=control("R") CL alpha=0.10;
run;

Descripción del código:

• proc glm: aplica el modelo general lineal; su equivalente en R es la función lm().
• data=replicate: importa la tabla de datos replicate.
• class ...: declara factores (variables categóricas) las variables formulation, subject, period y sequence. En R, la especificación del tipo de variable se realiza en la tabla de datos (objeto data.frame).
• model ...: especifica el modelo log-lineal (ver ecuación 1). El término subject(sequence) indica que el término subject está anidado dentro de sequence.
• estimate "test-ref" formulation -1+1: contraste para estimar la diferencia entre los promedios de los dos tratamientos.
• test h=sequence e=subject(sequence): determina la varianza inter-individual como término error para probar la hipótesis nula del efecto secuencia.

Esta sentencia es contradictoria con la descripción literal del método A, en la que se declara que todos los términos del modelo son efectos fijos y el modelo contiene sólo un término error, \(\sigma_W^2\). Ver la sección Modelo Estadístico más abajo.

• lsmeans ...: estimación de las medias mínimocuadráticas de ambos tratamientos y los respectivos intervalos de confianza 90% utilizando la distribución t-student.

Bioequivalencia de promedios

Modelo estadístico

El modelo estadístico utilizado descrito en EMA (2016) para evaluar la bioequivalencia en los ensayos cruzados e implícito en el método A del anexo I es el denominado log-lineal de efectos fijos,

\[\begin{equation} y_{ijk}~=~ln(x_{ijk})~=~\mu + G_k + S_{ik} + P_j + F_{(jk)} + \epsilon_{ijk} \tag{1} \end{equation}\]

donde:

• y_ijk es la observación realizada sobre el individuo i (i = 1 … n_k) de la secuencia k (k = 1, … K) en el período j (j = 1, … J).
• μ es la media general
• G_k es el efecto de la secuencia k (k = 1 … K).
• S_ik es el efecto del individuo i perteneciente a la secuencia k (i = 1 … n_k).
• P_j, el efecto del período j (j = 1 … J)
• F_(j,k) es el efecto de la formulación que recibe el individuo de la secuencia k durante el período j.
• ε_ijk es el término aleatorio del modelo que se asume,

\[\begin{equation} \epsilon \sim N(0~,~\mathbf{I}_{2~n}~\sigma_W^2) \tag{2} \end{equation}\]

\(\sigma_W^2\) (within-subject variance) recoje todas las fuentes de variación no explicadas. El principal componente es la variabilidad intra-individual, pero también hemos de prestar atención al error del método analítico y al error del método para calcular la variable analizada (por ejemplo, el error al calcular ABC utilizando la regla trapezoidal).

Los factores secuencia, individuo dentro de secuencia y período se consideran factores perturbadores, es decir, que no intervienen en la hipótesis de bioequivalencia y que se introducen en el modelo estadístico para reducir la varianza residual y aumentar la potencia de la prueba de bioequivalencia.

No pocos autores se decantan por interpretar el término S_ik como un término aleatorio con distribución independiente de ε y

\[\begin{equation} S_{ik} \sim N(0~,~\mathbf{I}_n~\sigma^2_W) \end{equation}\]

No obstante hemos de tener en cuenta:

• el resultado de la prueba de bioequivalencia de promedios es idéntica en ambos casos,
• la estimada de \(\sigma_B^2\) obtenida de la tabla del ANOVA no es insesgada si el diseño experimental no es balanceado (igual número de individuos por secuencia),
• la estimada de \(\sigma_B^2\) obtenida de la tabla del ANOVA puede ser incluso negativa si \(\sigma^2_B < \sigma^2_W\),

Si se desea estimar \(\sigma_W^2\) y \(\sigma_B^2\) es recomendable utilizar el método B citado más arriba.

Criterio de bioequivalencia

El criterio de bioequivalencia de la formulación T con la formulación R es:

\[\begin{equation} \theta_L \le \frac{\tilde{\mu}_T}{\tilde{\mu}_R} \le \theta_U \tag{1} \end{equation}\]

donde \(\tilde{\mu}_T\) y \(\tilde{\mu}_R\) son las medias geométricas del parámetro farmacocinético utilizado para evaluar la bioequivalencia de las formulaciones T y R respectivamente, y \([\theta_L~,~\theta_U]\) el intervalo de bioequivalencia. Dado que el modelo estadístico se basa en los logaritmos de las observaciones, el criterio de bioequivalencia también puede expresarse en la forma

\[\begin{equation} \delta_L \le \mu_T - \mu_R \le \delta_U \tag{2} \end{equation}\]

donde \(\mu_T\) y \(\mu_R\) son las medias de los logaritmos del parámetro utilizado para evaluar la bioequivalencia.

En la guía EMA se especifican los siguientes intervalos de bioequivalencia:

• Caso general:\([\theta_L~,~\theta_U]~=~[0.80~,~1.25]\), o bien, en escala logarítmica \([\delta_L~,~\delta_U]~=~[-0.2231~,~0.2231]\).
• Fármacos extrecho intervalo terapéutico:
• Fármacos con elevada variabilidad intra-individual:

La razón por la que el criterio de bioequivalencia se basa en las medias geométricas y no en la las medias aritméticas es la siguiente: sea la variable aleatoria Y ~ N(μ , σ²)), y X = e^Y; entonces la esperanza (valor medio) y la varianza de X son,

\[\begin{align} E(x)~=~e^{\left(\mu + \frac{\sigma^2}{2} \right)} \\ VAR(x)~=~\left(e^{\sigma^2} - 1 \right)~e^{2~\mu + \sigma^2} \tag{5} \end{align}\]

Ver que el valor esperado de X tiene un sesgo igua a σ²/2. No ocurre lo mismo con la mediana,

\[\begin{equation} MEDIANA(X)~=~e^{\mu} \tag{6} \end{equation}\]

Prueba de bioequivalencia de promedios

Sea \(\hat{\delta}~=~\hat{F}_T - \hat{F}_R\) la diferencia entre las estimadas mínimocuadráticas de los promedios del logaritmo del parámetro farmacocinético analizado de las formulaciones T y R. La formulación T se declara bioequivalente con la formulación R si

\[\begin{equation} prob \left[\hat{\delta} \in [\delta_L~,~\delta_U] \right] \ge 1 - \alpha \end{equation}\]

para \(\alpha = 0.10\). Es decir, si el intervalo de confianza 90% de la diferencia entre formulaciones se encuentra dentro del intervalo de bioequivalencia.

No hemos de olvidar que esta prueba, propuesta inicialmente por Metzler (xxxx), carece de fundamento teórico, y que su aplicación está justificada porque es operacionalmente equivalente con la prueba TOST propuesta por Shuirmann (xxxx).

   

Implementación del método A

El procedmiento que se describe en este documento está basado en la función aov() que realiza el análisis de varianza basado en el modelo general lineal. La figura 1 muestra el diagrama de flujo de esta función así como su aplicación a la estimación del intervalo de confianza de la biodisponibilidad relativa.

Obsérvaciones:

• La función aov() llama internamente a las funciones lm() y anova().
• La función lm() identifica si las variables independientes son variables categóricas (factor), en cuyo caso calcula la matriz de diseño y estima los parámetros del modelo de regresión.
• La etapa de estimación está representada en la figura por la ecuación para estimar los parámetros por el método de mínimos cuadráticos ordinarios (b = (X^T X)^-1) X^T y). Sin embargo en los modelos de ANOVA es frecuente que la matriz X^TX sea singular. La función lm() calcula las estimadas del modelo por el método qr; en la librería sasLM se calcula la inversa generalizada utilizando el operador SWEEP.

 

Figura 1. Diagrama de flujo de la función aov() y su aplicación a la evaluación de los ensayos de bioequivalencia

 

Función aov()

La función aov() realiza el análisis de la varianza en dos etapas:

1. Estimación de los parámetros del modelo general lineal utilizando la función lm().

   La función lm() identifica el tipo de variables en la tabla de datos, numéricas (numeric) o categóricas (factor ); en el segundo caso calcula la matriz de diseño.

   .
2. Cálculo de la tabla del modelo 1 del ANOVA utilizando la función anova()

   Para imprimir la tabla del análisis de la varianza utilizar la función genérica summary()

   .

La sintaxis de la función aov() es:

T <- aov(formula , data, ...)

que es equivalente a

T <- anova(lm(formula , data, ...))

Los argumentos de la función son:

• formula describe el modelo estadístico utilizado; los nombres de las variables deben coincidir con los de los campos de la tabla de datos. Para indicar la presencia de factores anidados se utiliza el operador %in%. La fórmula para el modelo estadístico utilizado en este documento (ecuación 1) es:

• data es la tabla de datos (objeto data.frame); debe contener los campos incluidos en el modelo estadístico.

• ... indica otros argumentos.

 

La función aov() devuelve un objeto tipo list que contiene los siguientes objetos generados internamente por la función lm():

 [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values"
 [6] "assign"        "qr"            "df.residual"   "contrasts"     "xlevels"      
[11] "call"          "terms"         "model"

Ejemplo de trabajo

Datos furosemida; data.frame:

str(furo.dat)

## 'data.frame':    32 obs. of  5 variables:
##  $ subj: Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 2 ...
##  $ per : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 2 2 ...
##  $ seq : Factor w/ 2 levels "RTRT","TRTR": 2 2 2 2 1 1 1 1 2 2 ...
##  $ trt : Factor w/ 2 levels "R","T": 2 2 2 2 1 1 1 1 1 1 ...
##  $ auc : num  164.9 268.9 95 62.3 184.3 ...

Medias período x secuencia

N1 <- with(furo.dat, tapply(auc, list(seq, per), length))
M1 <- round(with(furo.dat, tapply(log(auc), list(seq, per), mean)), 4)
T1 <- cbind(N1[,1], M1)
colnames(T1) <- c("n[k]", "per1", "per2","per3", "per.4")
T1

##      n[k]   per1   per2   per3  per.4
## RTRT    4 5.1659 5.0927 4.9454 5.1821
## TRTR    4 4.8464 4.9917 5.0203 5.1136

Modelo estadístico:

log(auc) ~ per + seq + subj %in% seq + trt

Matriz de diseño

Tal como se detalla más abajo, la función lm() genera internamente la matriz de diseño aplicando la función model.matrix() utilizando por defecto el contraste por tratamiento (contr.treatment), con el objetivo de interpretar el modelo estadístico del ANOVA (ecuación X) mediante un modelo de regresión. Para el ejemplo de trabajo, el modelo de regresión es:

\[\begin{align} y_i~=~\beta_0 + \beta_{trtT}~x_{1,i} + \beta_{per2}~x_{2,i} + \beta_{per3}~x_{3,i} + \beta_{per4}~x_{4,i} + \beta_{seqTRTR}~x_{5,i} \\ + \sum secuencia : individuos \end{align}\]

Para simplificar la exposición hemos agrupado en el término \(\sum secuencia:individuos\) los K x (n-1) términos correspondientes a individuo dentro de secuencia.

Para construir la matriz de diseño:

• Cada factor del modelo estadístico del ANOVA se interpreta en el modelo de regresión a través de un número de variables igual al número de niveles del factor menos uno.

  Por ejemplo, el factor trt tiene dos niveles y por lo tanto es necesario una variable. El factor per tiene 4 niveles y son necesarias 3 variables.

• Para cada factor, el nivel de referencia es el primero de los niveles ordenados alfabéticamente.

  El nivel de referencia del factor tratamiento es R, y en modelo de regresión incluye un parámetro que estima la diferencia entre el segundo nivel, T y el nivel de referencia \(\beta_{trtT}\). El factor período tiene 4 niveles y el nivel de referencia es el período etiquetado con 1; el modelo de regresión incluye entonces tres variables para estimar las diferencias entre los períodos 2, 3 y 4, y el período 1.

• La variable toma el valor 0 para el nivel de referencia, y el valor 1 para los restantes niveles.

La tabla siguiente muestra los valores de las variables del modelo de regresión para los factores no anidados para el diseño cruzado de dos secuencias y 4 períodos. La variables son:

• \(x_1\): factor tratamiento, nivel T; toma el valor 0 para el nivel de referencia (R), y 1 para el segundo nivel.
• \(x_2\): factor período, nivel 2; toma el valor 0 para el nivel de referencia (período 1), y 1 para el nivel 2.
• \(x_3\): factor período, nivel 3; toma el valor 0 para el nivel de referencia, y 1 para el nivel 3.
• \(x_4\): factor período, nivel 4; toma el valor 0 para el nivel de referencia, y 1 para el nivel 4.
• \(x_5\): factor secuencia, nivel TRTR; toma el valor 0 para la secuencia de referencia (RTRT), y 1 para el nivel 2 (secuencia TRTR).

Variables del modelo de regresión para el ensayo cruzado de 2 secuencias y 4 períodos

seq	per	trt	x1	x2	x3	x4	x5
RTRT	1	R	0	0	0	0	0
RTRT	2	T	1	1	0	0	0
RTRT	3	R	0	0	1	0	0
RTRT	4	T	1	0	0	1	0
TRTR	1	T	1	0	0	0	1
TRTR	2	R	0	1	0	0	1
TRTR	3	T	1	0	1	0	1
TRTR	4	R	0	0	0	1	1

Estimada de la biodisponibilidad relativa

El coeficiente \(\beta_{trtT}\) es la diferencia entre el efecto del tratamiento ensayado y el tratamiento de referencia, es decir,

\[\begin{equation} E[\beta_{trtT} \mid datos]~=~F_T - F_R = \delta \end{equation}\]

El vector coefficients tiene los elementos etiquetados (named vector), por lo que es inmediado identificar la estimada de \(\beta_{trtT}\). Los nombres de los elementos del vector de estimadas para el ejemplo de trabajo son:

furo.aov <- aov(log(auc) ~ trt + per + seq + subj %in% seq ,
                  furo.dat)
names(furo.aov$coefficients)

##  [1] "(Intercept)"   "trtT"          "per2"          "per3"         
##  [5] "per4"          "seqTRTR"       "seqRTRT:subj2" "seqTRTR:subj2"
##  [9] "seqRTRT:subj3" "seqTRTR:subj3" "seqRTRT:subj4" "seqTRTR:subj4"
## [13] "seqRTRT:subj5" "seqTRTR:subj5" "seqRTRT:subj6" "seqTRTR:subj6"
## [17] "seqRTRT:subj7" "seqTRTR:subj7" "seqRTRT:subj8" "seqTRTR:subj8"

Obsérvese que las etiquetas se forman por la unión de la etiqueta del factor y la etiqueta del nivel. Las estimadas de los coeficientes de los términos no anidados son:

furo.aov$coefficients[1:6]

## (Intercept)        trtT        per2        per3        per4     seqTRTR 
##  5.19153240 -0.01879736  0.03605287 -0.02329578  0.14170928 -0.26420059

Es decir, \(\hat{\delta}\) = -0.0188 y la estimada de la biodisponibilidad relativa \(\hat{\theta} = e^{\hat{\beta}_{trtT}}\) = 0.9814.

Intervalo de confianza de la biodisponibilidad relativa

El intervalo de confianza de trtT se estima utilizando la función confint(). Su sintásis es:

confint(objet , parm , level , ...)

Los argumentos de la función son:

• object, es el objeto devuelto por la función utilizada para estimar los parámetros,
• parm, la etiqueta del parámetro de interés,
• level, nivel probabilidad.

Los intervalos de confianza 90% de \(\hat{\delta}\) y de \(\hat{\theta}\) para el ejemplo de trabajo son:

furo.abe <- c(furo.aov$coefficients["trtT"],
             confint(furo.aov, "trtT", level = 0.90))
furo.abe <- rbind(round(furo.abe, 5), round(100 * exp(furo.abe), 2))
row.names(furo.abe) <- c("delta" , "theta")
colnames(furo.abe) <- c("estim", "L inf", "L sup")
furo.abe

##         estim    L inf     L sup
## delta -0.0188 -0.19351   0.15592
## theta 98.1400 82.41000 116.87000

El intervalo de confianza 90% de \(\hat{\theta}\) [82.41 - 116.87] se encuentra dentro del intervalo de bioequivalencia de promedios, [80.00 , 125.00] y por lo tanto se acepta la bioequivalencia del producto ensayado.

Resumen del procedimiento

El cuadro siguiente agrupa las sentencias R para evaluar el ensayo de bioequivalencia.

datos <- data.frame( ......)                             # construcción de la tabla de datos
datos.aov <- aov(modelo_estadístico , datos)             # análisis de la varianza
summary(datos.aov)                                       # tabla del análisis de la varianza
resultado <- c(gran.aov$coefficients["trtT"],            # estimada mínimocuadrática de la
             confint(gran.aov, "trtT", level = 0.90))    # diferencia de promedios y su
resultado <- rbind(resultado, exp(resultado))            # intervalo de confianza, y
                                                         # biodisponibilidad relativa y su
                                                         # intervalo de confianza

Validación

Para validar la aplicación de otras aplicaciones informáticas, el documento Statistical analysis methods compatible with CHMP guideline incluye los resultados obtenidos con el procedimiento SAS de dos series de datos recogidos en los anexos II (Annex II) y III (Annex III). La librería replicateBE incluye estas dos tablas junto con otras utilizadas con el mismo fin (Schütz y cols. 2020).

La tabla siguiente muestra los resultados recogidos en el anexo I para el método A.

Resultados datos validación; biodisponibilidad relativa estimada (%) e intervalo de confianza 90%

tabla	diseño	n[k]	estimada	lim. inf.	lim.sup
Set I (annex II)	RTRT/TRTR	38/39	115.16	107.11	124.89
Set II (annex III)	RRT/RTR/TRR	8/8/8	102.26	97.32	107.46

Datos Annex II

Estructura de la tabla de datos:

str(rds01)

## 'data.frame':    298 obs. of  6 variables:
##  $ subject  : Factor w/ 77 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ period   : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ sequence : Factor w/ 2 levels "RTRT","TRTR": 1 1 1 1 2 2 2 2 2 2 ...
##  $ treatment: Factor w/ 2 levels "R","T": 1 2 1 2 2 1 2 1 2 1 ...
##  $ PK       : num  2286 1956 1346 2856 3152 ...
##  $ logPK    : num  7.73 7.58 7.2 7.96 8.06 ...
##  - attr(*, "rset")= chr "rds01"

Número de individuos por secuencia y promedios de secuencias x períodos:

Nset1 <- with(rds01, tapply(PK, list(sequence, period), length))
Mset1 <- round(with(rds01, tapply(PK, list(sequence, period), mean)), 3)
Tset1 <- cbind(Nset1[,1], Mset1)
colnames(Tset1) <- c("n[k]", "per1", "per2","per3", "per.4")
Tset1

##      n[k]     per1     per2     per3    per.4
## RTRT   38 3193.080 3562.105 3417.279 3772.025
## TRTR   39 3923.641 3663.481 4017.557 3765.765

ANOVA:

rds01.aov <- aov(log(PK) ~ sequence + subject %in% sequence + period + treatment,
                  rds01)
summary(rds01.aov)

##                   Df Sum Sq Mean Sq F value  Pr(>F)    
## sequence           1   0.01  0.0077   0.048 0.82710    
## period             3   0.70  0.2328   1.455 0.22783    
## treatment          1   1.77  1.7681  11.051 0.00104 ** 
## sequence:subject  75 214.13  2.8551  17.845 < 2e-16 ***
## Residuals        217  34.72  0.1600                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de bioequivalencia:

rds01.ci <- c(rds01.aov$coefficients["treatmentT"],
             confint(rds01.aov, "treatmentT", level = 0.90))
rds01.ci <- rbind(round(rds01.ci, 5), round(100 * exp(rds01.ci), 2))
colnames(rds01.ci) <- c("estim", "Linf", "Lsup")
row.names(rds01.ci) <- c("delta" , "theta")
rds01.ci

##           estim      Linf     Lsup
## delta   0.14547   0.06865   0.2223
## theta 115.66000 107.11000 124.8900

Datos Annex III

Estructura de la tabla de datos:

str(rds02)

## 'data.frame':    72 obs. of  6 variables:
##  $ subject  : Factor w/ 24 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ period   : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
##  $ sequence : Factor w/ 3 levels "RRT","RTR","TRR": 2 2 2 2 2 2 1 1 1 3 ...
##  $ treatment: Factor w/ 2 levels "R","T": 1 2 1 1 2 1 1 1 2 2 ...
##  $ PK       : num  4054 3970 3749 2986 2379 ...
##  $ logPK    : num  8.31 8.29 8.23 8 7.77 ...
##  - attr(*, "rset")= chr "rds02"

Número de individuos por secuencia y promedios de secuencias x períodos:

Nset2 <- with(rds02, tapply(PK, list(sequence, period), length))
Mset2 <- round(with(rds02, tapply(PK, list(sequence, period), mean)), 3)
Tset2 <- cbind(Nset2[,1], Mset2)
colnames(Tset2) <- c("n[k]", "per1", "per2","per3")
Tset2

##     n[k]     per1     per2     per3
## RRT    8 2694.350 2822.800 3111.500
## RTR    8 2974.550 2824.250 2980.700
## TRR    8 3095.775 3043.675 2987.675

ANOVA:

rds02.aov <- aov(log(PK) ~ sequence + subject %in% sequence + period + treatment,
                  rds02)
summary(rds02.aov)

##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## sequence          2 0.0239 0.01197   0.857    0.431    
## period            2 0.0396 0.01982   1.420    0.252    
## treatment         1 0.0080 0.00802   0.575    0.452    
## sequence:subject 21 2.9490 0.14043  10.061 7.58e-11 ***
## Residuals        45 0.6281 0.01396                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prueba de bioequivalencia:

rds02.ci <- c(rds02.aov$coefficients["treatmentT"],
             confint(rds02.aov, "treatmentT", level = 0.90))
rds02.ci <- rbind(round(rds02.ci, 5)
                  , round(100 * exp(rds02.ci), 2))
colnames(rds02.ci) <- c("estim", "Linf", "Lsup")
row.names(rds02.ci) <- c("delta" , "theta")
rds02.ci

##           estim     Linf      Lsup
## delta   0.02239 -0.02721   0.07199
## theta 102.26000 97.32000 107.46000

Ejemplos y aplicaciones

Eritromicina

Ne <- with(eryth.dat, tapply(auc, list(seq, per), length))
Me <- round(with(eryth.dat, tapply(auc, list(seq, per), mean)), 3)
cbind(Ne[, 1], Me)

##          1     2
## RT 9 5.513 4.721
## TR 9 3.512 4.949

Prueba de bioequivalencia:

s2 <- aov(log(auc) ~ seq + subj %in% seq + per + trt,
                  eryth.dat)
r2 <- c(s2$coefficients["trtT"],            # estimada mínimocuadrática de la
             confint(s2, "trtT", level = 0.90))    # diferencia de promedios y su
r2 <- rbind(r2, 100 * exp(r2))
r2

##          trtT                       
## r2 -0.3410759 -0.6645745 -0.01757729
##    71.1004943 51.4492394 98.25762929

ANOVA

summary(s2)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## seq          1  1.305  1.3053   4.224 0.0566 .
## per          1  0.196  0.1959   0.634 0.4376  
## trt          1  1.047  1.0470   3.388 0.0843 .
## seq:subj    16  8.343  0.5215   1.688 0.1528  
## Residuals   16  4.944  0.3090                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Apéndice

Datos ejemplo de trabajo 1

Clayton y Leslie (1981)

eryth.tab <- reshape(eryth.dat, v.names = c("trt" , "auc"),
                    timevar = "per", idvar = "subj", direction = "wide")
eryth.tab <- eryth.tab[order(eryth.tab$seq, decreasing = TRUE),]
colnames(eryth.tab) <- c("individuo", "secuencia", "formulación", "período 1","formulación",   "período 2")
eryth.tab[, c(1, 2, 4, 6)]

##    individuo secuencia período 1 período 2
## 1          1        TR      2.52      5.47
## 2          2        TR      8.87      4.84
## 3          3        TR      0.79      2.25
## 4          4        TR      1.68      1.82
## 5          5        TR      6.95      7.87
## 6          6        TR      1.05      3.25
## 7          7        TR      0.99     12.39
## 8          8        TR      5.60      4.77
## 9          9        TR      3.16      1.88
## 10        10        RT      4.98      3.19
## 11        11        RT      7.14      9.83
## 12        12        RT      1.81      2.91
## 13        13        RT      7.34      4.58
## 14        14        RT      4.25      7.05
## 15        15        RT      6.66      3.41
## 16        16        RT      4.76      2.49
## 17        17        RT      7.16      6.18
## 18        18        RT      5.52      2.85

Ekbohm y Melander

furo.tab <- reshape(furo.dat, v.names = c("trt" , "auc"),
                    timevar = "per", idvar = "subj", direction = "wide")
furo.tab <- furo.tab[order(furo.tab$subj),]
# colnames(eryth.tab) <- c("individuo", "secuencia", "formulación", "período 1","formulación",   "período 2")
# eryth.tab[, c(1, 2, 4, 6)]
furo.tab

##   subj  seq trt.1 auc.1 trt.2 auc.2 trt.3 auc.3 trt.4 auc.4
## 1    1 TRTR     T 164.9     R 137.3     T 162.4     R 110.9
## 2    2 TRTR     T 268.9     R 166.7     T 234.9     R 209.4
## 3    3 TRTR     T  95.0     R 107.8     T 108.0     R 173.4
## 4    4 TRTR     T  62.3     R 190.2     T 127.7     R 189.8
## 5    5 RTRT     R 184.3     T 147.2     R  91.2     T 109.2
## 6    6 RTRT     R 173.5     T 175.7     R 146.5     T 183.7
## 7    7 RTRT     R 194.7     T 161.0     R 153.7     T 207.1
## 8    8 RTRT     R 151.3     T 168.8     R 189.9     T 241.9

 

Bibliografía

Chow, S.-C. y Liu, J-p (2009) Design and analysis of bioavailability and bioequivalence studies, 3th ed. CRC Press, Boca Raton

Schütz, H., Labes, D., Tomashevskiy, M., González-de la Parra, M., Shitova, A., Fuglsang, A. (2020). Reference Datasets for Studies in a Replicate Design Intended for Average Bioequivalence with Expanding Limits The APPS J 22:44