#### Cleaning data for new and clean data set along with calculating age and erasing patient info after
el <- EL_DATA %>%
  rename(
    DOS = DOS,
    Name = NAME,
    DOB = DOB,
    Procedure = PROCEDURE,
    Eye = EYE,
    Target = `TARGET (DISTANCE or NEAR)`,
    Pre_Sphere = `SPHERE (pre-op)`,
    Pre_Cyl = `CYL (pre-op)`,
    Pre_Axis = `AXIS (pre-op)`,
    Lasik_Type = `TYPE (Contoura vs WFO if LASIK)`,
    POD1_VA = `POD1 VA`,
    POW1_VA = `POW1 VA`,
    POM1_VA = `POM1 or 2 VA`,
    Post_Sphere = `SPHERE (post-op)`,
    Post_Cyl = `CYL (post-op)`,
    Post_Axis = `AXIS (post-op)`,
    Enhancement = `ENH? (YES or NO)`,
    Enh_Sphere = `ENH SPH`,
    Enh_Cyl = `ENH CYL`,
    Enh_Axis = `ENH AXIS`,
    ICL_Exchange_Reason = `Reason for Exchange (ICL)`
  )

el <- el %>%
  mutate(
    DOS = mdy(DOS),
    DOB = mdy(DOB),
    Age = as.numeric(DOS - DOB) / 365.25
  ) %>%
  select(-Name, -DOB, -DOS)
#### further cleaning more focused on prescription section will also be showing where NA values or miss inputs may be hidden 
el <- el %>%
  mutate(
    Pre_Cyl = ifelse(Pre_Cyl %in% c("Sph", "SPH", "sph", "Plano", "plano"), 0, Pre_Cyl),
    Post_Cyl = ifelse(Post_Cyl %in% c("Sph", "SPH", "sph", "Plano", "plano"), 0, Post_Cyl),
    
    Pre_Sphere = as.numeric(Pre_Sphere),
    Pre_Cyl = as.numeric(Pre_Cyl),
    Pre_Axis = as.numeric(Pre_Axis),
    
    Post_Sphere = as.numeric(Post_Sphere),
    Post_Cyl = as.numeric(Post_Cyl),
    Post_Axis = as.numeric(Post_Axis)
  )

EL_DATA %>%
  filter(is.na(as.numeric(`SPHERE (pre-op)`))) %>%
  select(`SPHERE (pre-op)`) %>%
  distinct()
## # A tibble: 3 × 1
##   `SPHERE (pre-op)`
##   <chr>            
## 1 <NA>             
## 2 -3/75            
## 3 --9.00
### stats for calculations residual
el <- el %>%
  mutate(
    Pre_SE = Pre_Sphere + (Pre_Cyl / 2),
    Post_SE = Post_Sphere + (Post_Cyl / 2),
    Residual_SE = abs(Post_SE),
    Change_SE = Post_SE - Pre_SE,
    Abs_Cyl_Post = abs(Post_Cyl)
  )
#### categorical variables 

el <- el %>%
  mutate(
    Procedure = str_trim(Procedure),
    Procedure = case_when(
      str_detect(Procedure, regex("EVO", ignore_case = TRUE)) ~ "EVO",
      str_detect(Procedure, regex("LASIK", ignore_case = TRUE)) ~ "LASIK",
      TRUE ~ Procedure
    ),
    
    Lasik_Type = str_trim(Lasik_Type),
    Lasik_Type = case_when(
      str_detect(Lasik_Type, regex("contoura", ignore_case = TRUE)) ~ "Contoura",
      str_detect(Lasik_Type, regex("WFO", ignore_case = TRUE)) ~ "WFO",
      TRUE ~ Lasik_Type
    ),
    
    Enhancement = str_trim(Enhancement),
    Enhancement = case_when(
      str_detect(Enhancement, regex("^yes$|^y$", ignore_case = TRUE)) ~ "Yes",
      str_detect(Enhancement, regex("^no$|^n$", ignore_case = TRUE)) ~ "No",
      TRUE ~ NA_character_
    ),
    
    Procedure = factor(Procedure),
    Lasik_Type = factor(Lasik_Type),
    Enhancement = factor(Enhancement, levels = c("No", "Yes"))
  )
### scope of NA values as suspected 
####Abs_Cyl_Post = amount of residual astigmatism  Change_SE = how much the patient's refractive error changed because of surgery 
####Residual_SE = absolute residual error Spherical Equivalent (SE) Pre_SE = (Pre-operative Spherical Equivalent) Post_SE = Post-operative Spherical Equivalent

