Overview

This analysis examines pulmonary function test (PFT) data across three time periods: - PRE: Pre-pandemic - PAN: During pandemic - POST: Post-pandemic

We will perform: 1. One-way ANOVA for continuous variables (FVC%, FVC Absolute, FEV1%, FEV1 Absolute, FEV1/FVC) 2. Chi-squared test for categorical variable (Sawtooth Pattern presence)

Load Required Libraries

library(tidyverse)
library(car)       
library(ggplot2)
library(gridExtra)  
library(knitr)      
library(broom)

Data Entry

FVC% Data

fvc_percent_pre <- c(1.125827815, 0.8934579439, 1.001751313, 1.138028169, 0.9933628319, 
          1.023622047, 1.015250545, 1.041436464, 1.202087287, 1.210268949, 
          1.130319149, 1.001312336, 1.15530303, 1.30945122, 1.080615942, 
          1.116972477, 1.136075949, 0.9304511278, 1.154022989, 0.9736842105, 
          0.8197969543, 1.063829787, 1.07745098, 1.178149606, 1.169340463, 
          1.186725664, 1.077108434, 1.23853211, 1.090196078, 1.248511905, 
          1.016055046, 1.041845494, 1.14, 1.167481663, 1.082236842, 
          0.9478527607, 0.9860426929, 1.2079566, 1.340753425, 1.017751479, 
          0.9667832168, 1.171428571, 1.173486088, 0.9443493151, 0.885286783, 
          1.05511811, 0.8910550459, 0.6631205674, 1.075875486, 1.112944162, 
          1.120171674)

fvc_percent_pan <- c(1.073563218, 1.050245098, 0.8611987382, 1.343010753, 1.018446602, 
          0.957337884, 1, 1.169054441, 0.8679245283, 1.042009885, 
          1.085280374, 0.9787018256, 1.036666667, 1.061611374, 0.9673659674, 
          1.098669623, 1.109289617, 0.9569160998, 0.8655172414, 1.102134146, 
          1.093478261, 0.9368932039, 0.8425287356, 0.8228643216, 0.4581632653, 
          1.040772532, 0.8959435626, 0.9973544974, 1.00776699, 0.9832935561, 
          1.073253833, 1.15, 1.03490566, 0.9513556619, 1.100543478, 
          0.8669656203, 1.003773585, 1.17414248, 1.081613508, 0.9090909091, 
          1.019699812, 1.067028986, 1.138086643, 1.457236842, 1.055405405, 
          1.037694013, 1.056872038, 1.050713154, 0.9340490798, 1.0225, 
          1.216066482)

fvc_percent_post <- c(1.032062392, 1.029484029, 0.9576659039, 1.214650767, 0.9807302231, 
           1.213207547, 0.9259259259, 1.03436019, 1.254452926, 1.081771721, 
           0.8295165394, 1.083333333, 1.031806616, 0.6995412844, 0.9435483871, 
           1.443064182, 1.191295547, 1.066157761, 0.9398034398, 1.107583774, 
           0.883248731, 1.224683544, 1.073015873, 1.106617647, 1.132379249, 
           0.9255150555, 0.9561403509, 1.239669421, 1.063451777, 1.094202899, 
           1.010141988, 0.9141791045, 1.086288416, 1, 1.069852941, 
           1.164691943, 0.9885496183, 0.8510802469, 0.9197080292, 1.23875, 
           1.081743869, 1.014251781, 0.9775725594, 1.075, 1.011930586, 
           1.092819615, 0.83847981, 1.190239044, 1.095823096, 0.7485928705, 
           1.15049505, 0.9677754678, 1.00525394, 1.119089317, 0.7622820919, 
           0.9567901235, 0.9557522124, 1.061165049, 0.9966159052, 0.8880070547, 
           1.003623188, 0.8348765432, 0.8955223881, 1.152061856, 1.067681895, 
           1.064367816, 0.9527145359, 0.8631840796, 1.158119658, 1.053333333, 
           1.088832487, 1.096226415, 0.8657894737, 0.7826086957, 1.203448276, 
           0.937037037, 1.03562341, 1.079422383, 1.019417476, 1.453125, 
           1.195911414, 0.9459102902, 1.190709046, 1.031301483, 0.7691588785, 
           0.8957783641, 0.7240527183, 0.92228739, 0.9637681159, 0.9796954315, 
           0.9565789474, 1.03785489, 0.9105263158)

