library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggpubr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(datarium)
df <- data.frame(patient=rep(1:5, each=4),
                 drug=rep(1:4, times=5),
                 response=c(30, 28, 16, 34,
                            14, 18, 10, 22,
                            24, 20, 18, 30,
                            38, 34, 20, 44,
                            26, 28, 14, 30))

df
##    patient drug response
## 1        1    1       30
## 2        1    2       28
## 3        1    3       16
## 4        1    4       34
## 5        2    1       14
## 6        2    2       18
## 7        2    3       10
## 8        2    4       22
## 9        3    1       24
## 10       3    2       20
## 11       3    3       18
## 12       3    4       30
## 13       4    1       38
## 14       4    2       34
## 15       4    3       20
## 16       4    4       44
## 17       5    1       26
## 18       5    2       28
## 19       5    3       14
## 20       5    4       30
df$patient=factor(df$patient)
df$drug=factor(df$drug)
df
##    patient drug response
## 1        1    1       30
## 2        1    2       28
## 3        1    3       16
## 4        1    4       34
## 5        2    1       14
## 6        2    2       18
## 7        2    3       10
## 8        2    4       22
## 9        3    1       24
## 10       3    2       20
## 11       3    3       18
## 12       3    4       30
## 13       4    1       38
## 14       4    2       34
## 15       4    3       20
## 16       4    4       44
## 17       5    1       26
## 18       5    2       28
## 19       5    3       14
## 20       5    4       30

Resumen Estadistico

df %>% group_by(drug) %>% get_summary_stats(response,type="mean_sd")
## # A tibble: 4 x 5
##   drug  variable     n  mean    sd
##   <fct> <chr>    <dbl> <dbl> <dbl>
## 1 1     response     5  26.4  8.76
## 2 2     response     5  25.6  6.54
## 3 3     response     5  15.6  3.85
## 4 4     response     5  32    8
df
##    patient drug response
## 1        1    1       30
## 2        1    2       28
## 3        1    3       16
## 4        1    4       34
## 5        2    1       14
## 6        2    2       18
## 7        2    3       10
## 8        2    4       22
## 9        3    1       24
## 10       3    2       20
## 11       3    3       18
## 12       3    4       30
## 13       4    1       38
## 14       4    2       34
## 15       4    3       20
## 16       4    4       44
## 17       5    1       26
## 18       5    2       28
## 19       5    3       14
## 20       5    4       30

VisualizaciĆ³n

bxp<-ggboxplot(df,x="drug",y="response", add="point")
bxp

Supuestos del Modelo

Outliers

df %>% group_by(drug) %>% identify_outliers(response)
## # A tibble: 3 x 5
##   drug  patient response is.outlier is.extreme
##   <fct> <fct>      <dbl> <lgl>      <lgl>     
## 1 1     2             14 TRUE       FALSE     
## 2 4     2             22 TRUE       FALSE     
## 3 4     4             44 TRUE       FALSE

Normalidad

df %>% group_by(drug) %>% shapiro_test(response)
## # A tibble: 4 x 4
##   drug  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 1     response     0.985 0.962
## 2 2     response     0.922 0.544
## 3 3     response     0.979 0.928
## 4 4     response     0.949 0.731
ggqqplot(df,"response",facet.by="drug")

Calculo de ANOVA

res.aov<-anova_test(data=df,dv=response,wid=patient,within=drug)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##   Effect DFn DFd      F        p p<.05   ges
## 1   drug   3  12 24.759 1.99e-05     * 0.468

Comparacion de Medias

pwc<-df %>% pairwise_t_test(response~drug,paired=TRUE,p.adjust.method="bonferroni")
pwc
## # A tibble: 6 x 10
##   .y.      group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
## * <chr>    <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
## 1 response 1      2          5     5     0.492     4 0.648 1     ns          
## 2 response 1      3          5     5     4.19      4 0.014 0.083 ns          
## 3 response 1      4          5     5    -7.48      4 0.002 0.01  *           
## 4 response 2      3          5     5     4.39      4 0.012 0.071 ns          
## 5 response 2      4          5     5    -4         4 0.016 0.097 ns          
## 6 response 3      4          5     5    -7.36      4 0.002 0.011 *

Reporte

pwc <- pwc %>% add_xy_position(x="drug")
  bxp+
  stat_pvalue_manual(pwc) +
    labs(
         subtitle = get_test_label(res.aov,detailed=TRUE),
         caption=get_pwc_label(pwc)
         )