table(el$Procedure, useNA = "ifany")
## 
##   EVO LASIK 
##   282   286
table(el$Eye, useNA = "ifany")
## 
##   OD   OS <NA> 
##  288  277    3
table(el$Procedure, el$Eye, useNA = "ifany")
##        
##          OD  OS <NA>
##   EVO   145 136    1
##   LASIK 143 141    2
table(el$Lasik_Type, useNA = "ifany")
## 
## Contoura      WFO     <NA> 
##       68      216      284
table(el$Enhancement, useNA = "ifany")
## 
##   No  Yes <NA> 
##  525    6   37
summary_by_procedure <- el %>%
  group_by(Procedure) %>%
  summarise(
    n = n(),
    Mean_Age = mean(Age, na.rm = TRUE),
    SD_Age = sd(Age, na.rm = TRUE),
    Mean_Pre_SE = mean(Pre_SE, na.rm = TRUE),
    SD_Pre_SE = sd(Pre_SE, na.rm = TRUE),
    Mean_Post_SE = mean(Post_SE, na.rm = TRUE),
    SD_Post_SE = sd(Post_SE, na.rm = TRUE),
    Mean_Residual_SE = mean(Residual_SE, na.rm = TRUE),
    SD_Residual_SE = sd(Residual_SE, na.rm = TRUE),
    Mean_Post_Cyl = mean(Abs_Cyl_Post, na.rm = TRUE),
    SD_Post_Cyl = sd(Abs_Cyl_Post, na.rm = TRUE)
  )

