Efecto del Viagra

Diseño experimental

Se quiere saber la influencia del uso del viagra en relación al libido que se evalúa en el individuo que la consume.

Se tomaron tres niveles de dosis de viagra (1=efecto placebo, 2=dosis baja, 3=dosis alta) y también se tomó en cosideración el grado de libido de la pareja, porque puede influenciar la respuesta del individuo que toma la pastilla.

El experimento tomo a 45 parejas y se distribuyó al azar, incialmente se propuso 15 parejas para cada nivel del experimento, pero dado que cada pareja era voluntaria varios participantes no terminaron el experimento y sólo 30 parejas lo completó.

Paquetes utilizados

library(ggplot2)
library(car)

## Loading required package: carData

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

Lectura de la base de datos

viagra <- read.table("ViagraCovariate.dat", header = TRUE)
viagra$dosis <- factor(viagra$dosis, labels = c("Placebo", "Dosis baja", "Dosis alta"))
viagra

##         dosis libido Libido_pareja
## 1     Placebo      3             4
## 2     Placebo      2             1
## 3     Placebo      5             5
## 4     Placebo      2             1
## 5     Placebo      2             2
## 6     Placebo      2             2
## 7     Placebo      7             7
## 8     Placebo      2             4
## 9     Placebo      4             5
## 10 Dosis baja      7             5
## 11 Dosis baja      5             3
## 12 Dosis baja      3             1
## 13 Dosis baja      4             2
## 14 Dosis baja      4             2
## 15 Dosis baja      7             6
## 16 Dosis baja      5             4
## 17 Dosis baja      4             2
## 18 Dosis alta      9             1
## 19 Dosis alta      2             3
## 20 Dosis alta      6             5
## 21 Dosis alta      3             4
## 22 Dosis alta      4             3
## 23 Dosis alta      4             3
## 24 Dosis alta      4             2
## 25 Dosis alta      6             0
## 26 Dosis alta      4             1
## 27 Dosis alta      6             3
## 28 Dosis alta      2             0
## 29 Dosis alta      8             1
## 30 Dosis alta      5             0

Análisis exploratorio

library(ggplot2)
g1 <- ggplot(viagra, aes(dosis, libido))
g1 + geom_jitter(width = 0.05) +
     stat_summary(fun.y = "mean", geom = "point", col = "red", shape = 15, size = 3)

Ajuste de un modelo incial

\[ y_i = \mu + \alpha_i + \epsilon_{ij} \\ \epsilon_{ij} \sim \mathcal{N}(0, \sigma^2) \quad i.i.d. \]

m1 <- aov(libido ~ dosis, data = viagra)
anova(m1)

## Analysis of Variance Table
## 
## Response: libido
##           Df Sum Sq Mean Sq F value Pr(>F)
## dosis      2 16.844  8.4219  2.4159 0.1083
## Residuals 27 94.123  3.4860

shapiro.test(residuals(m1))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m1)
## W = 0.9292, p-value = 0.04676

shapiro.test(rstandard(m1))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(m1)
## W = 0.92873, p-value = 0.04543

shapiro.test(rstudent(m1))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(m1)
## W = 0.92518, p-value = 0.03662

library(car)
leveneTest(m1)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.3256 0.7249
##       27

tk1 <- TukeyHSD(m1)

op <- par(no.readonly = TRUE)
par(mar = op$mar + c(0, 5, 0, 0))
plot(tk1, las = 1)

par(op)

Exploratorio con la covariable

library(ggplot2)
g2 <- ggplot(viagra, aes(Libido_pareja, libido,  col = dosis))
g2 + geom_jitter(width = 0.05, height = 0.05) +
     geom_smooth(method = "lm", aes(fill = dosis), alpha = 0.1) +
     scale_color_manual(values = c("blue", "darkgreen", "red")) +
     scale_fill_manual(values = c("blue", "darkgreen", "red")) 

Estimación de modelo con covariable

Nuevo modelo con covariable \[ y_i = \beta_0 + \beta_1 D_{2} + \beta_2 D_{3} + \beta_3 x_i + \epsilon_{ij} \\ \epsilon_{ij} \sim \mathcal{N}(0, \sigma^2) \quad i.i.d. \] Donde:

• \(D_2\) = 1 si dosis es baja, de lo contrario 0.
• \(D_3\) = 1 si dosis es alta, de lo contrario 0.
• Si tanto \(D_2\) como \(D_3\) son cero entonces se trata del placebo.
• \(x\) es la covariable líbido de la pareja.
• \(\epsilon\) = Error aleatorio.