fvc_percent <- data.frame(
  Period = c(rep("PRE", length(fvc_percent_pre)), 
             rep("PAN", length(fvc_percent_pan)), 
             rep("POST", length(fvc_percent_post))),
  Value = c(fvc_percent_pre, fvc_percent_pan, fvc_percent_post)
) %>%
  mutate(Period = factor(Period, levels = c("PRE", "PAN", "POST")))

FVC Absolute Average Data

fvc_absolute_pre <- c(3.4, 4.78, 5.72, 4.04, 4.49, 3.9, 4.66, 3.77, 6.335, 4.95, 
          4.25, 3.815, 4.575, 4.295, 5.965, 4.87, 5.385, 4.95, 5.02, 3.7, 
          3.23, 4, 5.495, 5.985, 6.56, 6.705, 4.47, 5.4, 5.56, 4.195, 
          4.43, 4.855, 3.99, 4.775, 6.58, 3.09, 6.005, 6.68, 7.83, 3.44, 
          5.53, 4.1, 7.17, 5.515, 3.55, 4.02, 3.885, 3.74, 5.53, 4.385, 
          5.22)

fvc_absolute_pan <- c(4.67, 4.285, 5.46, 6.245, 5.245, 5.61, 5.3, 4.08, 4.6, 6.325, 
          4.645, 4.825, 4.665, 4.48, 4.15, 4.955, 4.06, 4.22, 3.765, 3.615, 
          7.55, 4.825, 3.665, 3.275, 2.245, 4.85, 5.08, 3.77, 5.19, 4.12, 
          6.3, 6.095, 5.49, 5.965, 6.075, 5.8, 5.32, 4.45, 5.765, 2.9, 
          5.435, 5.89, 6.305, 4.43, 3.905, 4.68, 4.46, 6.63, 3.045, 4.09, 
          4.39)

fvc_absolute_post <- c(5.955, 4.19, 4.185, 7.13, 4.835, 6.43, 6, 4.365, 4.93, 6.35, 
           3.26, 3.965, 4.055, 3.05, 4.68, 6.97, 5.885, 4.19, 3.825, 6.28, 
           5.22, 5.805, 3.38, 4.515, 6.33, 5.84, 5.995, 4.5, 4.19, 6.04, 
           4.98, 3.675, 4.595, 3.41, 4.365, 4.915, 3.885, 5.515, 5.04, 4.955, 
           3.97, 4.27, 3.705, 6.45, 4.665, 6.24, 3.53, 5.975, 4.46, 3.99, 
           5.81, 4.655, 5.74, 6.39, 4.81, 5.425, 4.32, 5.465, 5.89, 5.035, 
           6.925, 5.41, 4.8, 6.705, 6.31, 4.63, 5.44, 3.47, 4.065, 4.74, 
           4.29, 5.81, 3.29, 4.32, 5.235, 5.06, 4.07, 5.98, 5.25, 6.51, 
           7.02, 3.585, 4.87, 6.26, 4.115, 3.395, 4.395, 3.145, 5.32, 
           5.79, 3.635, 6.58, 3.46)

fvc_absolute <- data.frame(
  Period = c(rep("PRE", length(fvc_absolute_pre)), 
             rep("PAN", length(fvc_absolute_pan)), 
             rep("POST", length(fvc_absolute_post))),
  Value = c(fvc_absolute_pre, fvc_absolute_pan, fvc_absolute_post)
) %>%
  mutate(Period = factor(Period, levels = c("PRE", "PAN", "POST")))

FEV1% Data

fev1_percent_pre <- c(0.9759259259, 0.8509933775, 0.9802904564, 0.9870967742, 1.037275064, 
          1.055891239, 0.9987531172, 0.9921875, 1.227171492, 1.143258427, 
          1.096969697, 1.027108434, 1.130813953, 1.227739726, 0.8683083512, 
          1.080687831, 1.153186275, 0.7605321508, 0.9498680739, 0.9234234234, 
          0.8226744186, 1.036363636, 1.104597701, 1.101149425, 1.180293501, 
          1.107740586, 1.02892562, 1.169312169, 1.049425287, 1.258333333, 
          1.005291005, 1.002487562, 1.211290323, 1.171875, 0.9774951076, 
          0.9597902098, 0.8857421875, 1.128479657, 1.076219512, 0.9598662207, 
          0.8226141079, 1.119354839, 1.153996101, 1.071283096, 0.9419263456, 
          1.063253012, 0.8902116402, 0.8221757322, 1.02283105, 0.8648255814, 
          1.012437811)

