Estudio de bioequivalencia

Los estudios de bioequivalencia consisten en experimentos que perimten establecer si son o no similares las reacciones a formas o presentaciones diferentes en cuanto a la concentración del compuesto activo en la sangre durante un periodo de tiempo.

Se suele utilizar el Área Bajo la Curva (AUC, por sus siglas en inglés) para evaluar el efecto de drogas sobre sujetos humanos o animales.

Después de suministrada la dosis se toma una muestra de sangre cada media hora durante cuatro horas y se analiza la concentración del elemento activo en la sangre. Seguidamente se calcula el Área Bajo la Curva (AUC) de ese resultado obtenido.

Debido a que se sabe que hay diferencia entre sujetos y además también hay diferencias en el tiempo se tomaron encuenta estas dos situaciones escogiendo a tres individuos y en tres momentos diferentes.

Entre cada periodo de prueba se realiza un proceso de “borrado” o de lavado del elemento activo en el sujeto.

Se quiere saber si hay diferencia entre tres presentaciones, A = Solución, B = Tableta, C = Cápsula de una determinada droga.

Se propone un diseño 3x3 de cudrado latino.

Programación del experimento

library(agricolae)

presentacion <- c("A=Solución", "B=Tableta", "C=Cápsula")
diseñoCuadLat <- design.lsd(presentacion , seed = 4)
programacion <- diseñoCuadLat$book
levels(programacion$row) <- c("Sujeto 1", "Sujeto 2", "Sujeto 3")
levels(programacion$col) <- c("Semana 1", "Semana 2", "Semana 3")

programacion
##   plots      row      col presentacion
## 1   101 Sujeto 1 Semana 1   A=Solución
## 2   102 Sujeto 1 Semana 2    B=Tableta
## 3   103 Sujeto 1 Semana 3    C=Cápsula
## 4   201 Sujeto 2 Semana 1    B=Tableta
## 5   202 Sujeto 2 Semana 2    C=Cápsula
## 6   203 Sujeto 2 Semana 3   A=Solución
## 7   301 Sujeto 3 Semana 1    C=Cápsula
## 8   302 Sujeto 3 Semana 2   A=Solución
## 9   303 Sujeto 3 Semana 3    B=Tableta
write.csv2(programacion, file = "programacion.csv",
           row.names = FALSE,
           fileEncoding = "latin1")

Resultados de un estudio de bioequvalencia

Tratamiento y resultados de Áreas Bajo la Curva (AUC, siglas en inglés) para un estudio de bioequivalentealencia.

library(readxl)
bioequivalente <- read_excel("resultados.xlsx")
bioequivalente
## # A tibble: 9 x 5
##   plots Sujeto   Periodo  Presentacion   AUC
##   <dbl> <chr>    <chr>    <chr>        <dbl>
## 1  101. Sujeto 1 Semana 1 A=Solución   1186.
## 2  102. Sujeto 1 Semana 2 B=Tableta     642.
## 3  103. Sujeto 1 Semana 3 C=Cápsula    1183.
## 4  201. Sujeto 2 Semana 1 B=Tableta     984.
## 5  202. Sujeto 2 Semana 2 C=Cápsula    1135.
## 6  203. Sujeto 2 Semana 3 A=Solución   1305.
## 7  301. Sujeto 3 Semana 1 C=Cápsula    1426.
## 8  302. Sujeto 3 Semana 2 A=Solución   1540.
## 9  303. Sujeto 3 Semana 3 B=Tableta     873.

Análisis exploratorio

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
bioequivalente  %>% 
  group_by(Presentacion) %>%
  summarise(promedio = mean(AUC),
            desvEst = sd(AUC))
## # A tibble: 3 x 3
##   Presentacion promedio desvEst
##   <chr>           <dbl>   <dbl>
## 1 A=Solución      1344.    180.
## 2 B=Tableta        833.    174.
## 3 C=Cápsula       1248.    156.
library(dplyr)
bioequivalente  %>% 
  group_by(Sujeto) %>%
  summarise(promedio = mean(AUC),
            desvEst = sd(AUC))
## # A tibble: 3 x 3
##   Sujeto   promedio desvEst
##   <chr>       <dbl>   <dbl>
## 1 Sujeto 1    1004.    313.
## 2 Sujeto 2    1141.    161.
## 3 Sujeto 3    1280.    357.
library(dplyr)
bioequivalente  %>% 
  group_by(Periodo) %>%
  summarise(promedio = mean(AUC),
            desvEst = sd(AUC))
## # A tibble: 3 x 3
##   Periodo  promedio desvEst
##   <chr>       <dbl>   <dbl>
## 1 Semana 1    1199.    221.
## 2 Semana 2    1106.    450.
## 3 Semana 3    1120.    223.
library(ggplot2)
ggplot(bioequivalente, aes(Sujeto, AUC, 
                           col = Periodo,
                           shape = Presentacion)) +
  geom_point(size = 4)

library(ggplot2)
ggplot(bioequivalente, aes(Presentacion, AUC, 
                           col = Sujeto,
                           shape = Periodo)) +
  geom_point(size = 4)

Modelo

Modelo

El modelo para este caso es:

\[ y_{ijk} = \mu + \alpha_i + \beta_j + \gamma {k} + \epsilon_{ijkm} \\ i = 1, 2, 3. \quad j = 1, 2, 3. \quad k = 1, 2, 3. \quad \textrm{ y } \quad m = 1, 2, 3\\ \epsilon_{ijkm} \sim \mathcal{N}(0, \sigma^2) i.i.d. \]

