R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

#install packages

library(tidyverse)
library(readxl)
library(ggpubr)
library(dplyr)

#import_data

df <- read_excel("data.xlsx")

#normmal distribution of my data

df$logFerritin <- log(df$Ferritin)
df$logVitD     <- log(df$VitaminD)

#check normal distribution

shapiro.test(df$logFerritin) #for ferritin
## 
##  Shapiro-Wilk normality test
## 
## data:  df$logFerritin
## W = 0.96744, p-value = 0.02463
qqnorm(df$logFerritin)

qqnorm(df$logVitD)

#Ferritin_vitD ratio

df$FD_ratio <- df$Ferritin / df$VitaminD

#CRP group distribution

df <- df %>%
  mutate(CRP_Group = ifelse(nCRP < 1,
                                    "Normal",
                                    "Subclinical Inflammation"))

#vitmin D group

df$VitD_group <- cut(df$VitaminD,
                       breaks = c(-Inf, 20, 30, Inf),
                       labels = c("Deficient", "Insufficient", "Sufficient"))

#Ferritin_median split

ferritin_median <- median(df$Ferritin, na.rm = TRUE)

df$Ferritin_group <- ifelse(df$Ferritin <= ferritin_median,
                              "Low Ferritin", "High Ferritin")

#Combined Ferritin–Vitamin D Interaction Groups

df$Interaction_group <- with(df,
  ifelse(Ferritin_group == "Low Ferritin" & VitD_group == "Sufficient", "G1",
  ifelse(Ferritin_group == "Low Ferritin" & VitD_group != "Sufficient", "G2",
  ifelse(Ferritin_group == "High Ferritin" & VitD_group == "Sufficient", "G3",
         "G4"))))

df$Interaction_group <- factor(df$Interaction_group,
                                 levels = c("G1","G2","G3","G4"))

#FD Ratio Quartiles

df$FD_quartile <- ntile(df$FD_ratio, 4)
df$FD_quartile <- factor(df$FD_quartile,
                            labels = c("Q1","Q2","Q3","Q4"))

#DESCRIPTIVE STATISTICS

describeBy(df$nCRP, df$Interaction_group)
## 
##  Descriptive statistics by group 
## group: G1
##    vars n mean   sd median trimmed mad min max range skew kurtosis   se
## X1    1 7 1.74 1.55    0.8    1.74 0.3 0.6 4.8   4.2 0.94     -0.7 0.59
## ------------------------------------------------------------ 
## group: G2
##    vars  n mean   sd median trimmed  mad min max range skew kurtosis  se
## X1    1 38 1.03 0.62    0.7    0.94 0.15 0.5 2.4   1.9 1.35     0.42 0.1
## ------------------------------------------------------------ 
## group: G3
##    vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## X1    1 7  3.1 3.48    0.7     3.1 0.15 0.6 9.6     9 0.76    -1.09 1.31
## ------------------------------------------------------------ 
## group: G4
##    vars  n mean  sd median trimmed  mad min max range skew kurtosis   se
## X1    1 37 2.02 1.9    1.2    1.74 0.89 0.5 9.6   9.1 1.98     4.61 0.31
describeBy(df$FD_ratio, df$FD_quartile)
## 
##  Descriptive statistics by group 
## group: Q1
##    vars  n mean  sd median trimmed  mad  min max range skew kurtosis   se
## X1    1 23 1.22 0.9   0.98    1.18 1.01 0.03 2.8  2.76 0.37    -1.39 0.19
## ------------------------------------------------------------ 
## group: Q2
##    vars  n mean   sd median trimmed  mad  min   max range skew kurtosis   se
## X1    1 22 6.34 2.26      6    6.16 2.21 3.13 11.03   7.9 0.48    -0.63 0.48
## ------------------------------------------------------------ 
## group: Q3
##    vars  n mean   sd median trimmed mad   min   max range skew kurtosis   se
## X1    1 22 20.6 6.01  19.85   19.91 5.4 12.18 36.81 24.63 0.96     0.39 1.28
## ------------------------------------------------------------ 
## group: Q4
##    vars  n  mean     sd median trimmed   mad   min    max  range skew kurtosis
## X1    1 22 222.4 227.01  97.32  189.36 83.45 39.91 714.29 674.38 1.13    -0.23
##      se
## X1 48.4

#GROUP COMPARISONS

pairwise.t.test(df$nCRP, df$Interaction_group, p.adjust.method = "holm",pool.sd = FALSE)
## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  df$nCRP and df$Interaction_group 
## 
##    G1    G2    G3   
## G2 1.000 -     -    
## G3 1.000 0.831 -    
## G4 1.000 0.025 1.000
## 
## P value adjustment method: holm

#correlation

pairwise.t.test(df$nCRP, df$Interaction_group, p.adjust.method = "holm")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$nCRP and df$Interaction_group 
## 
##    G1    G2    G3   
## G2 0.587 -     -    
## G3 0.461 0.018 -    
## G4 0.685 0.054 0.461
## 
## P value adjustment method: holm
pairwise.t.test(df$nCRP, df$FD_quartile, p.adjust.method = "holm")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$nCRP and df$FD_quartile 
## 
##    Q1    Q2    Q3   
## Q2 0.850 -     -    
## Q3 0.256 0.238 -    
## Q4 0.042 0.032 0.694
## 
## P value adjustment method: holm