fev1_percent_pan <- c(1.116094987, 0.9788732394, 0.8371212121, 0.8660049628, 0.997716895, 
          0.9716024341, 1.002222222, 1.035483871, 0.8155555556, 0.7941176471, 
          1.089333333, 0.9265402844, 0.9654731458, 0.8160762943, 0.9, 
          1.034615385, 0.9068322981, 1.058139535, 0.8562005277, 1.05862069, 
          0.8886956522, 0.7694063927, 0.7189973615, 0.8614285714, 0.3627684964, 
          1.099502488, 0.9125, 0.963963964, 0.9285714286, 1.021798365, 
          0.6232323232, 1.112222222, 1, 0.9058935361, 1.003211991, 
          0.8422939068, 1.006666667, 1.073573574, 1.064301552, 0.9642857143, 
          0.9723451327, 1.026766595, 1.048179872, 1.383333333, 1.043209877, 
          1.016666667, 0.9972752044, 0.8816287879, 0.9775862069, 0.9114285714, 
          1.174050633)

fev1_percent_post <- c(0.9630390144, 0.9225352113, 0.8422459893, 1.201010101, 0.9952606635, 
           1.213333333, 1.046125461, 1.054347826, 1.300872093, 1.062626263, 
           0.8880813953, 1.096273292, 1.049418605, 0.7552910053, 0.9549763033, 
           1.498783455, 1.04009434, 1.11627907, 0.9957746479, 1.114583333, 
           0.9436619718, 1.243872549, 1.041666667, 1.004225352, 1.053571429, 
           0.7755681818, 0.8108365019, 1.059375, 1.064841499, 0.9421841542, 
           0.990521327, 0.9670487106, 0.900273224, 0.8183333333, 1.009859155, 
           1.108991826, 0.7456395349, 0.8118081181, 0.9849462366, 1.144475921, 
           1.040372671, 1.032697548, 1.016516517, 1.015810277, 0.9409638554, 
           1.102697095, 0.89373297, 1.096736597, 1.177464789, 0.7212389381, 
           1.024249423, 0.9939759036, 1.031120332, 1.064315353, 0.7547348485, 
           0.9, 1.014248705, 0.997716895, 0.9839034205, 0.696875, 0.9426086957, 
           0.8256457565, 0.9303097345, 1.032064128, 1.020120724, 0.9406332454, 
           0.9139004149, 0.8567335244, 1.16828479, 1.034526854, 1.068313953, 
           1.104444444, 0.7717717718, 0.8286937901, 1.211081794, 0.9287280702, 
           1.066860465, 0.9646680942, 1.020547945, 1.3671875, 1.135353535, 
           0.9624624625, 1.110169492, 1.019607843, 0.7953539823, 0.9294294294, 
           0.8352941176, 0.9783333333, 1.027837259, 1.044265594, 0.966966967, 
           1.035984848, 0.987987988)

fev1_percent <- data.frame(
  Period = c(rep("PRE", length(fev1_percent_pre)), 
             rep("PAN", length(fev1_percent_pan)), 
             rep("POST", length(fev1_percent_post))),
  Value = c(fev1_percent_pre, fev1_percent_pan, fev1_percent_post)
) %>%
  mutate(Period = factor(Period, levels = c("PRE", "PAN", "POST")))

FEV1 Absolute Average Data

fev1_absolute_pre <- c(2.635, 3.855, 4.725, 3.06, 4.035, 3.495, 4.005, 3.175, 5.51, 4.07, 
          3.62, 3.41, 3.89, 3.585, 4.055, 4.085, 4.705, 3.43, 3.6, 3.075, 
          2.83, 3.42, 4.805, 4.79, 5.63, 5.295, 3.735, 4.42, 4.565, 3.775, 
          3.8, 4.03, 3.755, 4.125, 4.995, 2.745, 4.535, 5.27, 5.295, 2.87, 
          3.965, 3.47, 5.92, 5.26, 3.325, 3.53, 3.365, 3.93, 4.48, 2.975, 
          4.07)

