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
library(ggplot2)
library(ggpubr)
library(dlookr)
## Registered S3 methods overwritten by 'dlookr':
##   method          from  
##   plot.transform  scales
##   print.transform scales
## 
## Attaching package: 'dlookr'
## The following object is masked from 'package:base':
## 
##     transform
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(ggsignif)
library(dplyr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
#valores de pearsonCC para vehiculo (30 observacion)
pearson_vehicle <-c(
   0.69, 0.67, 0.81, 0.78, 0.68, 0.80, 0.91, 0.88, 0.88, 0.79,
  0.76, 0.64, 0.84, 0.79, 0.74, 0.92, 0.74, 0.71, 0.87, 0.64,
  0.84, 0.79, 0.74, 0.92, 0.74, 0.71, 0.87, 0.79, 0.80, 0.78)
length(pearson_vehicle)
## [1] 30
# Valores de PearsonCC para metformina (30 observaciones)
pearson_metfotmina <-c (
  0.16, 0.21, 0.31, 0.47, 0.68, 0.69, 0.37, 0.48, 0.33, 0.54,
  0.43, 0.25, 0.48, 0.35, 0.43, 0.59, 0.23, 0.37, 0.20, 0.59,
  0.36, 0.27, 0.24, 0.48, 0.46, 0.09, 0.57, 0.50, 0.28, 0.39
)

length(pearson_metfotmina)
## [1] 30
mett <-data.frame(
  Treatment = factor (
    c(rep("vehicle", length(pearson_vehicle)), rep("metformin", length(pearson_metfotmina))),
    levels = c("vehicle", "metformin")
    ), 
  PearsonCC = c(pearson_vehicle, pearson_metfotmina)
)
str(mett)
## 'data.frame':    60 obs. of  2 variables:
##  $ Treatment: Factor w/ 2 levels "vehicle","metformin": 1 1 1 1 1 1 1 1 1 1 ...
##  $ PearsonCC: num  0.69 0.67 0.81 0.78 0.68 0.8 0.91 0.88 0.88 0.79 ...
head(mett)
##   Treatment PearsonCC
## 1   vehicle      0.69
## 2   vehicle      0.67
## 3   vehicle      0.81
## 4   vehicle      0.78
## 5   vehicle      0.68
## 6   vehicle      0.80
nrow(mett)
## [1] 60
lower_ci <- function(mean, se, n, conf_level = 0.95) {
  mean - qt(1 - (1 - conf_level) / 2, df = n - 1) * se
}

upper_ci <- function(mean, se, n, conf_level = 0.95) {
  mean + qt(1 - (1 - conf_level) / 2, df = n - 1) * se
}
desc_mett <- mett %>%
  group_by(Treatment) %>%
  summarise(
    n    = n(),
    mean = mean(PearsonCC, na.rm = TRUE),
    sd   = sd(PearsonCC, na.rm = TRUE),
    .groups = "drop"
  ) %>%
 mutate(se = sd/sqrt(n),
       lower_ci = lower_ci(mean, se, n),
       upper_ci = upper_ci(mean, se, n))

desc_mett
## # A tibble: 2 × 7
##   Treatment     n  mean     sd     se lower_ci upper_ci
##   <fct>     <int> <dbl>  <dbl>  <dbl>    <dbl>    <dbl>
## 1 vehicle      30 0.784 0.0803 0.0147    0.754    0.814
## 2 metformin    30 0.393 0.153  0.0279    0.336    0.450
with(mett, shapiro.test(PearsonCC[Treatment == "vehicle"]))
## 
##  Shapiro-Wilk normality test
## 
## data:  PearsonCC[Treatment == "vehicle"]
## W = 0.96534, p-value = 0.4207
with(mett, shapiro.test(PearsonCC[Treatment == "metformin"]))
## 
##  Shapiro-Wilk normality test
## 
## data:  PearsonCC[Treatment == "metformin"]
## W = 0.98343, p-value = 0.9075

mett %>%
  group_by(Treatment) %>%
  plot_normality(PearsonCC)
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the dlookr package.
##   Please report the issue at <https://github.com/choonghyunryu/dlookr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

res_ftest <- var.test(PearsonCC ~ Treatment, data = mett)
res_ftest
## 
##  F test to compare two variances
## 
## data:  PearsonCC by Treatment
## F = 0.2758, num df = 29, denom df = 29, p-value = 0.0008667
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1312703 0.5794512
## sample estimates:
## ratio of variances 
##          0.2757984
# Ajustar var.equal según el resultado de la prueba F
res_t_mett <- t.test(
  PearsonCC ~ Treatment,
  data       = mett,
  alternative = "two.sided",
  var.equal   = FALSE   # modifique a TRUE si las varianzas son homogéneas
)

res_t_mett
## 
##  Welch Two Sample t-test
## 
## data:  PearsonCC by Treatment
## t = 12.385, df = 43.866, p-value = 6.455e-16
## alternative hypothesis: true difference in means between group vehicle and group metformin is not equal to 0
## 95 percent confidence interval:
##  0.3270915 0.4542418
## sample estimates:
##   mean in group vehicle mean in group metformin 
##               0.7840000               0.3933333
p_value_t <- res_t_mett$p.value
means_t   <- res_t_mett$estimate
ci_t      <- res_t_mett$conf.int

p_value_t
## [1] 6.455027e-16
means_t
##   mean in group vehicle mean in group metformin 
##               0.7840000               0.3933333
ci_t
## [1] 0.3270915 0.4542418
## attr(,"conf.level")
## [1] 0.95
max_y_mett <- max(mett$PearsonCC, na.rm = TRUE)

boxplot_mett <- ggboxplot(
  mett,
  x      = "Treatment",
  y      = "PearsonCC",
  color  = "Treatment",
  palette = c("#5A0E24", "#BF124D"),
  ylab   = "Índice de colocalización de Pearson (E vs NS4A)",
  xlab   = "Tratamiento",
  add    = "jitter"
)

boxplot_mett +
  stat_compare_means(
    method  = "t.test",
    label.y = max_y_mett * 1.05
  ) +
  stat_compare_means(
    aes(label = ..p.signif..),
    method  = "t.test",
    label.y = max_y_mett * 1.10
  )

barplot_mett <- ggbarplot(
  mett,
  x       = "Treatment",
  y       = "PearsonCC",
  color   = "Treatment",
  fill    = "Treatment",          # usar la variable de grupo como relleno
  palette = c("#5A0E24", "#BF124D"),
  ylab    = "Índice de colocalización de Pearson (E vs NS4A)",
  xlab    = "Tratamiento",
  add     = "mean_sd"             # barras de error con media ± SD
)

barplot_mett +
  stat_compare_means(
    method  = "t.test",
    label.y = max_y_mett * 1.05
  ) +
  stat_compare_means(
    aes(label = ..p.signif..),
    method  = "t.test",
    label.y = max_y_mett * 1.10
  )

PRACTICA 2

nup98 <- data.frame(
  nup98 = factor(
    c(rep("Zika24h", 3),
      rep("NS2B3Mut", 3),
      rep("NS2B3WT", 3)),
    levels  = c("Zika24h", "NS2B3Mut", "NS2B3WT"),
    ordered = TRUE
  ),
  RBI = c(
    0.547, 0.52, 0.65,
    0.988, 1.09, 0.927,
    0.72, 0.739, 0.82
  )
)

str(nup98)
## 'data.frame':    9 obs. of  2 variables:
##  $ nup98: Ord.factor w/ 3 levels "Zika24h"<"NS2B3Mut"<..: 1 1 1 2 2 2 3 3 3
##  $ RBI  : num  0.547 0.52 0.65 0.988 1.09 0.927 0.72 0.739 0.82
nup98
##      nup98   RBI
## 1  Zika24h 0.547
## 2  Zika24h 0.520
## 3  Zika24h 0.650
## 4 NS2B3Mut 0.988
## 5 NS2B3Mut 1.090
## 6 NS2B3Mut 0.927
## 7  NS2B3WT 0.720
## 8  NS2B3WT 0.739
## 9  NS2B3WT 0.820
desc_nup98 <- nup98 %>%
  group_by(nup98) %>%
  summarise(
    n    = n(),
    mean = mean(RBI, na.rm = TRUE),
    sd   = sd(RBI, na.rm = TRUE),
    .groups = "drop"
  )

desc_nup98
## # A tibble: 3 × 4
##   nup98        n  mean     sd
##   <ord>    <int> <dbl>  <dbl>
## 1 Zika24h      3 0.572 0.0686
## 2 NS2B3Mut     3 1.00  0.0824
## 3 NS2B3WT      3 0.760 0.0531
resultado_anova <- aov(RBI ~ nup98, data = nup98)
summary(resultado_anova)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## nup98        2 0.27798 0.13899   29.14 0.000813 ***
## Residuals    6 0.02862 0.00477                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(resultado_anova, which = 1)

leveneTest(RBI ~ nup98, data = nup98)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.1364 0.8751
##        6
resid_anova <- residuals(resultado_anova)

# Q-Q plot
plot(resultado_anova, which = 2)

# Prueba de Shapiro–Wilk sobre los residuos
shapiro.test(resid_anova)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_anova
## W = 0.88378, p-value = 0.172
summary(resultado_anova)
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## nup98        2 0.27798 0.13899   29.14 0.000813 ***
## Residuals    6 0.02862 0.00477                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(resultado_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = RBI ~ nup98, data = nup98)
## 
## $nup98
##                        diff         lwr         upr     p adj
## NS2B3Mut-Zika24h  0.4293333  0.25631482  0.60235185 0.0006545
## NS2B3WT-Zika24h   0.1873333  0.01431482  0.36035185 0.0367107
## NS2B3WT-NS2B3Mut -0.2420000 -0.41501852 -0.06898148 0.0121671
summary(
  glht(
    resultado_anova,
    linfct = mcp(nup98 = "Tukey")
  )
)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = RBI ~ nup98, data = nup98)
## 
## Linear Hypotheses:
##                         Estimate Std. Error t value Pr(>|t|)    
## NS2B3Mut - Zika24h == 0  0.42933    0.05639   7.614   <0.001 ***
## NS2B3WT - Zika24h == 0   0.18733    0.05639   3.322   0.0369 *  
## NS2B3WT - NS2B3Mut == 0 -0.24200    0.05639  -4.292   0.0123 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
max_y_nup <- max(nup98$RBI, na.rm = TRUE)

