Introduction

  • EI Min (Elongation Index Minimum)
  • EI Max (Elongation Index Maximum)
  • EI Delta (Elongation Index Delta)
  • POS (Point of Sickling)

The analysis focuses on evaluating these variables across different groupings such as transfusion status, Lorrca instrument ID, and genotype. The main objective is to establish expected ranges for each variable within these groups and understand how they vary over time.

Data Loading and Cleaning

Load Libraries and Data

In this section, the necessary libraries for data manipulation, statistical analysis, and visualization are loaded. The data is then read from a CSV file and prepared for analysis by cleaning up column names and converting relevant columns to numeric.

# Load required libraries
library(tidyverse)
library(plotly)
library(viridis)
library(knitr)
library(kableExtra)
library(janitor)
library(moments)
library(lmerTest)
library(emmeans)

# Load and clean data
data <- read_csv("Omics_reformatted_2024-09-26.csv") %>%
  clean_names()

# Convert relevant columns to numeric
data <- data %>%
  mutate(across(c("ei_min_lorrca", "ei_max_lorrca", "ei_delta_lorrca", "pos_lorrca"), as.numeric))

# Filter samples based on quality control
data <- data %>%
  filter(qc_pass_lorrca == 'Yes' | 
         qc_pass_lorrca == 'Yes, PoS (fixed range)' | 
         qc_pass_lorrca == 'Yes, (fixed range)')

# 1. Data Preparation
data <- data %>%
  mutate(
    transfusion_status = case_when(
      transfusion_status %in% c("Yes", "Yes, chronic TXN", "Yes, acutely transfused", 
                                "Yes, acutely transfused?", "Yes, acutely transfused or may be chronic") ~ "Yes",
      transfusion_status == "No" ~ "No",
      is.na(transfusion_status) ~ "No",
      TRUE ~ transfusion_status
    ),
    date_of_collection = as.Date(date_of_collection),
    genotype = factor(genotype),
    transfusion_status = factor(transfusion_status)
  )

Statistical and Visualization Functions

Summary Statistics Calculation

A function is defined to compute summary statistics for each Lorrca variable, grouped by factors like genotype, transfusion status, and Lorrca instrument ID. The function provides mean, standard deviation, quartiles, skewness, kurtosis, and confidence intervals for the selected variable.

calculate_stats <- function(data, var, group_vars) {
  data %>%
    filter(!is.na(.data[[var]])) %>%
    group_by(across(all_of(group_vars))) %>%
    summarise(
      Mean = mean(.data[[var]], na.rm = TRUE),
      SD = sd(.data[[var]], na.rm = TRUE),
      Median = median(.data[[var]], na.rm = TRUE),
      Q1 = quantile(.data[[var]], 0.25, na.rm = TRUE),
      Q3 = quantile(.data[[var]], 0.75, na.rm = TRUE),
      IQR = IQR(.data[[var]], na.rm = TRUE),
      Min = min(.data[[var]], na.rm = TRUE),
      Max = max(.data[[var]], na.rm = TRUE),
      Skewness = skewness(.data[[var]], na.rm = TRUE),
      Kurtosis = kurtosis(.data[[var]], na.rm = TRUE),
      `Range (Mean ± 2SD)` = paste(round(Mean - 2*SD, 2), "to", round(Mean + 2*SD, 2)),
      `Range (Median ± 1.5IQR)` = paste(round(Median - 1.5*IQR, 2), "to", round(Median + 1.5*IQR, 2)),
      `95% CI Lower` = quantile(.data[[var]], 0.025, na.rm = TRUE),
      `95% CI Upper` = quantile(.data[[var]], 0.975, na.rm = TRUE),
      N = n(),
      .groups = 'drop'
    )
}

Visualization of Lorrca Variables

Another function generates box plots and histograms to visualize the distribution of Lorrca variables across the specified groupings.

plot_box_hist_plotly <- function(data, var, group_vars) {
  # Create a combined grouping variable
  data <- data %>%
    mutate(Group = interaction(across(all_of(group_vars)), sep = " - "))

  # Filter out NA values
  data_filtered <- data %>%
    filter(!is.na(Group), !is.na(.data[[var]]))

  unique_groups <- unique(data_filtered$Group)
  num_groups <- length(unique_groups)

  # Box plot
  p1 <- plot_ly(
    data = data_filtered, 
    x = ~Group, 
    y = ~.data[[var]], 
    type = 'box', 
    color = ~Group,
    colors = viridis::viridis_pal()(num_groups)
  ) %>%
    layout(
      title = paste("Box Plot of", var, "by", paste(group_vars, collapse = " & ")),
      xaxis = list(title = paste(group_vars, collapse = " & ")),
      yaxis = list(title = var),
      height = 600,
      width = 1000
    )

  # Histogram
  p2 <- plot_ly(
    data = data_filtered, 
    x = ~.data[[var]], 
    type = 'histogram', 
    color = ~Group,
    colors = viridis::viridis_pal()(num_groups),
    nbinsx = 30
  ) %>%
    layout(
      title = paste("Histogram of", var, "by", paste(group_vars, collapse = " & ")),
      xaxis = list(title = var),
      yaxis = list(title = "Count"),
      height = 600,
      width = 1000
    )

  return(list(boxplot = p1, histogram = p2))
}

Analysis Results

EI Min Analysis

The analysis of EI Min is performed by grouping the data by transfusion status, Lorrca instrument ID, and genotype. The results are presented as box plots and histograms, along with a table of summary statistics.