fev1_absolute_pan <- c(4.23, 3.475, 4.42, 3.49, 4.37, 4.79, 4.51, 3.21, 3.67, 4.05, 
          4.085, 3.91, 3.775, 2.995, 3.375, 4.035, 2.92, 4.095, 3.245, 3.07, 
          5.11, 3.37, 2.725, 3.015, 1.52, 4.42, 4.38, 3.21, 4.095, 3.75, 
          3.085, 5.005, 4.5, 4.765, 4.685, 4.7, 4.53, 3.575, 4.8, 2.565, 
          4.395, 4.795, 4.895, 3.735, 3.38, 3.965, 3.66, 4.655, 2.835, 3.19, 
          3.71)

fev1_absolute_post <- c(4.69, 3.275, 3.15, 5.945, 4.2, 5.46, 5.67, 3.88, 4.475, 5.26, 
           3.055, 3.53, 3.61, 2.855, 4.03, 6.16, 4.41, 3.84, 3.535, 5.35, 
           4.69, 5.075, 2.875, 3.565, 5.015, 4.095, 4.265, 3.39, 3.695, 4.4, 
           4.18, 3.375, 3.295, 2.455, 3.585, 4.07, 2.565, 4.4, 4.58, 4.04, 
           3.35, 3.79, 3.385, 5.14, 3.905, 5.315, 3.28, 4.705, 4.18, 3.26, 
           4.435, 4.125, 4.97, 5.13, 3.985, 4.32, 3.915, 4.37, 4.89, 3.345, 
           5.42, 4.475, 4.205, 5.15, 5.07, 3.565, 4.405, 2.99, 3.61, 4.045, 
           3.675, 4.97, 2.57, 3.87, 4.59, 4.235, 3.67, 4.505, 4.47, 5.25, 
           5.62, 3.205, 3.93, 5.2, 3.595, 3.095, 4.26, 2.935, 4.8, 5.19, 
           3.22, 5.47, 3.29)

fev1_absolute <- data.frame(
  Period = c(rep("PRE", length(fev1_absolute_pre)), 
             rep("PAN", length(fev1_absolute_pan)), 
             rep("POST", length(fev1_absolute_post))),
  Value = c(fev1_absolute_pre, fev1_absolute_pan, fev1_absolute_post)
) %>%
  mutate(Period = factor(Period, levels = c("PRE", "PAN", "POST")))

FEV1/FVC Ratio Data

fev1_fvc_ratio_pre <- c(87, 81, 83, 76, 90, 90, 87, 84, 87, 82, 85, 89.5, 85, 83, 68, 84, 
          87.5, 69.5, 81, 84, 87.5, 86, 87.5, 80, 86, 82, 83.5, 83, 82, 90, 
          86, 82, 94, 86.5, 76, 88.5, 75.5, 79, 72, 83, 71.5, 85, 82.5, 95.5, 
          94.5, 88, 86.5, 68, 81, 69.5, 80)

fev1_fvc_ratio_pan <- c(91, 89, 83, 86, 83.5, 85.5, 78, 78.5, 81, 64, 88, 81, 81, 67, 81.5, 
          81.5, 72, 97, 86.5, 85, 67.5, 69.5, 74.5, 92, 72.5, 91, 86, 85, 79, 
          91, 49, 82, 82, 79.5, 77, 81, 85.5, 80.5, 83, 88.5, 81, 81.5, 78, 
          84, 87, 84.5, 82.5, 70.5, 93, 78, 84.5)

fev1_fvc_ratio_post <- c(79, 78, 75.5, 83.5, 87, 85, 95, 88.5, 90.5, 83, 94, 94.5, 88.5, 93.5, 
           86, 88, 75, 92, 92.5, 85, 89, 87, 85, 79, 79, 70, 71, 75, 87.5, 73, 
           84, 91.5, 71.5, 72, 82, 83, 66, 80, 91, 81.5, 84, 89, 91.5, 79.5, 
           84, 85.5, 92.5, 79, 94, 81.5, 76.5, 88.5, 87, 80, 83, 79.5, 90.5, 
           80, 83, 66.5, 78.5, 82.5, 87.5, 77, 80.5, 77, 81, 86, 88.5, 85.5, 
           85.5, 86, 78.5, 90, 88, 84, 90, 75.5, 85, 81, 80, 90, 81, 83, 87, 
           91.5, 97, 93.5, 90, 90, 88.5, 83, 95)