Hipótesis a probar

Factor Presentación.

\[ H_0: \alpha_1 = \alpha_2 = \alpha_3 = 0 \\ H_1: \textrm{Alguna } \alpha \textrm{ diferente de 0} \]

Hipótesis de verificación

Bloque Sujeto.

\[ H_0: \beta_1 = \beta_2 = \beta_3 = 0 \\ H_1: \textrm{Algún } \beta \textrm{ diferente de 0} \]

Bloque Periodo.

\[ H_0: \gamma_{1} = \gamma_{2} = \gamma_{3} = 0 \\ H_1: \textrm{Algún } \gamma \textrm{ diferente de 0} \]

Modelización

modelo1 <- aov(AUC ~ Presentacion + Sujeto + Periodo, data = bioequivalente)
summary(modelo1)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## Presentacion  2 442158  221079   9.783 0.0927 .
## Sujeto        2 114264   57132   2.528 0.2834  
## Periodo       2  15000    7500   0.332 0.7508  
## Residuals     2  45196   22598                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
modelo1lm <- lm(AUC ~ Presentacion + Sujeto + Periodo, data = bioequivalente)
summary(modelo1lm)
## 
## Call:
## lm(formula = AUC ~ Presentacion + Sujeto + Periodo, data = bioequivalente)
## 
## Residuals:
##      1      2      3      4      5      6      7      8      9 
## -76.89 -17.22  94.11  94.11 -76.89 -17.22 -17.22  94.11 -76.89 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            1262.89     132.58   9.526   0.0108 *
## PresentacionB=Tableta  -510.67     122.74  -4.161   0.0532 .
## PresentacionC=Cápsula   -95.67     122.74  -0.779   0.5173  
## SujetoSujeto 2          137.67     122.74   1.122   0.3786  
## SujetoSujeto 3          276.00     122.74   2.249   0.1535  
## PeriodoSemana 2         -93.00     122.74  -0.758   0.5277  
## PeriodoSemana 3         -78.33     122.74  -0.638   0.5887  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 150.3 on 2 degrees of freedom
## Multiple R-squared:  0.9267, Adjusted R-squared:  0.7068 
## F-statistic: 4.214 on 6 and 2 DF,  p-value: 0.2042
model.tables(modelo1, type = "means", se = TRUE)
## Tables of means
## Grand mean
##          
## 1141.556 
## 
##  Presentacion 
## Presentacion
## A=Solución  B=Tableta  C=Cápsula 
##     1343.7      833.0     1248.0 
## 
##  Sujeto 
## Sujeto
## Sujeto 1 Sujeto 2 Sujeto 3 
##   1003.7   1141.3   1279.7 
## 
##  Periodo 
## Periodo
## Semana 1 Semana 2 Semana 3 
##   1198.7   1105.7   1120.3 
## 
## Standard errors for differences of means
##         Presentacion Sujeto Periodo
##                122.7  122.7   122.7
## replic.            3      3       3
model.tables(modelo1, type = "effects", se = TRUE)
## Tables of effects
## 
##  Presentacion 
## Presentacion
## A=Solución  B=Tableta  C=Cápsula 
##     202.11    -308.56     106.44 
## 
##  Sujeto 
## Sujeto
## Sujeto 1 Sujeto 2 Sujeto 3 
##  -137.89    -0.22   138.11 
## 
##  Periodo 
## Periodo
## Semana 1 Semana 2 Semana 3 
##    57.11   -35.89   -21.22 
## 
## Standard errors of effects
##         Presentacion Sujeto Periodo
##                86.79  86.79   86.79
## replic.            3      3       3

Diagnósticos del modelo.

Normalidad.

shapiro.test(residuals(modelo1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo1)
## W = 0.79873, p-value = 0.01975
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
id <- qqPlot(modelo1lm, id.n = 12, reps = 20000)

hist(residuals(modelo1), freq = FALSE,
     main = "histograma de valore residuales", las = 1,
     xlab = "Residuales",
     ylab = "Densidad")
lines(density(residuals(modelo1)), col = "red", lwd = 2)
rug(residuals(modelo1))

Homocedasticidad.

library(car)
leveneTest(AUC ~ Presentacion, data = bioequivalente)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.0256 0.9748
##        6
leveneTest(AUC ~ Periodo, data = bioequivalente)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.6433 0.5583
##        6
leveneTest(AUC ~ Sujeto, data = bioequivalente)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.1611 0.8548
##        6

Independencia

library(ggplot2)
ggplot(bioequivalente, aes(Periodo, AUC, 
                           col = Sujeto,
                           shape = Presentacion)) +
  geom_point(size = 4)

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
dwtest(modelo1lm, order.by = ~ Periodo, 
       data = bioequivalente, 
       alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  modelo1lm
## DW = 1.706, p-value = 0.2164
## alternative hypothesis: true autocorrelation is not 0

Comparaciones pareadas

TukeyHSD(modelo1, "Presentacion")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = AUC ~ Presentacion + Sujeto + Periodo, data = bioequivalente)
## 
## $Presentacion
##                            diff        lwr       upr     p adj
## B=Tableta-A=Solución -510.66667 -1233.7049  212.3715 0.0954131
## C=Cápsula-A=Solución  -95.66667  -818.7049  627.3715 0.7494600
## C=Cápsula-B=Tableta   415.00000  -308.0382 1138.0382 0.1376819