library(car)
m2 <- aov(libido ~ dosis + Libido_pareja, data = viagra)
Anova(m2)

## Anova Table (Type II tests)
## 
## Response: libido
##               Sum Sq Df F value  Pr(>F)  
## dosis         25.185  2  4.1419 0.02745 *
## Libido_pareja 15.076  1  4.9587 0.03483 *
## Residuals     79.047 26                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m2l <- lm(libido ~ dosis + Libido_pareja, data = viagra)
summary(m2l)

## 
## Call:
## lm(formula = libido ~ dosis + Libido_pareja, data = viagra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2622 -0.7899 -0.3230  0.8811  4.5699 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       1.7892     0.8671   2.063   0.0492 *
## dosisDosis baja   1.7857     0.8494   2.102   0.0454 *
## dosisDosis alta   2.2249     0.8028   2.771   0.0102 *
## Libido_pareja     0.4160     0.1868   2.227   0.0348 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.744 on 26 degrees of freedom
## Multiple R-squared:  0.2876, Adjusted R-squared:  0.2055 
## F-statistic:   3.5 on 3 and 26 DF,  p-value: 0.02954

Validación de supuestos

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(m2l)

par(op)

shapiro.test(residuals(m2l))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m2l)
## W = 0.94327, p-value = 0.1114

shapiro.test(rstandard(m2l))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(m2l)
## W = 0.9446, p-value = 0.121

shapiro.test(rstudent(m2l))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(m2l)
## W = 0.92734, p-value = 0.04175

library(car)
library(lmtest)

bartlett.test(libido ~ dosis, data = viagra)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  libido by dosis
## Bartlett's K-squared = 1.1104, df = 2, p-value = 0.574

car::ncvTest(m2l)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.153817    Df = 1     p = 0.2827515

lmtest::bptest(m2l, ~ Libido_pareja, studentize = FALSE, data = viagra)

## 
##  Breusch-Pagan test
## 
## data:  m2l
## BP = 1.046, df = 1, p-value = 0.3064

lmtest::bptest(m2l, ~ Libido_pareja, studentize = TRUE, data = viagra)

## 
##  studentized Breusch-Pagan test
## 
## data:  m2l
## BP = 0.6881, df = 1, p-value = 0.4068

library(car)
car::influenceIndexPlot(m2l)

library(car)
car::influencePlot(m2l)

##      StudRes       Hat      CookD
## 7  1.5711252 0.2562590 0.20126041
## 15 0.5956225 0.2199011 0.02563739
## 18 3.1940828 0.0884045 0.18268737
## 29 2.3176228 0.0884045 0.11148247

Modelo final con covariable.

viagra <- data.frame(viagra, predict(m2l, interval = "confidence"))

g3 <- ggplot(viagra, aes(Libido_pareja, libido, col = dosis, label = rownames(viagra)))
g3 + geom_line(aes(y = fit)) +
     geom_ribbon(aes(ymin = lwr, ymax = upr, fill = dosis), alpha = 0.1) +
     scale_color_manual(values = c("blue", "darkgreen", "red")) +
     scale_fill_manual(values = c("blue", "darkgreen", "red")) +
     geom_text(size = 3)

Estimación del modelo con covariable e interacción.

Nuevo modelo con covariable \[ y_i = \beta_0 + \beta_1 D_{2} + \beta_2 D_{3} + \beta_3 x_i + \beta_4 D_{2} x_i + \beta_5 D_{3} x_i + \epsilon_{ij} \\ \epsilon_{ij} \sim \mathcal{N}(0, \sigma^2) \quad i.i.d. \] Donde:

• \(D_2\) = 1 si dosis es baja, de lo contrario 0.
• \(D_3\) = 1 si dosis es alta, de lo contrario 0.
• Si tanto \(D_2\) como \(D_3\) son cero entonces se trata del placebo.
• \(x\) es la covariable líbido de la pareja.
• \(\epsilon\) = Error aleatorio.

library(car)
m3 <- aov(libido ~ dosis + Libido_pareja + dosis:Libido_pareja, data = viagra)
Anova(m3)

## Anova Table (Type II tests)
## 
## Response: libido
##                     Sum Sq Df F value  Pr(>F)  
## dosis               25.185  2  5.1556 0.01372 *
## Libido_pareja       15.076  1  6.1722 0.02035 *
## dosis:Libido_pareja 20.427  2  4.1815 0.02767 *
## Residuals           58.621 24                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m3l <- lm(libido ~ dosis + Libido_pareja + dosis:Libido_pareja, data = viagra)
summary(m3l)