fev1_fvc_ratio <- data.frame(
  Period = c(rep("PRE", length(fev1_fvc_ratio_pre)), 
             rep("PAN", length(fev1_fvc_ratio_pan)), 
             rep("POST", length(fev1_fvc_ratio_post))),
  Value = c(fev1_fvc_ratio_pre, fev1_fvc_ratio_pan, fev1_fvc_ratio_post)
) %>%
  mutate(Period = factor(Period, levels = c("PRE", "PAN", "POST")))

Sawtooth Pattern Data

sawtooth_data <- data.frame(
  Period = c("PRE", "PAN", "POST"),
  Present = c(12, 30, 41),
  Total = c(51, 51, 93),
  Absent = c(51-12, 51-30, 93-41)
)

sawtooth_data$Percent <- (sawtooth_data$Present / sawtooth_data$Total) * 100

Descriptive Statistics

Summary Statistics Function

get_summary_stats <- function(data, variable_name) {
  
  summary_stats <- data %>%
    group_by(Period) %>%
    summarise(
      N = n(),
      Mean = mean(Value),
      SD = sd(Value),
      Median = median(Value),
      Min = min(Value),
      Max = max(Value),
      .groups = 'drop'
    )
  
  return(summary_stats)
}

FVC% Summary

fvc_pct_summary <- get_summary_stats(fvc_percent, "FVC%")
kable(fvc_pct_summary, digits = 3, caption = "Summary Statistics for FVC%")
Summary Statistics for FVC%
Period N Mean SD Median Min Max
PRE 51 1.074 0.125 1.081 0.663 1.341
PAN 51 1.022 0.142 1.037 0.458 1.457
POST 93 1.023 0.139 1.029 0.700 1.453

FVC Absolute Summary

fvc_abs_summary <- get_summary_stats(fvc_absolute, "FVC Absolute")
kable(fvc_abs_summary, digits = 3, caption = "Summary Statistics for FVC Absolute Average")
Summary Statistics for FVC Absolute Average
Period N Mean SD Median Min Max
PRE 51 4.878 1.097 4.775 3.090 7.83
PAN 51 4.847 1.053 4.680 2.245 7.55
POST 93 4.972 1.064 4.870 3.050 7.13

FEV1% Summary

fev1_pct_summary <- get_summary_stats(fev1_percent, "FEV1%")
kable(fev1_pct_summary, digits = 3, caption = "Summary Statistics for FEV1%")
Summary Statistics for FEV1%
Period N Mean SD Median Min Max
PRE 51 1.029 0.119 1.029 0.761 1.258
PAN 51 0.953 0.149 0.972 0.363 1.383
POST 93 0.998 0.137 1.010 0.697 1.499

FEV1 Absolute Summary

fev1_abs_summary <- get_summary_stats(fev1_absolute, "FEV1 Absolute")
kable(fev1_abs_summary, digits = 3, caption = "Summary Statistics for FEV1 Absolute Average")
Summary Statistics for FEV1 Absolute Average
Period N Mean SD Median Min Max
PRE 51 4.020 0.805 3.930 2.635 5.92
PAN 51 3.858 0.754 3.910 1.520 5.11
POST 93 4.154 0.834 4.125 2.455 6.16

FEV1/FVC Ratio Summary

ratio_summary <- get_summary_stats(fev1_fvc_ratio, "FEV1/FVC")
kable(ratio_summary, digits = 3, caption = "Summary Statistics for FEV1/FVC Ratio")
Summary Statistics for FEV1/FVC Ratio
Period N Mean SD Median Min Max
PRE 51 83.049 6.483 84 68 95.5
PAN 51 81.186 8.264 82 49 97.0
POST 93 84.129 6.676 85 66 97.0

Sawtooth Pattern Summary

kable(sawtooth_data, digits = 2, caption = "Sawtooth Pattern Presence by Period")
Sawtooth Pattern Presence by Period
Period Present Total Absent Percent
PRE 12 51 39 23.53
PAN 30 51 21 58.82
POST 41 93 52 44.09

One-Way ANOVA Tests

ANOVA Function

perform_anova <- function(data, test_name) {
  anova_result <- aov(Value ~ Period, data = data)
  levene_test <- leveneTest(Value ~ Period, data = data)
  
 
  tukey_result <- TukeyHSD(anova_result)
  
  
  results <- list(
    test_name = test_name,
    data = data,
    anova = anova_result,
    anova_summary = summary(anova_result),
    levene = levene_test,
    tukey = tukey_result
  )
  
  return(results)
}