# Grouped by Transfusion Status
plots <- plot_box_hist_plotly(data, "ei_min_lorrca", "transfusion_status")
subplot(plots$boxplot, plots$histogram, nrows = 2)
# Summary statistics by transfusion status
transfusion_stats <- calculate_stats(data, "ei_min_lorrca", "transfusion_status")
kable(transfusion_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
transfusion_status Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
No 0.1558606 0.1249924 0.1250 0.069000 0.21025 0.141250 -0.030 0.582 1.1779131 4.072112 -0.09 to 0.41 -0.09 to 0.34 -0.007000 0.499000 312
Yes 0.2902024 0.1361678 0.3015 0.166375 0.38100 0.214625 0.082 0.577 0.1589615 2.020257 0.02 to 0.56 -0.02 to 0.62 0.092125 0.536525 42 
# Grouped by Lorrca ID
plots <- plot_box_hist_plotly(data, "ei_min_lorrca", "instrument_lorrca")
subplot(plots$boxplot, plots$histogram, nrows = 2)
lorrca_stats <- calculate_stats(data, "ei_min_lorrca", "instrument_lorrca")
kable(lorrca_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
instrument_lorrca Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
BO 0.2365417 0.1521385 0.1845 0.146125 0.347250 0.201125 0.0230 0.5095 0.6695694 2.291054 -0.07 to 0.54 -0.12 to 0.49 0.048025 0.5034500 12
Bilbo 0.1685769 0.0967565 0.1520 0.108000 0.209500 0.101500 0.0425 0.3990 1.0592580 3.617386 -0.02 to 0.36 0 to 0.3 0.051350 0.3691500 13
Frodo 0.1259375 0.1383934 0.0805 0.051000 0.133750 0.082750 -0.0110 0.4990 1.7533991 5.034599 -0.15 to 0.4 -0.04 to 0.2 0.000250 0.4630000 16
Gollum 0.2027165 0.1303335 0.1695 0.105250 0.287000 0.181750 -0.0015 0.5390 0.7750760 2.848659 -0.06 to 0.46 -0.1 to 0.44 0.007350 0.5067000 127
Review 0.3900000 NA 0.3900 0.390000 0.390000 0.000000 0.3900 0.3900 NaN NaN NA to NA 0.39 to 0.39 0.390000 0.3900000 1
SM 0.1796000 0.0949720 0.1550 0.121500 0.247500 0.126000 0.0700 0.3040 0.2318770 1.601000 -0.01 to 0.37 -0.03 to 0.34 0.075150 0.2983500 5
Smeagol 0.1482886 0.1333169 0.1090 0.058000 0.204000 0.146000 -0.0300 0.5820 1.2504439 4.118843 -0.12 to 0.41 -0.11 to 0.33 -0.008650 0.4965500 175
Two Machines 0.1660000 NA 0.1660 0.166000 0.166000 0.000000 0.1660 0.1660 NaN NaN NA to NA 0.17 to 0.17 0.166000 0.1660000 1
Woodruff 0.1556250 0.0648734 0.1520 0.105500 0.202125 0.096625 0.0920 0.2265 0.0845843 1.207531 0.03 to 0.29 0.01 to 0.3 0.093350 0.2240625 4 
# Grouped by Genotype
plots <- plot_box_hist_plotly(data, "ei_min_lorrca", "genotype")
subplot(plots$boxplot, plots$histogram, nrows = 2)
genotype_stats <- calculate_stats(data, "ei_min_lorrca", "genotype")
kable(genotype_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
genotype Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
SBet0/HPFH 0.3990000 NA 0.399 0.39900 0.39900 0.00000 0.399 0.399 NaN NaN NA to NA 0.4 to 0.4 0.3990 0.3990000 1
Sbeta+ 0.2510588 0.1461941 0.206 0.14700 0.30800 0.16100 0.086 0.582 1.0267240 2.928335 -0.04 to 0.54 -0.04 to 0.45 0.0956 0.5488000 17
Sbeta0 0.1262000 0.0990994 0.113 0.10000 0.16400 0.06400 -0.009 0.263 0.0351838 2.225624 -0.07 to 0.32 0.02 to 0.21 0.0019 0.2531000 5
SC 0.1586923 0.1080480 0.140 0.07450 0.22925 0.15475 -0.011 0.436 0.8070388 2.852243 -0.06 to 0.37 -0.09 to 0.37 0.0205 0.4132500 91
SS 0.1711583 0.1397933 0.136 0.06975 0.24125 0.17150 -0.030 0.577 0.9982879 3.246943 -0.11 to 0.45 -0.12 to 0.39 -0.0070 0.5097125 240

EI Max Analysis

Similarly, the same analysis is performed for EI Max.

# Grouped by Transfusion Status
plots <- plot_box_hist_plotly(data, "ei_max_lorrca", "transfusion_status")
subplot(plots$boxplot, plots$histogram, nrows = 2)
# Summary statistics by transfusion status
transfusion_stats <- calculate_stats(data, "ei_max_lorrca", "transfusion_status")
kable(transfusion_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
transfusion_status Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
No 0.4853113 0.0813055 0.4995 0.436 0.54725 0.11125 0.2030379 0.628 -0.7945631 3.272139 0.32 to 0.65 0.33 to 0.67 0.28765 0.600900 312
Yes 0.5132619 0.0842757 0.5255 0.471 0.58075 0.10975 0.3300000 0.630 -0.6900097 2.486487 0.34 to 0.68 0.36 to 0.69 0.34310 0.625825 42 
# Grouped by Lorrca ID
plots <- plot_box_hist_plotly(data, "ei_max_lorrca", "instrument_lorrca")
subplot(plots$boxplot, plots$histogram, nrows = 2)
lorrca_stats <- calculate_stats(data, "ei_max_lorrca", "instrument_lorrca")
kable(lorrca_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
instrument_lorrca Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
BO 0.5056948 0.1045175 0.5271319 0.4966152 0.5714553 0.0748401 0.2030379 0.5935962 -2.1525256 7.106316 0.3 to 0.71 0.41 to 0.64 0.2744441 0.5933476 12
Bilbo 0.5303077 0.0443366 0.5370000 0.5110000 0.5690000 0.0580000 0.4340000 0.5850000 -0.8002512 2.898024 0.44 to 0.62 0.45 to 0.62 0.4433000 0.5823000 13
Frodo 0.4587500 0.0860926 0.4905000 0.3667500 0.5172500 0.1505000 0.2890000 0.5860000 -0.5262666 2.085245 0.29 to 0.63 0.26 to 0.72 0.3118750 0.5736250 16
Gollum 0.4653780 0.0749548 0.4730000 0.4040000 0.5190000 0.1150000 0.2510000 0.5930000 -0.4970847 2.774362 0.32 to 0.62 0.3 to 0.65 0.3097000 0.5838000 127
Review 0.6260000 NA 0.6260000 0.6260000 0.6260000 0.0000000 0.6260000 0.6260000 NaN NaN NA to NA 0.63 to 0.63 0.6260000 0.6260000 1
SM 0.4867369 0.0814469 0.5285625 0.4553333 0.5438125 0.0884792 0.3574286 0.5485476 -0.8783490 2.193105 0.32 to 0.65 0.4 to 0.66 0.3672190 0.5480741 5
Smeagol 0.5038634 0.0825683 0.5220000 0.4580000 0.5635000 0.1055000 0.2520000 0.6300000 -0.9275870 3.405831 0.34 to 0.67 0.36 to 0.68 0.2935875 0.6186500 175
Two Machines 0.5660000 NA 0.5660000 0.5660000 0.5660000 0.0000000 0.5660000 0.5660000 NaN NaN NA to NA 0.57 to 0.57 0.5660000 0.5660000 1
Woodruff 0.4417500 0.0755000 0.4715000 0.4305000 0.4827500 0.0522500 0.3300000 0.4940000 -1.0644488 2.264777 0.29 to 0.59 0.39 to 0.55 0.3400500 0.4928750 4 
# Grouped by Genotype
plots <- plot_box_hist_plotly(data, "ei_max_lorrca", "genotype")
subplot(plots$boxplot, plots$histogram, nrows = 2)
genotype_stats <- calculate_stats(data, "ei_max_lorrca", "genotype")
kable(genotype_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
genotype Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
SBet0/HPFH 0.5760000 NA 0.576 0.576 0.57600 0.00000 0.5760000 0.576 NaN NaN NA to NA 0.58 to 0.58 0.576000 0.5760 1
Sbeta+ 0.5716290 0.0417604 0.573 0.553 0.60800 0.05500 0.4760000 0.628 -0.7072841 2.846051 0.49 to 0.66 0.49 to 0.66 0.489200 0.6240 17
Sbeta0 0.5066000 0.1013721 0.547 0.503 0.54700 0.04400 0.3360000 0.600 -1.0782357 2.747608 0.3 to 0.71 0.48 to 0.61 0.352700 0.5947 5
SC 0.4756949 0.0655605 0.486 0.433 0.52550 0.09250 0.3250000 0.608 -0.1874182 2.367124 0.34 to 0.61 0.35 to 0.62 0.340250 0.5885 91
SS 0.4869133 0.0862485 0.507 0.433 0.55025 0.11725 0.2030379 0.630 -0.8713696 3.234079 0.31 to 0.66 0.33 to 0.68 0.280825 0.6040 240

POS Analysis

The analysis for POS is performed similarly, grouped by transfusion status, Lorrca ID, and genotype.

# Grouped by Transfusion Status
plots <- plot_box_hist_plotly(data, "pos_lorrca", "transfusion_status")
subplot(plots$boxplot, plots$histogram, nrows = 2)
# Summary statistics by transfusion status
transfusion_stats <- calculate_stats(data, "pos_lorrca", "transfusion_status")
kable(transfusion_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
transfusion_status Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
No 39.66717 18.41171 38.4870 26.51850 50.1985 23.68000 5.620 161.6717 1.1717899 8.329228 2.84 to 76.49 2.97 to 74.01 9.906825 78.27890 310
Yes 38.46728 14.52254 39.7095 29.62275 45.3710 15.74825 6.429 75.6410 0.3929657 3.534145 9.42 to 67.51 16.09 to 63.33 14.469825 74.51682 40 
# Grouped by Lorrca ID
plots <- plot_box_hist_plotly(data, "pos_lorrca", "instrument_lorrca")
subplot(plots$boxplot, plots$histogram, nrows = 2)
lorrca_stats <- calculate_stats(data, "pos_lorrca", "instrument_lorrca")
kable(lorrca_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
instrument_lorrca Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
BO 51.30214 38.23327 40.31323 31.14589 60.25847 29.11258 15.05363 161.67175 2.1371527 7.015830 -25.16 to 127.77 -3.36 to 83.98 18.37973 135.72043 12
Bilbo 40.61815 12.75774 44.88600 34.39300 50.98300 16.59000 12.61200 54.82000 -0.7650845 2.698680 15.1 to 66.13 20 to 69.77 16.39980 54.54220 13
Frodo 57.61280 20.28234 58.51700 50.54700 76.19850 25.65150 17.18200 83.32900 -0.6354102 2.410566 17.05 to 98.18 20.04 to 96.99 19.77795 81.64165 15
Gollum 33.33913 17.64062 31.31700 21.26000 42.24300 20.98300 5.62000 91.76500 0.9018957 3.754757 -1.94 to 68.62 -0.16 to 62.79 8.26360 75.25030 125
Review 32.44400 NA 32.44400 32.44400 32.44400 0.00000 32.44400 32.44400 NaN NaN NA to NA 32.44 to 32.44 32.44400 32.44400 1
SM 41.21129 15.16564 41.48864 30.82780 50.04025 19.21245 22.70143 60.99833 0.0780604 1.703614 10.88 to 71.54 12.67 to 70.31 23.51407 59.90253 5
Smeagol 41.32234 14.50343 43.01900 30.58850 51.37725 20.78875 6.42900 86.96800 0.0887463 2.973130 12.32 to 70.33 11.84 to 74.2 14.33000 70.22190 174
Two Machines 46.94700 NA 46.94700 46.94700 46.94700 0.00000 46.94700 46.94700 NaN NaN NA to NA 46.95 to 46.95 46.94700 46.94700 1
Woodruff 46.18375 19.19058 38.84700 35.57550 49.45525 13.87975 32.55300 74.48800 1.0416910 2.246289 7.8 to 84.56 18.03 to 59.67 32.85525 71.98472 4 
# Grouped by Genotype
plots <- plot_box_hist_plotly(data, "pos_lorrca", "genotype")
subplot(plots$boxplot, plots$histogram, nrows = 2)
genotype_stats <- calculate_stats(data, "pos_lorrca", "genotype")
kable(genotype_stats) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
genotype Mean SD Median Q1 Q3 IQR Min Max Skewness Kurtosis Range (Mean ± 2SD) Range (Median ± 1.5IQR) 95% CI Lower 95% CI Upper N
SBet0/HPFH 12.61200 NA 12.612 12.61200 12.612 0.00000 12.612 12.6120 NaN NaN NA to NA 12.61 to 12.61 12.612000 12.61200 1
Sbeta+ 24.81592 10.20173 24.486 21.70225 29.543 7.84075 7.628 48.3290 0.2413311 3.336184 4.41 to 45.22 12.72 to 36.25 7.857875 43.68612 16
Sbeta0 37.55100 19.81498 29.813 25.51900 35.018 9.49900 25.136 72.2690 1.3548575 3.043054 -2.08 to 77.18 15.56 to 44.06 25.174300 68.54390 5
SC 27.81876 10.99707 27.240 19.12750 35.369 16.24150 8.467 55.0980 0.2168359 2.343858 5.82 to 49.81 2.88 to 51.6 9.813250 48.22400 91
SS 45.17547 17.79901 43.955 33.70100 54.586 20.88500 5.620 161.6717 1.2203838 9.896227 9.58 to 80.77 12.63 to 75.28 14.753400 79.44280 237

Summary of Statistical Interpretation

  • Range (Mean ± 2SD) assumes a normal distribution and covers approximately 95% of the data. This is useful when the data is symmetrically distributed.
  • Range (Median ± 1.5IQR) is more robust to outliers and represents the central 50% of the data, expanded by 1.5 times the interquartile range. This is ideal for skewed distributions.
  • Skewness indicates the asymmetry of the distribution. A positive value means the distribution is right-skewed (longer tail on the right), while a negative value suggests left skewness.
  • Kurtosis measures the “tailedness” of the distribution. High kurtosis indicates the presence of outliers or heavy tails in the data.

Modeling and Predictions Over Time

This section models EI Min, EI Max, and POS over time, considering interactions between genotype, transfusion status, and time. The models are fitted using mixed-effects regression, controlling for instrument variability (random effect on Lorrca instrument).

# 1. Data Preparation
data <- data %>%
  mutate(
    transfusion_status = case_when(
      transfusion_status %in% c("Yes", "Yes, chronic TXN", "Yes, acutely transfused", 
                                "Yes, acutely transfused?", "Yes, acutely transfused or may be chronic") ~ "Yes",
      transfusion_status == "No" ~ "No",
      is.na(transfusion_status) ~ "No",
      TRUE ~ transfusion_status
    ),
    date_of_collection = as.Date(date_of_collection),
    genotype = factor(genotype),
    transfusion_status = factor(transfusion_status)
  )

# Create time variable using the earliest date as reference
data <- data %>%
  group_by(subject_id) %>%
  mutate(
    reference_date = min(date_of_collection, na.rm = TRUE),
    time = as.numeric(difftime(date_of_collection, reference_date, units = "days"))
  ) %>%
  ungroup()

# Check the created time variable
summary(data$time)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0     0.0    68.6     0.0   678.0
# 2. Normalize Lorrca variables
lorrca_vars <- c("ei_min_lorrca", "ei_max_lorrca", "ei_delta_lorrca", "pos_lorrca")
data <- data %>%
  mutate(across(all_of(lorrca_vars), ~ (. - mean(., na.rm = TRUE)), .names = "{col}"))

# 3. Fit mixed-effects models
model_ei_min <- lmer(ei_min_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca), data = data)
model_ei_max <- lmer(ei_max_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca), data = data)
model_pos <- lmer(pos_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca), data = data)

# 4. Summarize models
summary(model_ei_min)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ei_min_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca)
##    Data: data
## 
## REML criterion at convergence: -384.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7523 -0.6867 -0.2138  0.4091  3.7334 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  instrument_lorrca (Intercept) 0.001164 0.03412 
##  Residual                      0.014637 0.12098 
## Number of obs: 354, groups:  instrument_lorrca, 9
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                 2.360e-01  1.237e-01  3.277e+02   1.908   0.0573
## genotypeSbeta+             -1.348e-01  1.281e-01  3.387e+02  -1.052   0.2934
## genotypeSbeta0             -2.761e-01  1.379e-01  3.410e+02  -2.003   0.0460
## genotypeSC                 -2.357e-01  1.244e-01  3.396e+02  -1.895   0.0589
## genotypeSS                 -2.589e-01  1.239e-01  3.396e+02  -2.090   0.0373
## transfusion_statusYes       1.526e-01  2.325e-02  3.410e+02   6.565 1.94e-10
## time                        1.368e-05  6.878e-05  3.426e+02   0.199   0.8424
## genotypeSbeta+:time        -1.792e-04  2.651e-04  3.388e+02  -0.676   0.4994
## genotypeSbeta0:time         1.397e-04  3.474e-04  3.383e+02   0.402   0.6878
## genotypeSC:time            -1.202e-04  1.061e-04  3.387e+02  -1.133   0.2580
## transfusion_statusYes:time -3.289e-05  1.532e-04  3.424e+02  -0.215   0.8301
##                               
## (Intercept)                .  
## genotypeSbeta+                
## genotypeSbeta0             *  
## genotypeSC                 .  
## genotypeSS                 *  
## transfusion_statusYes      ***
## time                          
## genotypeSbeta+:time           
## genotypeSbeta0:time           
## genotypeSC:time               
## transfusion_statusYes:time    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) gntyS+ gntyS0 gntySC gntySS trns_Y time   gntS+: gntS0:
## genotypSbt+ -0.958                                                        
## genotypSbt0 -0.890  0.867                                                 
## genotypeSC  -0.985  0.959  0.892                                          
## genotypeSS  -0.989  0.963  0.895  0.990                                   
## trnsfsn_stY -0.002  0.006  0.008  0.006 -0.026                            
## time        -0.008  0.015  0.019  0.012 -0.019  0.188                     
## gntypSbt+:t  0.003 -0.118 -0.002 -0.002  0.006 -0.043 -0.246              
## gntypSbt0:t  0.001 -0.003 -0.196 -0.003  0.004 -0.036 -0.193  0.048       
## gentypSC:tm  0.001  0.000  0.000 -0.046  0.020 -0.109 -0.608  0.159  0.122
## trnsfsn_sY: -0.001 -0.004 -0.005 -0.005  0.010 -0.457 -0.430  0.110  0.085
##             gntSC:
## genotypSbt+       
## genotypSbt0       
## genotypeSC        
## genotypeSS        
## trnsfsn_stY       
## time              
## gntypSbt+:t       
## gntypSbt0:t       
## gentypSC:tm       
## trnsfsn_sY:  0.277
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 9 columns / coefficients
summary(model_ei_max)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## ei_max_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca)
##    Data: data
## 
## REML criterion at convergence: -686.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7041 -0.5795  0.1646  0.7370  1.6689 
## 
## Random effects:
##  Groups            Name        Variance  Std.Dev.
##  instrument_lorrca (Intercept) 0.0004978 0.02231 
##  Residual                      0.0060597 0.07784 
## Number of obs: 354, groups:  instrument_lorrca, 9
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)                 6.678e-02  7.963e-02  3.244e+02   0.839   0.4023  
## genotypeSbeta+              1.773e-02  8.246e-02  3.381e+02   0.215   0.8299  
## genotypeSbeta0             -7.481e-02  8.874e-02  3.408e+02  -0.843   0.3998  
## genotypeSC                 -7.588e-02  8.006e-02  3.391e+02  -0.948   0.3439  
## genotypeSS                 -7.151e-02  7.971e-02  3.391e+02  -0.897   0.3703  
## transfusion_statusYes       2.653e-02  1.496e-02  3.406e+02   1.773   0.0771 .
## time                        2.141e-06  4.426e-05  3.424e+02   0.048   0.9614  
## genotypeSbeta+:time        -1.435e-05  1.706e-04  3.378e+02  -0.084   0.9330  
## genotypeSbeta0:time         2.722e-04  2.235e-04  3.372e+02   1.218   0.2240  
## genotypeSC:time            -8.705e-06  6.825e-05  3.378e+02  -0.128   0.8986  
## transfusion_statusYes:time  4.184e-05  9.857e-05  3.424e+02   0.424   0.6715  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) gntyS+ gntyS0 gntySC gntySS trns_Y time   gntS+: gntS0:
## genotypSbt+ -0.958                                                        
## genotypSbt0 -0.890  0.868                                                 
## genotypeSC  -0.985  0.959  0.892                                          
## genotypeSS  -0.989  0.963  0.895  0.990                                   
## trnsfsn_stY -0.002  0.006  0.008  0.006 -0.026                            
## time        -0.008  0.015  0.020  0.012 -0.019  0.188                     
## gntypSbt+:t  0.003 -0.118 -0.002 -0.002  0.006 -0.042 -0.246              
## gntypSbt0:t  0.002 -0.003 -0.196 -0.003  0.004 -0.036 -0.193  0.048       
## gentypSC:tm  0.002  0.000  0.000 -0.046  0.020 -0.109 -0.608  0.159  0.122
## trnsfsn_sY: -0.001 -0.004 -0.005 -0.005  0.010 -0.457 -0.430  0.110  0.085
##             gntSC:
## genotypSbt+       
## genotypSbt0       
## genotypeSC        
## genotypeSS        
## trnsfsn_stY       
## time              
## gntypSbt+:t       
## gntypSbt0:t       
## gentypSC:tm       
## trnsfsn_sY:  0.277
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 9 columns / coefficients
summary(model_pos)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## pos_lorrca ~ genotype * transfusion_status * time + (1 | instrument_lorrca)
##    Data: data
## 
## REML criterion at convergence: 2889.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6710 -0.5596 -0.0014  0.4392  7.1607 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  instrument_lorrca (Intercept)  52.89    7.273  
##  Residual                      222.39   14.913  
## Number of obs: 350, groups:  instrument_lorrca, 9
## 
## Fixed effects:
##                              Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)                -25.628016  15.535598 319.285077  -1.650  0.10000   
## genotypeSbeta+              12.813122  15.980305 338.979131   0.802  0.42323   
## genotypeSbeta0              30.816204  17.127922 338.757414   1.799  0.07288 . 
## genotypeSC                  16.809607  15.467143 338.953587   1.087  0.27790   
## genotypeSS                  36.373570  15.400722 338.953886   2.362  0.01875 * 
## transfusion_statusYes       -8.286976   2.946876 335.653057  -2.812  0.00521 **
## time                         0.002569   0.008547 335.888391   0.301  0.76388   
## genotypeSbeta+:time          0.012334   0.032940 333.727185   0.374  0.70832   
## genotypeSbeta0:time         -0.046277   0.042826 333.637056  -1.081  0.28066   
## genotypeSC:time              0.014886   0.013097 333.983158   1.137  0.25652   
## transfusion_statusYes:time   0.002992   0.019099 338.976430   0.157  0.87559   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) gntyS+ gntyS0 gntySC gntySS trns_Y time   gntS+: gntS0:
## genotypSbt+ -0.946                                                        
## genotypSbt0 -0.882  0.868                                                 
## genotypeSC  -0.975  0.957  0.893                                          
## genotypeSS  -0.979  0.961  0.897  0.990                                   
## trnsfsn_stY -0.002  0.006  0.009  0.007 -0.024                            
## time        -0.011  0.020  0.024  0.016 -0.014  0.183                     
## gntypSbt+:t  0.004 -0.127 -0.004 -0.004  0.004 -0.041 -0.246              
## gntypSbt0:t  0.002 -0.004 -0.195 -0.003  0.003 -0.035 -0.194  0.049       
## gentypSC:tm  0.003 -0.001  0.000 -0.046  0.019 -0.106 -0.607  0.159  0.122
## trnsfsn_sY: -0.005 -0.004 -0.005 -0.005  0.009 -0.464 -0.426  0.109  0.085
##             gntSC:
## genotypSbt+       
## genotypSbt0       
## genotypeSC        
## genotypeSS        
## trnsfsn_stY       
## time              
## gntypSbt+:t       
## gntypSbt0:t       
## gentypSC:tm       
## trnsfsn_sY:  0.276
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 9 columns / coefficients

Predictions Over Time

We create a prediction grid to simulate the Lorrca variables’ expected values over a span of two years (730 days), divided into time intervals of 6 months. Predictions are made for each combination of genotype and transfusion status using the fitted models.

# 5. Create prediction grid
genotypes <- levels(data$genotype)
transfusion_statuses <- levels(data$transfusion_status)
timepoints <- seq(0, 730, by = 182.5)  # Every 6 months over 2 years

pred_grid <- expand.grid(
  genotype = genotypes,
  transfusion_status = transfusion_statuses,
  time = timepoints,
  instrument_lorrca = factor("default_instrument", levels = levels(data$instrument_lorrca))
)

# 6. Make predictions
pred_grid$pred_ei_min <- predict(model_ei_min, newdata = pred_grid, re.form = NA)
pred_grid$pred_ei_max <- predict(model_ei_max, newdata = pred_grid, re.form = NA)
pred_grid$pred_pos <- predict(model_pos, newdata = pred_grid, re.form = NA)

# 7. Reshape data for plotting
pred_long <- pred_grid %>%
  pivot_longer(cols = c(pred_ei_min, pred_ei_max, pred_pos),
               names_to = "measure", 
               values_to = "value") %>%
  mutate(measure = factor(measure, 
                          levels = c("pred_ei_min", "pred_ei_max", "pred_pos"),
                          labels = c("EI Min", "EI Max", "POS")))

# 8. Create the plot
plot <- ggplot(pred_long, aes(x = time, y = value, color = genotype, 
                              linetype = transfusion_status, group = interaction(genotype, transfusion_status))) +
  geom_line() +
  geom_point() +
  facet_wrap(~ measure, scales = "free_y", ncol = 1) +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Predicted Lorrca Values Over Time",
       subtitle = "By Genotype and Transfusion Status",
       x = "Time (days)", 
       y = "Predicted Value",
       color = "Genotype",
       linetype = "Transfusion Status") +
  theme_minimal() +
  theme(legend.position = "right",
        strip.background = element_rect(fill = "lightgrey"),
        strip.text = element_text(face = "bold"))

# 9. Save and display the plot
ggsave("predicted_lorrca_values_plot.png", plot, width = 12, height = 15, dpi = 300)
print(plot)

Additional Statistical Tests

To further evaluate the models, Type III ANOVA is conducted on the fixed effects, and estimated marginal means (EMMs) are calculated for genotype, transfusion status, and their interactions. These tests help to assess the effects of each variable and how they interact with time.

analyze_fixed_effects <- function(model, model_name) {
  cat("\n\n===", model_name, "===\n")
  
  cat("\nType III ANOVA of Fixed Effects:\n")
  print(anova(model, type = 3))
  
  cat("\nEstimated Marginal Means for Genotype:\n")
  print(emmeans(model, ~ genotype))
  
  cat("\nPairwise Comparisons of Genotypes:\n")
  print(emmeans(model, ~ genotype) |> pairs())
  
  cat("\nEstimated Marginal Means for Transfusion Status:\n")
  print(emmeans(model, ~ transfusion_status))
  
  cat("\nInteraction between Genotype and Transfusion Status:\n")
  print(emmeans(model, ~ genotype * transfusion_status) |> pairs())
  
  cat("\nEffect of Time:\n")
  print(summary(model)$coefficients["time",])
}

# Run additional tests for each model
analyze_fixed_effects(model_ei_min, "EI Min Model")
## 
## 
## === EI Min Model ===
## 
## Type III ANOVA of Fixed Effects:
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)
## genotype                         0.26130 0.06532     4 339.66  4.4629 0.001581
## transfusion_status               0.63084 0.63084     1 341.04 43.0985 1.94e-10
## time                             0.00131 0.00131     1 341.28  0.0896 0.764862
## genotype:time                    0.02706 0.00902     3 338.62  0.6162 0.604915
## transfusion_status:time          0.00067 0.00067     1 342.45  0.0461 0.830138
## genotype:transfusion_status                                                   
## genotype:transfusion_status:time                                              
##                                     
## genotype                         ** 
## transfusion_status               ***
## time                                
## genotype:time                       
## transfusion_status:time             
## genotype:transfusion_status         
## genotype:transfusion_status:time    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimated Marginal Means for Genotype:
##  genotype   emmean   SE   df lower.CL upper.CL
##  SBet0/HPFH nonEst   NA   NA       NA       NA
##  Sbeta+     nonEst   NA   NA       NA       NA
##  Sbeta0     nonEst   NA   NA       NA       NA
##  SC         nonEst   NA   NA       NA       NA
##  SS         0.0532 0.02 7.41  0.00642      0.1
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Pairwise Comparisons of Genotypes:
##  contrast                estimate SE df z.ratio p.value
##  (SBet0/HPFH) - (Sbeta+)   nonEst NA NA      NA      NA
##  (SBet0/HPFH) - Sbeta0     nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SC         nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SS         nonEst NA NA      NA      NA
##  (Sbeta+) - Sbeta0         nonEst NA NA      NA      NA
##  (Sbeta+) - SC             nonEst NA NA      NA      NA
##  (Sbeta+) - SS             nonEst NA NA      NA      NA
##  Sbeta0 - SC               nonEst NA NA      NA      NA
##  Sbeta0 - SS               nonEst NA NA      NA      NA
##  SC - SS                   nonEst NA NA      NA      NA
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 2 estimates 
## 
## Estimated Marginal Means for Transfusion Status:
##  transfusion_status emmean SE df asymp.LCL asymp.UCL
##  No                 nonEst NA NA        NA        NA
##  Yes                nonEst NA NA        NA        NA
## 
## Results are averaged over the levels of: genotype 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Interaction between Genotype and Transfusion Status:
##  contrast                           estimate     SE  df t.ratio p.value
##  (SBet0/HPFH No) - (Sbeta+ No)        nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 No          nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - SC No              nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - SS No              nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - (SBet0/HPFH Yes)   nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - (Sbeta+ Yes)       nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 Yes         nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - SC Yes             nonEst     NA  NA      NA      NA
##  (SBet0/HPFH No) - SS Yes             nonEst     NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 No             0.11942 0.0619 339   1.931  0.3030
##  (Sbeta+ No) - SC No                 0.09684 0.0323 340   2.998  0.0242
##  (Sbeta+ No) - SS No                 0.11180 0.0308 340   3.628  0.0030
##  (Sbeta+ No) - (SBet0/HPFH Yes)       nonEst     NA  NA      NA      NA
##  (Sbeta+ No) - (Sbeta+ Yes)           nonEst     NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 Yes             nonEst     NA  NA      NA      NA
##  (Sbeta+ No) - SC Yes                 nonEst     NA  NA      NA      NA
##  (Sbeta+ No) - SS Yes               -0.03859 0.0350 340  -1.102  0.8055
##  Sbeta0 No - SC No                  -0.02258 0.0559 338  -0.404  0.9944
##  Sbeta0 No - SS No                  -0.00762 0.0551 338  -0.138  0.9999
##  Sbeta0 No - (SBet0/HPFH Yes)         nonEst     NA  NA      NA      NA
##  Sbeta0 No - (Sbeta+ Yes)             nonEst     NA  NA      NA      NA
##  Sbeta0 No - Sbeta0 Yes               nonEst     NA  NA      NA      NA
##  Sbeta0 No - SC Yes                   nonEst     NA  NA      NA      NA
##  Sbeta0 No - SS Yes                 -0.15801 0.0574 338  -2.751  0.0489
##  SC No - SS No                       0.01496 0.0155 340   0.968  0.8694
##  SC No - (SBet0/HPFH Yes)             nonEst     NA  NA      NA      NA
##  SC No - (Sbeta+ Yes)                 nonEst     NA  NA      NA      NA
##  SC No - Sbeta0 Yes                   nonEst     NA  NA      NA      NA
##  SC No - SC Yes                       nonEst     NA  NA      NA      NA
##  SC No - SS Yes                     -0.13543 0.0227 340  -5.958  <.0001
##  SS No - (SBet0/HPFH Yes)             nonEst     NA  NA      NA      NA
##  SS No - (Sbeta+ Yes)                 nonEst     NA  NA      NA      NA
##  SS No - Sbeta0 Yes                   nonEst     NA  NA      NA      NA
##  SS No - SC Yes                       nonEst     NA  NA      NA      NA
##  SS No - SS Yes                     -0.15039 0.0207 341  -7.261  <.0001
##  (SBet0/HPFH Yes) - (Sbeta+ Yes)      nonEst     NA  NA      NA      NA
##  (SBet0/HPFH Yes) - Sbeta0 Yes        nonEst     NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SC Yes            nonEst     NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SS Yes            nonEst     NA  NA      NA      NA
##  (Sbeta+ Yes) - Sbeta0 Yes            nonEst     NA  NA      NA      NA
##  (Sbeta+ Yes) - SC Yes                nonEst     NA  NA      NA      NA
##  (Sbeta+ Yes) - SS Yes                nonEst     NA  NA      NA      NA
##  Sbeta0 Yes - SC Yes                  nonEst     NA  NA      NA      NA
##  Sbeta0 Yes - SS Yes                  nonEst     NA  NA      NA      NA
##  SC Yes - SS Yes                      nonEst     NA  NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## 
## Effect of Time:
##     Estimate   Std. Error           df      t value     Pr(>|t|) 
## 1.368223e-05 6.877785e-05 3.425877e+02 1.989337e-01 8.424326e-01
analyze_fixed_effects(model_ei_max, "EI Max Model")
## 
## 
## === EI Max Model ===
## 
## Type III ANOVA of Fixed Effects:
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq   Mean Sq NumDF  DenDF F value
## genotype                         0.108230 0.0270576     4 338.98  4.4652
## transfusion_status               0.019051 0.0190514     1 340.56  3.1440
## time                             0.006453 0.0064533     1 340.87  1.0650
## genotype:time                    0.009520 0.0031735     3 337.65  0.5237
## transfusion_status:time          0.001092 0.0010918     1 342.35  0.1802
## genotype:transfusion_status                                             
## genotype:transfusion_status:time                                        
##                                    Pr(>F)   
## genotype                         0.001575 **
## transfusion_status               0.077102 . 
## time                             0.302820   
## genotype:time                    0.666261   
## transfusion_status:time          0.671496   
## genotype:transfusion_status                 
## genotype:transfusion_status:time            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimated Marginal Means for Genotype:
##  genotype   emmean    SE   df lower.CL upper.CL
##  SBet0/HPFH nonEst    NA   NA       NA       NA
##  Sbeta+     nonEst    NA   NA       NA       NA
##  Sbeta0     nonEst    NA   NA       NA       NA
##  SC         nonEst    NA   NA       NA       NA
##  SS         0.0101 0.013 7.44  -0.0202   0.0404
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Pairwise Comparisons of Genotypes:
##  contrast                estimate SE df z.ratio p.value
##  (SBet0/HPFH) - (Sbeta+)   nonEst NA NA      NA      NA
##  (SBet0/HPFH) - Sbeta0     nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SC         nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SS         nonEst NA NA      NA      NA
##  (Sbeta+) - Sbeta0         nonEst NA NA      NA      NA
##  (Sbeta+) - SC             nonEst NA NA      NA      NA
##  (Sbeta+) - SS             nonEst NA NA      NA      NA
##  Sbeta0 - SC               nonEst NA NA      NA      NA
##  Sbeta0 - SS               nonEst NA NA      NA      NA
##  SC - SS                   nonEst NA NA      NA      NA
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 2 estimates 
## 
## Estimated Marginal Means for Transfusion Status:
##  transfusion_status emmean SE df asymp.LCL asymp.UCL
##  No                 nonEst NA NA        NA        NA
##  Yes                nonEst NA NA        NA        NA
## 
## Results are averaged over the levels of: genotype 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Interaction between Genotype and Transfusion Status:
##  contrast                           estimate      SE  df t.ratio p.value
##  (SBet0/HPFH No) - (Sbeta+ No)        nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 No          nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - SC No              nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - SS No              nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - (SBet0/HPFH Yes)   nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - (Sbeta+ Yes)       nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 Yes         nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - SC Yes             nonEst      NA  NA      NA      NA
##  (SBet0/HPFH No) - SS Yes             nonEst      NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 No             0.07288 0.03980 338   1.831  0.3573
##  (Sbeta+ No) - SC No                 0.09322 0.02078 340   4.486  0.0001
##  (Sbeta+ No) - SS No                 0.08825 0.01983 340   4.451  0.0001
##  (Sbeta+ No) - (SBet0/HPFH Yes)       nonEst      NA  NA      NA      NA
##  (Sbeta+ No) - (Sbeta+ Yes)           nonEst      NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 Yes             nonEst      NA  NA      NA      NA
##  (Sbeta+ No) - SC Yes                 nonEst      NA  NA      NA      NA
##  (Sbeta+ No) - SS Yes                0.05885 0.02253 340   2.612  0.0704
##  Sbeta0 No - SC No                   0.02034 0.03598 338   0.565  0.9799
##  Sbeta0 No - SS No                   0.01537 0.03546 338   0.433  0.9926
##  Sbeta0 No - (SBet0/HPFH Yes)         nonEst      NA  NA      NA      NA
##  Sbeta0 No - (Sbeta+ Yes)             nonEst      NA  NA      NA      NA
##  Sbeta0 No - Sbeta0 Yes               nonEst      NA  NA      NA      NA
##  Sbeta0 No - SC Yes                   nonEst      NA  NA      NA      NA
##  Sbeta0 No - SS Yes                 -0.01403 0.03696 338  -0.380  0.9956
##  SC No - SS No                      -0.00497 0.00994 340  -0.500  0.9873
##  SC No - (SBet0/HPFH Yes)             nonEst      NA  NA      NA      NA
##  SC No - (Sbeta+ Yes)                 nonEst      NA  NA      NA      NA
##  SC No - Sbeta0 Yes                   nonEst      NA  NA      NA      NA
##  SC No - SC Yes                       nonEst      NA  NA      NA      NA
##  SC No - SS Yes                     -0.03437 0.01463 340  -2.350  0.1320
##  SS No - (SBet0/HPFH Yes)             nonEst      NA  NA      NA      NA
##  SS No - (Sbeta+ Yes)                 nonEst      NA  NA      NA      NA
##  SS No - Sbeta0 Yes                   nonEst      NA  NA      NA      NA
##  SS No - SC Yes                       nonEst      NA  NA      NA      NA
##  SS No - SS Yes                     -0.02940 0.01333 341  -2.206  0.1800
##  (SBet0/HPFH Yes) - (Sbeta+ Yes)      nonEst      NA  NA      NA      NA
##  (SBet0/HPFH Yes) - Sbeta0 Yes        nonEst      NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SC Yes            nonEst      NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SS Yes            nonEst      NA  NA      NA      NA
##  (Sbeta+ Yes) - Sbeta0 Yes            nonEst      NA  NA      NA      NA
##  (Sbeta+ Yes) - SC Yes                nonEst      NA  NA      NA      NA
##  (Sbeta+ Yes) - SS Yes                nonEst      NA  NA      NA      NA
##  Sbeta0 Yes - SC Yes                  nonEst      NA  NA      NA      NA
##  Sbeta0 Yes - SS Yes                  nonEst      NA  NA      NA      NA
##  SC Yes - SS Yes                      nonEst      NA  NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## 
## Effect of Time:
##     Estimate   Std. Error           df      t value     Pr(>|t|) 
## 2.141303e-06 4.425801e-05 3.424158e+02 4.838228e-02 9.614398e-01
analyze_fixed_effects(model_pos, "POS Model")
## 
## 
## === POS Model ===
## 
## Type III ANOVA of Fixed Effects:
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## genotype                         23011.6  5752.9     4 335.62 25.8680 < 2.2e-16
## transfusion_status                1758.7  1758.7     1 335.65  7.9080  0.005211
## time                                 0.2     0.2     1 336.69  0.0007  0.978582
## genotype:time                      635.2   211.7     3 333.80  0.9521  0.415587
## transfusion_status:time              5.5     5.5     1 338.98  0.0245  0.875589
## genotype:transfusion_status                                                    
## genotype:transfusion_status:time                                               
##                                     
## genotype                         ***
## transfusion_status               ** 
## time                                
## genotype:time                       
## transfusion_status:time             
## genotype:transfusion_status         
## genotype:transfusion_status:time    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Estimated Marginal Means for Genotype:
##  genotype   emmean   SE   df lower.CL upper.CL
##  SBet0/HPFH nonEst   NA   NA       NA       NA
##  Sbeta+     nonEst   NA   NA       NA       NA
##  Sbeta0     nonEst   NA   NA       NA       NA
##  SC         nonEst   NA   NA       NA       NA
##  SS           6.88 3.32 7.92   -0.781     14.5
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Pairwise Comparisons of Genotypes:
##  contrast                estimate SE df z.ratio p.value
##  (SBet0/HPFH) - (Sbeta+)   nonEst NA NA      NA      NA
##  (SBet0/HPFH) - Sbeta0     nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SC         nonEst NA NA      NA      NA
##  (SBet0/HPFH) - SS         nonEst NA NA      NA      NA
##  (Sbeta+) - Sbeta0         nonEst NA NA      NA      NA
##  (Sbeta+) - SC             nonEst NA NA      NA      NA
##  (Sbeta+) - SS             nonEst NA NA      NA      NA
##  Sbeta0 - SC               nonEst NA NA      NA      NA
##  Sbeta0 - SS               nonEst NA NA      NA      NA
##  SC - SS                   nonEst NA NA      NA      NA
## 
## Results are averaged over the levels of: transfusion_status 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 2 estimates 
## 
## Estimated Marginal Means for Transfusion Status:
##  transfusion_status emmean SE df asymp.LCL asymp.UCL
##  No                 nonEst NA NA        NA        NA
##  Yes                nonEst NA NA        NA        NA
## 
## Results are averaged over the levels of: genotype 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## Interaction between Genotype and Transfusion Status:
##  contrast                           estimate   SE  df t.ratio p.value
##  (SBet0/HPFH No) - (Sbeta+ No)        nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 No          nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - SC No              nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - SS No              nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - (SBet0/HPFH Yes)   nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - (Sbeta+ Yes)       nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - Sbeta0 Yes         nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - SC Yes             nonEst   NA  NA      NA      NA
##  (SBet0/HPFH No) - SS Yes             nonEst   NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 No             -13.973 7.67 333  -1.823  0.3622
##  (Sbeta+ No) - SC No                  -4.172 4.07 334  -1.024  0.8440
##  (Sbeta+ No) - SS No                 -22.712 3.90 333  -5.830  <.0001
##  (Sbeta+ No) - (SBet0/HPFH Yes)       nonEst   NA  NA      NA      NA
##  (Sbeta+ No) - (Sbeta+ Yes)           nonEst   NA  NA      NA      NA
##  (Sbeta+ No) - Sbeta0 Yes             nonEst   NA  NA      NA      NA
##  (Sbeta+ No) - SC Yes                 nonEst   NA  NA      NA      NA
##  (Sbeta+ No) - SS Yes                -14.631 4.43 334  -3.300  0.0094
##  Sbeta0 No - SC No                     9.801 6.89 333   1.422  0.6139
##  Sbeta0 No - SS No                    -8.740 6.79 333  -1.287  0.6997
##  Sbeta0 No - (SBet0/HPFH Yes)         nonEst   NA  NA      NA      NA
##  Sbeta0 No - (Sbeta+ Yes)             nonEst   NA  NA      NA      NA
##  Sbeta0 No - Sbeta0 Yes               nonEst   NA  NA      NA      NA
##  Sbeta0 No - SC Yes                   nonEst   NA  NA      NA      NA
##  Sbeta0 No - SS Yes                   -0.658 7.10 333  -0.093  1.0000
##  SC No - SS No                       -18.540 1.91 334  -9.711  <.0001
##  SC No - (SBet0/HPFH Yes)             nonEst   NA  NA      NA      NA
##  SC No - (Sbeta+ Yes)                 nonEst   NA  NA      NA      NA
##  SC No - Sbeta0 Yes                   nonEst   NA  NA      NA      NA
##  SC No - SC Yes                       nonEst   NA  NA      NA      NA
##  SC No - SS Yes                      -10.459 2.85 335  -3.664  0.0027
##  SS No - (SBet0/HPFH Yes)             nonEst   NA  NA      NA      NA
##  SS No - (Sbeta+ Yes)                 nonEst   NA  NA      NA      NA
##  SS No - Sbeta0 Yes                   nonEst   NA  NA      NA      NA
##  SS No - SC Yes                       nonEst   NA  NA      NA      NA
##  SS No - SS Yes                        8.081 2.61 337   3.091  0.0183
##  (SBet0/HPFH Yes) - (Sbeta+ Yes)      nonEst   NA  NA      NA      NA
##  (SBet0/HPFH Yes) - Sbeta0 Yes        nonEst   NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SC Yes            nonEst   NA  NA      NA      NA
##  (SBet0/HPFH Yes) - SS Yes            nonEst   NA  NA      NA      NA
##  (Sbeta+ Yes) - Sbeta0 Yes            nonEst   NA  NA      NA      NA
##  (Sbeta+ Yes) - SC Yes                nonEst   NA  NA      NA      NA
##  (Sbeta+ Yes) - SS Yes                nonEst   NA  NA      NA      NA
##  Sbeta0 Yes - SC Yes                  nonEst   NA  NA      NA      NA
##  Sbeta0 Yes - SS Yes                  nonEst   NA  NA      NA      NA
##  SC Yes - SS Yes                      nonEst   NA  NA      NA      NA
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## 
## Effect of Time:
##     Estimate   Std. Error           df      t value     Pr(>|t|) 
## 2.569461e-03 8.546990e-03 3.358884e+02 3.006275e-01 7.638845e-01