boxplot_nup <- ggboxplot(
  nup98,
  x      = "nup98",
  y      = "RBI",
  color  = "nup98",                         # color por grupo
  fill   = "nup98",                         # relleno por grupo
  palette = c("#1B211A", "#628141", "#8BAE66"),
  order  = c("Zika24h", "NS2B3Mut", "NS2B3WT"),
  ylab   = "Intensidad relativa (unidades arbitrarias)",
  xlab   = "Condiciones",
  add    = "jitter"
)

boxplot_nup +
  stat_compare_means(
    method  = "anova",
    label.y = max_y_nup * 1.05
  )

my_comparisons <- list(
  c("Zika24h",  "NS2B3Mut"),
  c("Zika24h",  "NS2B3WT"),
  c("NS2B3Mut", "NS2B3WT")
)
ggplot(nup98, aes(x = nup98, y = RBI)) +
  geom_boxplot() +
  geom_signif(
    comparisons      = my_comparisons,
    map_signif_level = TRUE,
    step_increase    = 0.1
  ) +
  ylab("Intensidad relativa (unidades arbitrarias)") +
  xlab("Condiciones")

barplot_nup <- ggbarplot(
  nup98,
  x       = "nup98",
  y       = "RBI",
  color   = "nup98",                         # contorno por grupo
  fill    = "nup98",                         # relleno por grupo
  palette = c("#540863", "#92487A", "#E49BA6"),
  order   = c("Zika24h", "NS2B3Mut", "NS2B3WT"),
  ylab    = "Intensidad relativa (unidades arbitrarias)",
  xlab    = "Condiciones",
  add     = "mean_sd"
)

max_y_nup <- max(nup98$RBI, na.rm = TRUE)

barplot_nup +
  stat_compare_means(
    method  = "anova",
    label.y = max_y_nup * 1.05
  )

 ggline(
  nup98,
  x     = "nup98",
  y     = "RBI",
  add   = c("mean_sd", "jitter"),
  order = c("Zika24h", "NS2B3Mut", "NS2B3WT"),
  ylab  = "Intensidad relativa (unidades arbitrarias)",
  xlab  = "Condiciones"
)