FVC% ANOVA

fvc_pct_anova <- perform_anova(fvc_percent, "FVC%")

cat("\n=== FVC% ANOVA Results ===\n")
## 
## === FVC% ANOVA Results ===
print(fvc_pct_anova$anova_summary)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## Period        2  0.101 0.05048   2.707 0.0693 .
## Residuals   192  3.581 0.01865                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
## 
## --- Levene's Test for Homogeneity of Variance ---
print(fvc_pct_anova$levene)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.4401 0.6446
##       192
cat("\n--- Tukey HSD Post-Hoc Test ---\n")
## 
## --- Tukey HSD Post-Hoc Test ---
print(fvc_pct_anova$tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Value ~ Period, data = data)
## 
## $Period
##                   diff         lwr         upr     p adj
## PAN-PRE  -0.0522654627 -0.11614857 0.011617650 0.1323690
## POST-PRE -0.0515032364 -0.10771288 0.004706408 0.0800368
## POST-PAN  0.0007622263 -0.05544742 0.056971871 0.9994345

FVC Absolute ANOVA

fvc_abs_anova <- perform_anova(fvc_absolute, "FVC Absolute")

cat("\n=== FVC Absolute ANOVA Results ===\n")
## 
## === FVC Absolute ANOVA Results ===
print(fvc_abs_anova$anova_summary)
##              Df Sum Sq Mean Sq F value Pr(>F)
## Period        2    0.6  0.3016   0.263  0.769
## Residuals   192  219.9  1.1452
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
## 
## --- Levene's Test for Homogeneity of Variance ---
print(fvc_abs_anova$levene)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.1934 0.8243
##       192
cat("\n--- Tukey HSD Post-Hoc Test ---\n")
## 
## --- Tukey HSD Post-Hoc Test ---
print(fvc_abs_anova$tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Value ~ Period, data = data)
## 
## $Period
##                 diff        lwr       upr     p adj
## PAN-PRE  -0.03137255 -0.5319338 0.4691887 0.9879912
## POST-PRE  0.09332385 -0.3471113 0.5337590 0.8711591
## POST-PAN  0.12469639 -0.3157388 0.5651316 0.7819119

FEV1% ANOVA

fev1_pct_anova <- perform_anova(fev1_percent, "FEV1%")

cat("\n=== FEV1% ANOVA Results ===\n")
## 
## === FEV1% ANOVA Results ===
print(fev1_pct_anova$anova_summary)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## Period        2  0.152 0.07623   4.113 0.0178 *
## Residuals   192  3.559 0.01853                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
## 
## --- Levene's Test for Homogeneity of Variance ---
print(fev1_pct_anova$levene)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.0883 0.9156
##       192
cat("\n--- Tukey HSD Post-Hoc Test ---\n")
## 
## --- Tukey HSD Post-Hoc Test ---
print(fev1_pct_anova$tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Value ~ Period, data = data)
## 
## $Period
##                 diff         lwr         upr     p adj
## PAN-PRE  -0.07661438 -0.14029502 -0.01293374 0.0137172
## POST-PRE -0.03073542 -0.08676692  0.02529607 0.3993469
## POST-PAN  0.04587896 -0.01015254  0.10191045 0.1319423

FEV1 Absolute ANOVA

fev1_abs_anova <- perform_anova(fev1_absolute, "FEV1 Absolute")

cat("\n=== FEV1 Absolute ANOVA Results ===\n")
## 
## === FEV1 Absolute ANOVA Results ===
print(fev1_abs_anova$anova_summary)
##              Df Sum Sq Mean Sq F value Pr(>F)
## Period        2   2.93   1.464   2.253  0.108
## Residuals   192 124.80   0.650
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
## 
## --- Levene's Test for Homogeneity of Variance ---
print(fev1_abs_anova$levene)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.3855 0.6807
##       192
cat("\n--- Tukey HSD Post-Hoc Test ---\n")
## 
## --- Tukey HSD Post-Hoc Test ---
print(fev1_abs_anova$tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Value ~ Period, data = data)
## 
## $Period
##                diff         lwr       upr     p adj
## PAN-PRE  -0.1618627 -0.53898240 0.2152569 0.5690565
## POST-PRE  0.1346300 -0.19719107 0.4664510 0.6040750
## POST-PAN  0.2964927 -0.03532833 0.6283138 0.0903550