summary_by_procedure
## # A tibble: 2 × 12
##   Procedure     n Mean_Age SD_Age Mean_Pre_SE SD_Pre_SE Mean_Post_SE SD_Post_SE
##   <fct>     <int>    <dbl>  <dbl>       <dbl>     <dbl>        <dbl>      <dbl>
## 1 EVO         282     34.1  11.0        -6.69      2.92       -0.216      0.512
## 2 LASIK       286     33.8   9.17       -2.70      1.85       -0.167      0.490
## # ℹ 4 more variables: Mean_Residual_SE <dbl>, SD_Residual_SE <dbl>,
## #   Mean_Post_Cyl <dbl>, SD_Post_Cyl <dbl>
#### Compare residual refractive error t test
t.test(Residual_SE ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  Residual_SE by Procedure
## t = 0.91321, df = 427.7, p-value = 0.3616
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -0.04003627  0.10952378
## sample estimates:
##   mean in group EVO mean in group LASIK 
##           0.3796512           0.3449074
wilcox.test(Residual_SE ~ Procedure, data = el)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Residual_SE by Procedure
## W = 24346, p-value = 0.3772
## alternative hypothesis: true location shift is not equal to 0

Wilcoxon rank-sum test is a nonparametric alternative that compares the distributions without assuming normality.

#### LASIK comparison 
lasik_only <- el %>%
  filter(Procedure == "LASIK",
         Lasik_Type %in% c("Contoura", "WFO"))

table(lasik_only$Lasik_Type, useNA = "ifany")
## 
## Contoura      WFO 
##       68      216
lasik_summary <- lasik_only %>%
  group_by(Lasik_Type) %>%
  summarise(
    n = n(),
    Mean_Pre_SE = mean(Pre_SE, na.rm = TRUE),
    SD_Pre_SE = sd(Pre_SE, na.rm = TRUE),
    Mean_Post_SE = mean(Post_SE, na.rm = TRUE),
    SD_Post_SE = sd(Post_SE, na.rm = TRUE),
    Mean_Residual_SE = mean(Residual_SE, na.rm = TRUE),
    SD_Residual_SE = sd(Residual_SE, na.rm = TRUE),
    Mean_Post_Cyl = mean(Abs_Cyl_Post, na.rm = TRUE),
    SD_Post_Cyl = sd(Abs_Cyl_Post, na.rm = TRUE)
  )

lasik_summary
## # A tibble: 2 × 10
##   Lasik_Type     n Mean_Pre_SE SD_Pre_SE Mean_Post_SE SD_Post_SE
##   <fct>      <int>       <dbl>     <dbl>        <dbl>      <dbl>
## 1 Contoura      68       -3.44      1.54       -0.152      0.418
## 2 WFO          216       -2.47      1.89       -0.171      0.510
## # ℹ 4 more variables: Mean_Residual_SE <dbl>, SD_Residual_SE <dbl>,
## #   Mean_Post_Cyl <dbl>, SD_Post_Cyl <dbl>
t.test(Residual_SE ~ Lasik_Type, data = lasik_only)
## 
##  Welch Two Sample t-test
## 
## data:  Residual_SE by Lasik_Type
## t = -0.42911, df = 109.25, p-value = 0.6687
## alternative hypothesis: true difference in means between group Contoura and group WFO is not equal to 0
## 95 percent confidence interval:
##  -0.12726679  0.08196559
## sample estimates:
## mean in group Contoura      mean in group WFO 
##              0.3275000              0.3501506
wilcox.test(Residual_SE ~ Lasik_Type, data = lasik_only)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Residual_SE by Lasik_Type
## W = 4354, p-value = 0.592
## alternative hypothesis: true location shift is not equal to 0
#### Plots 
ggplot(el, aes(x = Procedure, y = Pre_SE, fill = Procedure)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(
    title = "Pre-operative Spherical Equivalent by Procedure",
    x = "Procedure",
    y = "Pre-op Spherical Equivalent (D)"
  ) +
  theme_minimal()

ggplot(el, aes(x = Procedure, y = Residual_SE, fill = Procedure)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(
    title = "Residual Refractive Error by Procedure",
    x = "Procedure",
    y = "Absolute Post-op SE (D)"
  ) +
  theme_minimal()

ggplot(lasik_only, aes(x = Lasik_Type, y = Residual_SE, fill = Lasik_Type)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(
    title = "Residual SE: Contoura vs WFO LASIK",
    x = "LASIK Type",
    y = "Absolute Post-op SE (D)"
  ) +
  theme_minimal()

ggplot(el, aes(x = Pre_SE, y = Post_SE, color = Procedure)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(
    title = "Pre-op SE vs Post-op SE",
    x = "Pre-op Spherical Equivalent (D)",
    y = "Post-op Spherical Equivalent (D)"
  ) +
  theme_minimal()

unique(el$POD1_VA)
##  [1] "20/20"     "20/20-1"   "20/30-1"   "20/25-1"   "20/40"     "20/15-2"  
##  [7] "20/30-2"   "20/15-1"   "20/20+2"   "20/20-2"   "20/15"     "20/12.5"  
## [13] "20/25+2"   "20/15+2"   "20/25"     "20/25+1"   "20/12.5-1" "20/20+1"  
## [19] NA          "20/30"     "20/40+1"   "20/50"     "20/40+2"   "20/30+2"  
## [25] "20/50+2"   "20/25-2"   "20/12.5+1" "20/12.5-2" "20/15+"    "20/15+1"  
## [31] "20/100"    "20/15-21"  "20/40-2"   "20/30+1"   "20/12.5+2"
unique(el$POW1_VA)
##  [1] NA          "20/40"     "20/25"     "20/20"     "20/30"     "20/20-1"  
##  [7] "20/15+1"   "20?15+2"   "20/12.5-1" "20/15-1"   "20/12.5"   "20/20+1"  
## [13] "20/12.5-2" "20/30-1"   "20/30+2"   "20/15-2"   "20/15"     "20/20+2"  
## [19] "20/25-1"   "20/20-2"   "20/30-2"   "20/25-2"   "20/15+2"   "10/15-2"  
## [25] "20/"       "20/25+2"   "20/12.5+1" "20/30+1"   "20/40+2"   "20/25+1"  
## [31] "20/50"     "20/60-1"
unique(el$POM1_VA)
##  [1] "20/20"     "20/20-2"   "20/40"     "20/25-2"   NA          "20/12.5-2"
##  [7] "20/15"     "20/20-1"   "20/25-1"   "20/25"     "20/15-1"   "20/15-2"  
## [13] "20/15+1"   "20/20+2"   "20/25+2"   "20/30-2"   "20/12.5"   "20/12.5-1"
## [19] "20/30-1"   "20/12.5+1" "20/30"     "20/20+1"   "20/15+2"   "20/60"    
## [25] "20/40-"    "20/30+2"   "20/-15-2"  "20/12.5+2" "20/25+1"
#fix errors on data set
el <- el %>%
  mutate(
    POD1_VA = case_when(
      POD1_VA == "20/15-21" ~ "20/15-2",
      TRUE ~ POD1_VA
    ),

    POW1_VA = case_when(
      POW1_VA == "20?15+2" ~ "20/15+2",
      POW1_VA == "10/15-2" ~ "20/15-2",
      POW1_VA == "20/" ~ NA_character_,
      TRUE ~ POW1_VA
    ),

    POM1_VA = case_when(
      POM1_VA == "20/-15-2" ~ "20/15-2",
      POM1_VA == "20/40-" ~ "20/40",
      TRUE ~ POM1_VA
    )
  )
#conversion to logmar
convert_va_logmar <- function(x) {
  
  # Convert to character and remove spaces
  x <- as.character(x)
  x <- str_trim(x)
  x <- str_replace_all(x, " ", "")
  
  # Keep missing values as missing
  x[x == ""] <- NA_character_
  
  # Extract numerator
  numerator <- str_extract(x, "^[0-9]+\\.?[0-9]*")
  
  # Extract denominator after /
  denominator <- str_extract(
    x,
    "(?<=/)[0-9]+\\.?[0-9]*"
  )
  
  numerator <- as.numeric(numerator)
  denominator <- as.numeric(denominator)
  
  # Base Snellen logMAR
  base_logmar <- log10(denominator / numerator)
  
  # Extract letter adjustment at end of value
  adjustment_text <- str_extract(
    x,
    "[+-][0-9]+$"
  )
  
  adjustment <- as.numeric(adjustment_text)
  
  # Values without + or - letters receive adjustment of 0
  adjustment[is.na(adjustment)] <- 0
  
  # Plus letters improve vision, so subtract 0.02 per letter
  # Minus letters worsen vision, so subtracting a negative adds 0.02
  logmar <- base_logmar - (adjustment * 0.02)
  
  return(logmar)
}
#conversion to logmar
el <- el %>%
  mutate(
    POD1_LogMAR = convert_va_logmar(POD1_VA),
    POW1_LogMAR = convert_va_logmar(POW1_VA),
    POM1_LogMAR = convert_va_logmar(POM1_VA)
  )
el %>%
  select(
    POD1_VA,
    POD1_LogMAR,
    POW1_VA,
    POW1_LogMAR,
    POM1_VA,
    POM1_LogMAR
  ) %>%
  head(25)
## # A tibble: 25 × 6
##    POD1_VA POD1_LogMAR POW1_VA POW1_LogMAR POM1_VA POM1_LogMAR
##    <chr>         <dbl> <chr>         <dbl> <chr>         <dbl>
##  1 20/20        0      <NA>        NA      20/20         0    
##  2 20/20-1      0.02   <NA>        NA      20/20-2       0.04 
##  3 20/30-1      0.196  20/40        0.301  20/40         0.301
##  4 20/30-1      0.196  20/25        0.0969 20/25-2       0.137
##  5 20/20-1      0.02   20/20        0      20/20         0    
##  6 20/20        0      20/20        0      20/20         0    
##  7 20/25-1      0.117  <NA>        NA      <NA>         NA    
##  8 20/40        0.301  <NA>        NA      <NA>         NA    
##  9 20/15-2     -0.0849 <NA>        NA      <NA>         NA    
## 10 20/30-2      0.216  <NA>        NA      <NA>         NA    
## # ℹ 15 more rows
summary(el$POD1_LogMAR)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -0.24412 -0.10494  0.00000 -0.01009  0.04000  0.69897       21
summary(el$POW1_LogMAR)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -0.22412 -0.12494 -0.02000 -0.04606  0.00000  0.49712      189
summary(el$POM1_LogMAR)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
## -0.24412 -0.12494 -0.04000 -0.05328  0.00000  0.47712      187
#how much are missing value wise
el %>%
  group_by(Procedure) %>%
  summarise(
    Total = n(),
    POD1_Available = sum(!is.na(POD1_LogMAR)),
    POD1_Missing = sum(is.na(POD1_LogMAR)),
    POW1_Available = sum(!is.na(POW1_LogMAR)),
    POW1_Missing = sum(is.na(POW1_LogMAR)),
    POM1_Available = sum(!is.na(POM1_LogMAR)),
    POM1_Missing = sum(is.na(POM1_LogMAR))
  )
## # A tibble: 2 × 8
##   Procedure Total POD1_Available POD1_Missing POW1_Available POW1_Missing
##   <fct>     <int>          <int>        <int>          <int>        <int>
## 1 EVO         282            272           10            189           93
## 2 LASIK       286            275           11            190           96
## # ℹ 2 more variables: POM1_Available <int>, POM1_Missing <int>
el %>%
  group_by(Procedure) %>%
  summarise(
    Total = n(),
    POD1_Percent = mean(!is.na(POD1_LogMAR)) * 100,
    POW1_Percent = mean(!is.na(POW1_LogMAR)) * 100,
    POM1_Percent = mean(!is.na(POM1_LogMAR)) * 100
  )
## # A tibble: 2 × 5
##   Procedure Total POD1_Percent POW1_Percent POM1_Percent
##   <fct>     <int>        <dbl>        <dbl>        <dbl>
## 1 EVO         282         96.5         67.0         63.8
## 2 LASIK       286         96.2         66.4         70.3
#logmar results
va_summary_procedure <- el %>%
  group_by(Procedure) %>%
  summarise(
    n_POD1 = sum(!is.na(POD1_LogMAR)),
    Mean_POD1 = mean(POD1_LogMAR, na.rm = TRUE),
    SD_POD1 = sd(POD1_LogMAR, na.rm = TRUE),
    Median_POD1 = median(POD1_LogMAR, na.rm = TRUE),
    IQR_POD1 = IQR(POD1_LogMAR, na.rm = TRUE),

    n_POW1 = sum(!is.na(POW1_LogMAR)),
    Mean_POW1 = mean(POW1_LogMAR, na.rm = TRUE),
    SD_POW1 = sd(POW1_LogMAR, na.rm = TRUE),
    Median_POW1 = median(POW1_LogMAR, na.rm = TRUE),
    IQR_POW1 = IQR(POW1_LogMAR, na.rm = TRUE),

    n_POM1 = sum(!is.na(POM1_LogMAR)),
    Mean_POM1 = mean(POM1_LogMAR, na.rm = TRUE),
    SD_POM1 = sd(POM1_LogMAR, na.rm = TRUE),
    Median_POM1 = median(POM1_LogMAR, na.rm = TRUE),
    IQR_POM1 = IQR(POM1_LogMAR, na.rm = TRUE)
  )

va_summary_procedure
## # A tibble: 2 × 16
##   Procedure n_POD1 Mean_POD1 SD_POD1 Median_POD1 IQR_POD1 n_POW1 Mean_POW1
##   <fct>      <int>     <dbl>   <dbl>       <dbl>    <dbl>  <int>     <dbl>
## 1 EVO          272    0.0101   0.127           0    0.142    189   -0.0442
## 2 LASIK        275   -0.0300   0.120           0    0.145    190   -0.0479
## # ℹ 8 more variables: SD_POW1 <dbl>, Median_POW1 <dbl>, IQR_POW1 <dbl>,
## #   n_POM1 <int>, Mean_POM1 <dbl>, SD_POM1 <dbl>, Median_POM1 <dbl>,
## #   IQR_POM1 <dbl>
#t test for comparsion
t.test(POD1_LogMAR ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  POD1_LogMAR by Procedure
## t = 3.7975, df = 542.39, p-value = 0.0001627
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.01936493 0.06086746
## sample estimates:
##   mean in group EVO mean in group LASIK 
##          0.01007390         -0.03004229
t.test(POW1_LogMAR ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  POW1_LogMAR by Procedure
## t = 0.32441, df = 367.75, p-value = 0.7458
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -0.01881613  0.02625098
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.04419274         -0.04791016
t.test(POM1_LogMAR ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  POM1_LogMAR by Procedure
## t = 2.2932, df = 347.7, p-value = 0.02243
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.003321369 0.043354766
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.04096843         -0.06430650
el <- el %>%
  mutate(
    POD1_20_20 = case_when(
      is.na(POD1_LogMAR) ~ NA_character_,
      POD1_LogMAR <= 0 ~ "Yes",
      TRUE ~ "No"
    ),
    POW1_20_20 = case_when(
      is.na(POW1_LogMAR) ~ NA_character_,
      POW1_LogMAR <= 0 ~ "Yes",
      TRUE ~ "No"
    ),
    POM1_20_20 = case_when(
      is.na(POM1_LogMAR) ~ NA_character_,
      POM1_LogMAR <= 0 ~ "Yes",
      TRUE ~ "No"
    )
  )
#percentages of model
el %>%
  group_by(Procedure) %>%
  summarise(
    POD1_n = sum(!is.na(POD1_20_20)),
    POD1_20_20_Percent =
      mean(POD1_20_20 == "Yes", na.rm = TRUE) * 100,

    POW1_n = sum(!is.na(POW1_20_20)),
    POW1_20_20_Percent =
      mean(POW1_20_20 == "Yes", na.rm = TRUE) * 100,

    POM1_n = sum(!is.na(POM1_20_20)),
    POM1_20_20_Percent =
      mean(POM1_20_20 == "Yes", na.rm = TRUE) * 100
  )
## # A tibble: 2 × 7
##   Procedure POD1_n POD1_20_20_Percent POW1_n POW1_20_20_Percent POM1_n
##   <fct>      <int>              <dbl>  <int>              <dbl>  <int>
## 1 EVO          272               55.9    189               77.2    180
## 2 LASIK        275               72.4    190               73.7    201
## # ℹ 1 more variable: POM1_20_20_Percent <dbl>
#compare AGE, SE, CYL
t.test(Age ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  Age by Procedure
## t = 0.38863, df = 546.09, p-value = 0.6977
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -1.337412  1.997135
## sample estimates:
##   mean in group EVO mean in group LASIK 
##            34.08820            33.75834
t.test(Pre_SE ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  Pre_SE by Procedure
## t = -19.313, df = 470.18, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -4.393111 -3.581723
## sample estimates:
##   mean in group EVO mean in group LASIK 
##           -6.687276           -2.699859
t.test(abs(Pre_Cyl) ~ Procedure, data = el)
## 
##  Welch Two Sample t-test
## 
## data:  abs(Pre_Cyl) by Procedure
## t = 2.01, df = 562, p-value = 0.04491
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.004141993 0.359634022
## sample estimates:
##   mean in group EVO mean in group LASIK 
##            1.297153            1.115265
#compare age, SE, and CYL
el %>%
 group_by(Procedure) %>%
 summarise(
   Mean_Age = mean(Age, na.rm=TRUE),
   Mean_Pre_SE = mean(Pre_SE, na.rm=TRUE),
   Mean_Pre_Cyl = mean(abs(Pre_Cyl), na.rm=TRUE)
 )
## # A tibble: 2 × 4
##   Procedure Mean_Age Mean_Pre_SE Mean_Pre_Cyl
##   <fct>        <dbl>       <dbl>        <dbl>
## 1 EVO           34.1       -6.69         1.30
## 2 LASIK         33.8       -2.70         1.12
#regression model
summary(
  lm(
    POM1_LogMAR ~ Procedure + Age + Pre_SE + abs(Pre_Cyl),
    data = el
  )
)
## 
## Call:
## lm(formula = POM1_LogMAR ~ Procedure + Age + Pre_SE + abs(Pre_Cyl), 
##     data = el)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22579 -0.07159 -0.00317  0.06595  0.50443 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -0.0522188  0.0233991  -2.232 0.026228 *  
## ProcedureLASIK -0.0041953  0.0119571  -0.351 0.725888    
## Age            -0.0013078  0.0005087  -2.571 0.010531 *  
## Pre_SE         -0.0053861  0.0020227  -2.663 0.008081 ** 
## abs(Pre_Cyl)    0.0193557  0.0049437   3.915 0.000107 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09443 on 375 degrees of freedom
##   (188 observations deleted due to missingness)
## Multiple R-squared:  0.0952, Adjusted R-squared:  0.08554 
## F-statistic: 9.864 on 4 and 375 DF,  p-value: 1.347e-07
#Baseline characteristics differed substantially between treatment groups. EVO ICL patients were significantly more myopic preoperatively than LASIK patients (mean spherical equivalent: -6.69 D vs. -2.70 D, p < 0.001) and also exhibited greater preoperative astigmatism (1.30 D vs. 1.12 D, p = 0.045). Age did not differ significantly between groups (34.1 vs. 33.8 years, p = 0.698).

#In unadjusted analyses, LASIK eyes demonstrated slightly better postoperative month 1 visual acuity than EVO ICL eyes (p = 0.022). However, after adjustment for age, preoperative spherical equivalent, and preoperative cylinder using multivariable linear regression, procedure type was no longer associated with month 1 visual acuity (β = -0.004, p = 0.726). Greater preoperative astigmatism was independently associated with worse postoperative visual acuity (β = 0.019, p < 0.001), while preoperative spherical equivalent also remained a significant predictor (β = -0.005, p = 0.008).
# ============================================================
# Figure 1: Visual Recovery Curve
# ============================================================

library(tidyverse)

va_long <- el %>%
  select(
    Procedure,
    POD1_LogMAR,
    POW1_LogMAR,
    POM1_LogMAR
  ) %>%
  pivot_longer(
    cols = c(POD1_LogMAR, POW1_LogMAR, POM1_LogMAR),
    names_to = "Visit",
    values_to = "LogMAR"
  ) %>%
  mutate(
    Visit = factor(
      Visit,
      levels = c(
        "POD1_LogMAR",
        "POW1_LogMAR",
        "POM1_LogMAR"
      ),
      labels = c(
        "Post-op Day 1",
        "Post-op Week 1",
        "Post-op Month 1"
      )
    )
  )

va_summary <- va_long %>%
  group_by(Procedure, Visit) %>%
  summarise(
    Mean = mean(LogMAR, na.rm = TRUE),
    SE = sd(LogMAR, na.rm = TRUE) /
      sqrt(sum(!is.na(LogMAR))),
    .groups = "drop"
  )

ggplot(
  va_summary,
  aes(
    x = Visit,
    y = Mean,
    color = Procedure,
    group = Procedure
  )
) +
  geom_line(linewidth = 1.4) +
  geom_point(size = 4) +
  geom_errorbar(
    aes(
      ymin = Mean - SE,
      ymax = Mean + SE
    ),
    width = 0.1,
    linewidth = 0.8
  ) +
  scale_color_manual(
    values = c(
      "EVO" = "#1F4E79",
      "LASIK" = "#C00000"
    )
  ) +
  labs(
    title = "Visual Recovery Following EVO ICL and LASIK",
    subtitle = "Mean LogMAR Visual Acuity ± Standard Error",
    x = "",
    y = "LogMAR Visual Acuity",
    color = "Procedure"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

ggplot(
  el,
  aes(
    x = Procedure,
    y = POM1_LogMAR,
    fill = Procedure
  )
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.35,
    size = 1.8
  ) +
  scale_fill_manual(
    values = c(
      "EVO" = "#1F4E79",
      "LASIK" = "#C00000"
    )
  ) +
  labs(
    title = "Month 1 Visual Acuity by Procedure",
    subtitle = "Distribution of LogMAR Outcomes",
    x = "",
    y = "Month 1 LogMAR"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none"
  )

ggplot(
  el,
  aes(
    x = Pre_SE,
    fill = Procedure
  )
) +
  geom_density(
    alpha = 0.4,
    linewidth = 1
  ) +
  scale_fill_manual(
    values = c(
      "EVO" = "#1F4E79",
      "LASIK" = "#C00000"
    )
  ) +
  labs(
    title = "Distribution of Preoperative Spherical Equivalent",
    subtitle = "EVO Patients Were Significantly More Myopic",
    x = "Preoperative Spherical Equivalent (D)",
    y = "Density"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

va20_summary <- el %>%
  group_by(Procedure) %>%
  summarise(
    POD1 = mean(POD1_LogMAR <= 0, na.rm = TRUE) * 100,
    POW1 = mean(POW1_LogMAR <= 0, na.rm = TRUE) * 100,
    POM1 = mean(POM1_LogMAR <= 0, na.rm = TRUE) * 100
  ) %>%
  pivot_longer(
    cols = POD1:POM1,
    names_to = "Visit",
    values_to = "Percent"
  )


ggplot(
  va20_summary,
  aes(
    x = Visit,
    y = Percent,
    fill = Procedure
  )
) +
  geom_col(
    position = position_dodge(0.7),
    width = 0.65
  ) +
  scale_fill_manual(
    values = c(
      "EVO" = "#1F4E79",
      "LASIK" = "#C00000"
    )
  ) +
  labs(
    title = "Eyes Achieving 20/20 or Better Vision",
    subtitle = "By Postoperative Follow-up Visit",
    x = "",
    y = "Percent (%)"
  ) +
  theme_classic(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

el <- el %>%
  mutate(
    Myopia_Group = case_when(
      Pre_SE <= -6.00 ~ "High Myopia (≤ -6 D)",
      Pre_SE > -6.00  ~ "Low/Moderate Myopia (> -6 D)"
    )
  )

el$Myopia_Group <- factor(
  el$Myopia_Group,
  levels = c(
    "Low/Moderate Myopia (> -6 D)",
    "High Myopia (≤ -6 D)"
  )
)
low_myopia <- el %>%
  filter(Myopia_Group == "Low/Moderate Myopia (> -6 D)")

t.test(POD1_LogMAR ~ Procedure, data = low_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POD1_LogMAR by Procedure
## t = 2.061, df = 258.48, p-value = 0.0403
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.001127145 0.049466515
## sample estimates:
##   mean in group EVO mean in group LASIK 
##        -0.004206706        -0.029503536
t.test(POW1_LogMAR ~ Procedure, data = low_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POW1_LogMAR by Procedure
## t = -1.2687, df = 203.67, p-value = 0.206
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -0.043531160  0.009443035
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.06279981         -0.04575575
t.test(POM1_LogMAR ~ Procedure, data = low_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POM1_LogMAR by Procedure
## t = 0.93981, df = 160.05, p-value = 0.3487
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -0.01204484  0.03391710
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.05246145         -0.06339758
high_myopia <- el %>%
  filter(Myopia_Group == "High Myopia (≤ -6 D)")

t.test(POD1_LogMAR ~ Procedure, data = high_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POD1_LogMAR by Procedure
## t = 2.1675, df = 18.467, p-value = 0.04349
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.002008555 0.121682326
## sample estimates:
##   mean in group EVO mean in group LASIK 
##          0.02175920         -0.04008624
t.test(POW1_LogMAR ~ Procedure, data = high_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POW1_LogMAR by Procedure
## t = 2.6382, df = 4.736, p-value = 0.04863
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  0.0008829715 0.1941727096
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.03009565         -0.12762349
t.test(POM1_LogMAR ~ Procedure, data = high_myopia)
## 
##  Welch Two Sample t-test
## 
## data:  POM1_LogMAR by Procedure
## t = 1.7116, df = 25.914, p-value = 0.09891
## alternative hypothesis: true difference in means between group EVO and group LASIK is not equal to 0
## 95 percent confidence interval:
##  -0.008678006  0.094976028
## sample estimates:
##   mean in group EVO mean in group LASIK 
##         -0.03166684         -0.07481585
interaction_model <- lm(
  POM1_LogMAR ~ Procedure * Myopia_Group + Age + abs(Pre_Cyl),
  data = el
)

summary(interaction_model)
## 
## Call:
## lm(formula = POM1_LogMAR ~ Procedure * Myopia_Group + Age + abs(Pre_Cyl), 
##     data = el)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.23590 -0.06963  0.00147  0.06560  0.52149 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                     -0.0295885  0.0214436  -1.380
## ProcedureLASIK                                  -0.0119376  0.0125749  -0.949
## Myopia_GroupHigh Myopia (≤ -6 D)                 0.0166812  0.0143504   1.162
## Age                                             -0.0012802  0.0005132  -2.495
## abs(Pre_Cyl)                                     0.0210940  0.0049946   4.223
## ProcedureLASIK:Myopia_GroupHigh Myopia (≤ -6 D) -0.0397986  0.0286500  -1.389
##                                                 Pr(>|t|)    
## (Intercept)                                        0.168    
## ProcedureLASIK                                     0.343    
## Myopia_GroupHigh Myopia (≤ -6 D)                   0.246    
## Age                                                0.013 *  
## abs(Pre_Cyl)                                    3.03e-05 ***
## ProcedureLASIK:Myopia_GroupHigh Myopia (≤ -6 D)    0.166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09517 on 374 degrees of freedom
##   (188 observations deleted due to missingness)
## Multiple R-squared:  0.08356,    Adjusted R-squared:  0.07131 
## F-statistic:  6.82 on 5 and 374 DF,  p-value: 4.19e-06
ggplot(
  el,
  aes(
    x = Procedure,
    y = POM1_LogMAR,
    fill = Procedure
  )
) +
  geom_boxplot(alpha = 0.8) +
  geom_jitter(width = 0.15, alpha = 0.3) +
  facet_wrap(~Myopia_Group) +
  scale_fill_manual(
    values = c(
      "EVO" = "#1F4E79",
      "LASIK" = "#C00000"
    )
  ) +
  labs(
    title = "Month 1 Visual Acuity by Procedure and Myopia Group",
    x = "",
    y = "Month 1 LogMAR"
  ) +
  theme_classic(base_size = 15) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "none",
    strip.text = element_text(face = "bold")
  )

ggplot(el,
       aes(Procedure,
           POM1_LogMAR,
           fill = Procedure)) +
  geom_violin(trim = FALSE,
              alpha = 0.7) +
  geom_boxplot(width = 0.12,
               fill = "white") +
  facet_wrap(~Myopia_Group)