## 
## Call:
## lm(formula = libido ~ dosis + Libido_pareja + dosis:Libido_pareja, 
##     data = viagra)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2837 -0.6274 -0.1201  0.6299  3.9351 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                    0.59416    1.05743   0.562  0.57940   
## dosisDosis baja                1.71722    1.60191   1.072  0.29439   
## dosisDosis alta                4.68950    1.26940   3.694  0.00114 **
## Libido_pareja                  0.76299    0.26716   2.856  0.00872 **
## dosisDosis baja:Libido_pareja  0.05737    0.43403   0.132  0.89594   
## dosisDosis alta:Libido_pareja -0.98174    0.38432  -2.554  0.01740 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.563 on 24 degrees of freedom
## Multiple R-squared:  0.4717, Adjusted R-squared:  0.3617 
## F-statistic: 4.286 on 5 and 24 DF,  p-value: 0.006295

Validación de supuestos

op <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(m3l)

par(op)

shapiro.test(residuals(m3l))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(m3l)
## W = 0.94754, p-value = 0.1453

shapiro.test(rstandard(m3l))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(m3l)
## W = 0.95489, p-value = 0.2281

shapiro.test(rstudent(m3l))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(m3l)
## W = 0.93528, p-value = 0.06793

library(car)
library(lmtest)

bartlett.test(libido ~ dosis, data = viagra)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  libido by dosis
## Bartlett's K-squared = 1.1104, df = 2, p-value = 0.574

car::ncvTest(m3l)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 2.057503    Df = 1     p = 0.1514582

lmtest::bptest(m3l, ~ Libido_pareja,  studentize = FALSE, data = viagra)

## 
##  Breusch-Pagan test
## 
## data:  m3l
## BP = 3.6786, df = 1, p-value = 0.05512

lmtest::bptest(m3l, ~ Libido_pareja,  studentize = TRUE, data = viagra)

## 
##  studentized Breusch-Pagan test
## 
## data:  m3l
## BP = 2.1416, df = 1, p-value = 0.1434

library(car)
car::influenceIndexPlot(m3l)

library(car)
car::influencePlot(m3l)

##       StudRes       Hat       CookD
## 7   0.9432328 0.4805195 0.137793595
## 15 -0.2115547 0.5209581 0.008448164
## 18  3.1112043 0.1081731 0.143707367
## 20  1.4812086 0.3581731 0.194388549
## 28 -2.6245943 0.2019231 0.233249917

Modelo final con covariable e interacción.

viagra <- data.frame(viagra, predict(m3l, interval = "confidence"))

g3 <- ggplot(viagra, aes(Libido_pareja, libido, col = dosis, label = rownames(viagra)))
g3 + geom_line(aes(y = fit.1)) +
     geom_ribbon(aes(ymin = lwr.1, ymax = upr.1, fill = dosis), alpha = 0.1) +
     scale_color_manual(values = c("blue", "darkgreen", "red")) +
     scale_fill_manual(values = c("blue", "darkgreen", "red")) +
     geom_text(size = 3)

Consecuencia de quitar valores influenciales.

m4l <- lm(libido ~ dosis + Libido_pareja + dosis:Libido_pareja, data = viagra,
          subset = -c(28))
car::Anova(m4l)

## Anova Table (Type II tests)
## 
## Response: libido
##                     Sum Sq Df F value   Pr(>F)   
## dosis               26.732  2  6.8148 0.004740 **
## Libido_pareja       10.925  1  5.5701 0.027128 * 
## dosis:Libido_pareja 29.312  2  7.4727 0.003159 **
## Residuals           45.110 23                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m4l)

## 
## Call:
## lm(formula = libido ~ dosis + Libido_pareja + dosis:Libido_pareja, 
##     data = viagra, subset = -c(28))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6867 -0.6867 -0.1201  0.5909  3.3614 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    0.59416    0.94756   0.627 0.536810    
## dosisDosis baja                1.71722    1.43547   1.196 0.243774    
## dosisDosis alta                5.52030    1.18072   4.675 0.000105 ***
## Libido_pareja                  0.76299    0.23940   3.187 0.004103 ** 
## dosisDosis baja:Libido_pareja  0.05737    0.38893   0.148 0.884012    
## dosisDosis alta:Libido_pareja -1.23889    0.35805  -3.460 0.002125 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.4 on 23 degrees of freedom
## Multiple R-squared:  0.5711, Adjusted R-squared:  0.4778 
## F-statistic: 6.125 on 5 and 23 DF,  p-value: 0.0009554

library(car)
influenceIndexPlot(m4l)