FEV1/FVC Ratio ANOVA

ratio_anova <- perform_anova(fev1_fvc_ratio, "FEV1/FVC Ratio")

cat("\n=== FEV1/FVC Ratio ANOVA Results ===\n")
## 
## === FEV1/FVC Ratio ANOVA Results ===
print(ratio_anova$anova_summary)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## Period        2    285  142.64   2.848 0.0604 .
## Residuals   192   9617   50.09                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Levene's Test for Homogeneity of Variance ---\n")
## 
## --- Levene's Test for Homogeneity of Variance ---
print(ratio_anova$levene)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2  0.4863 0.6156
##       192
cat("\n--- Tukey HSD Post-Hoc Test ---\n")
## 
## --- Tukey HSD Post-Hoc Test ---
print(ratio_anova$tukey)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Value ~ Period, data = data)
## 
## $Period
##               diff         lwr      upr     p adj
## PAN-PRE  -1.862745 -5.17310481 1.447615 0.3807819
## POST-PRE  1.080013 -1.83271556 3.992741 0.6561971
## POST-PAN  2.942758  0.03002954 5.855486 0.0470644

Chi-Squared Test for Sawtooth Pattern

sawtooth_matrix <- matrix(c(sawtooth_data$Present, sawtooth_data$Absent), 
                          nrow = 3, byrow = FALSE,
                          dimnames = list(Period = c("PRE", "PAN", "POST"),
                                        Pattern = c("Present", "Absent")))

cat("\n=== Contingency Table ===\n")
## 
## === Contingency Table ===
print(sawtooth_matrix)
##       Pattern
## Period Present Absent
##   PRE       12     39
##   PAN       30     21
##   POST      41     52
# Chi-squared test
chi_sq_result <- chisq.test(sawtooth_matrix)

cat("\n=== Chi-Squared Test Results ===\n")
## 
## === Chi-Squared Test Results ===
print(chi_sq_result)
## 
##  Pearson's Chi-squared test
## 
## data:  sawtooth_matrix
## X-squared = 13.162, df = 2, p-value = 0.001387
# Expected frequencies
cat("\n--- Expected Frequencies ---\n")
## 
## --- Expected Frequencies ---
print(chi_sq_result$expected)
##       Pattern
## Period  Present   Absent
##   PRE  21.70769 29.29231
##   PAN  21.70769 29.29231
##   POST 39.58462 53.41538
# Standardized residuals
cat("\n--- Standardized Residuals ---\n")
## 
## --- Standardized Residuals ---
print(chi_sq_result$stdres)
##       Pattern
## Period    Present     Absent
##   PRE  -3.1992894  3.1992894
##   PAN   2.7328320 -2.7328320
##   POST  0.4104278 -0.4104278

Visualizations

FVC% Visualization

p1 <- ggplot(fvc_pct_anova$data, aes(x = Period, y = Value, fill = Period)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "FVC% by Period", y = "FVC%", x = "Period") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))

p2 <- ggplot(fvc_pct_summary, aes(x = Period, y = Mean, group = 1)) +
  geom_line(color = "#2E86AB", size = 1.2) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Mean FVC% with Standard Deviation", y = "Mean FVC%", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

FVC Absolute Visualization

p1 <- ggplot(fvc_abs_anova$data, aes(x = Period, y = Value, fill = Period)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "FVC Absolute by Period", y = "FVC (L)", x = "Period") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))

p2 <- ggplot(fvc_abs_summary, aes(x = Period, y = Mean, group = 1)) +
  geom_line(color = "#2E86AB", size = 1.2) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Mean FVC Absolute with Standard Deviation", y = "Mean FVC (L)", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

FEV1% Visualization

p1 <- ggplot(fev1_pct_anova$data, aes(x = Period, y = Value, fill = Period)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "FEV1% by Period", y = "FEV1%", x = "Period") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))

p2 <- ggplot(fev1_pct_summary, aes(x = Period, y = Mean, group = 1)) +
  geom_line(color = "#2E86AB", size = 1.2) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Mean FEV1% with Standard Deviation", y = "Mean FEV1%", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

FEV1 Absolute Visualization

p1 <- ggplot(fev1_abs_anova$data, aes(x = Period, y = Value, fill = Period)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "FEV1 Absolute by Period", y = "FEV1 (L)", x = "Period") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))

p2 <- ggplot(fev1_abs_summary, aes(x = Period, y = Mean, group = 1)) +
  geom_line(color = "#2E86AB", size = 1.2) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Mean FEV1 Absolute with Standard Deviation", y = "Mean FEV1 (L)", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

FEV1/FVC Ratio Visualization

p1 <- ggplot(ratio_anova$data, aes(x = Period, y = Value, fill = Period)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "FEV1/FVC Ratio by Period", y = "FEV1/FVC (%)", x = "Period") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))

p2 <- ggplot(ratio_summary, aes(x = Period, y = Mean, group = 1)) +
  geom_line(color = "#2E86AB", size = 1.2) +
  geom_point(size = 3, color = "#2E86AB") +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Mean FEV1/FVC Ratio with Standard Deviation", y = "Mean FEV1/FVC (%)", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

Sawtooth Pattern Visualization

# Bar plot with percentages
p1 <- ggplot(sawtooth_data, aes(x = Period, y = Percent, fill = Period)) +
  geom_bar(stat = "identity", alpha = 0.7) +
  geom_text(aes(label = sprintf("%.1f%%", Percent)), vjust = -0.5, size = 4) +
  scale_fill_manual(values = c("PRE" = "#2E86AB", "PAN" = "#A23B72", "POST" = "#F18F01")) +
  labs(title = "Sawtooth Pattern Prevalence by Period", 
       y = "Percentage with Pattern (%)", x = "Period") +
  ylim(0, 70) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold"))


sawtooth_long <- sawtooth_data %>%
  select(Period, Present, Absent) %>%
  pivot_longer(cols = c(Present, Absent), names_to = "Pattern", values_to = "Count")

p2 <- ggplot(sawtooth_long, aes(x = Period, y = Count, fill = Pattern)) +
  geom_bar(stat = "identity", position = "fill") +
  scale_fill_manual(values = c("Present" = "#E63946", "Absent" = "#A8DADC")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Sawtooth Pattern Distribution", 
       y = "Proportion", x = "Period") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

grid.arrange(p1, p2, ncol = 1)

Summary of Findings

ANOVA Results Summary

anova_summary_table <- data.frame(
  Variable = c("FVC%", "FVC Absolute", "FEV1%", "FEV1 Absolute", "FEV1/FVC Ratio"),
  F_Statistic = c(
    summary(fvc_pct_anova$anova)[[1]]$`F value`[1],
    summary(fvc_abs_anova$anova)[[1]]$`F value`[1],
    summary(fev1_pct_anova$anova)[[1]]$`F value`[1],
    summary(fev1_abs_anova$anova)[[1]]$`F value`[1],
    summary(ratio_anova$anova)[[1]]$`F value`[1]
  ),
  P_Value = c(
    summary(fvc_pct_anova$anova)[[1]]$`Pr(>F)`[1],
    summary(fvc_abs_anova$anova)[[1]]$`Pr(>F)`[1],
    summary(fev1_pct_anova$anova)[[1]]$`Pr(>F)`[1],
    summary(fev1_abs_anova$anova)[[1]]$`Pr(>F)`[1],
    summary(ratio_anova$anova)[[1]]$`Pr(>F)`[1]
  ),
  Significant = c(
    summary(fvc_pct_anova$anova)[[1]]$`Pr(>F)`[1] < 0.05,
    summary(fvc_abs_anova$anova)[[1]]$`Pr(>F)`[1] < 0.05,
    summary(fev1_pct_anova$anova)[[1]]$`Pr(>F)`[1] < 0.05,
    summary(fev1_abs_anova$anova)[[1]]$`Pr(>F)`[1] < 0.05,
    summary(ratio_anova$anova)[[1]]$`Pr(>F)`[1] < 0.05
  )
)

kable(anova_summary_table, digits = 4, 
      caption = "Summary of One-Way ANOVA Results (α = 0.05)")
Summary of One-Way ANOVA Results (α = 0.05)
Variable F_Statistic P_Value Significant
FVC% 2.7066 0.0693 FALSE
FVC Absolute 0.2633 0.7688 FALSE
FEV1% 4.1131 0.0178 TRUE
FEV1 Absolute 2.2526 0.1079 FALSE
FEV1/FVC Ratio 2.8480 0.0604 FALSE