Data

This is the dataset of women where we measure equilibrium values and body composition values.

dato=read.csv("Data_Gervazoni_Natacha.csv",sep=";",dec=",",stringsAsFactors = TRUE,row.names = "ID") %>%
  dplyr::select(-"X")

nutrition=colnames(dato)[1:15]
equilibrium=colnames(dato)[16:33]

tutte.variabili=c(nutrition,equilibrium)

dat=dato %>%
  na.omit()

Equilibrium variables:

equilibrium
##  [1] "FT.DELTA.desl"   "FT.DELTA.deslap" "FT.DELTA.deslml" "FT.DELTA.vel"   
##  [5] "FT.DELTA.vel.ap" "FT.DELTA.vel.ml" "FT.DELTA.area"   "FT.DELTA.ampap" 
##  [9] "FT.DELTA.ampml"  "ST.DELTA.desl"   "ST.DELTA.deslap" "ST.DELTA.deslml"
## [13] "ST.DELTA.vel"    "ST.DELTA.vel.ap" "ST.DELTA.vel.ml" "ST.DELTA.area"  
## [17] "ST.DELTA.ampap"  "ST.DELTA.ampml"

Nutritional variables:

nutrition
##  [1] "Age"               "Weight"            "Height"           
##  [4] "BMI"               "Calf.C"            "Arm.C"            
##  [7] "Waist.C"           "Handgrip"          "PhA"              
## [10] "Z"                 "Zsp"               "Flexion.Highest"  
## [13] "Flexion.Media"     "Extension.Highest" "Ext.Media"

Here is a summary of the data set removing missing data on the variables of interest. Values are reported for all visits.

variabili_numeriche <- dat %>%
  dplyr::select(where(is.numeric)) %>%
  names()

# Overall summary for numerical variables
totali_num <- dat %>%
  summarise(across(all_of(variabili_numeriche),
                   list(mean = ~mean(.x, na.rm = TRUE),
                        sd   = ~sd(.x, na.rm = TRUE),
                        n    = ~sum(!is.na(.x)),
                        min = ~min(.x, na.rm = TRUE),
                        max=~max(.x, na.rm = TRUE)))) %>%
  pivot_longer(cols = everything(),
               names_to = c("var", ".value"),
               names_sep = "_")
# Set options to avoid scientific notation
options(scipen = 999)
tabella_finale_num=totali_num %>% 
  mutate(across(c(mean, sd,min,max), ~round(., 1))) %>% 
  mutate(across(c(n), as.integer)) %>% 
  dplyr::select(var, mean, sd, min,max,n)
# Display the final table
kableExtra::kable(tabella_finale_num)
var mean sd min max n
Age 69.5 5.9 60.1 87.6 112
Weight 65.5 11.2 40.7 95.6 112
Height 156.9 6.4 140.0 176.0 112
BMI 26.6 4.0 17.9 37.8 112
Calf.C 35.8 3.2 29.0 46.0 112
Arm.C 30.6 3.7 23.0 44.0 112
Waist.C 90.3 10.0 57.0 116.5 112
Handgrip 20.3 5.0 2.0 38.0 112
PhA 5.9 0.8 3.6 8.4 112
Z 387.4 42.0 265.7 493.5 112
Zsp 508.4 72.4 357.8 714.9 112
Flexion.Highest 58.5 14.9 21.7 103.7 112
Flexion.Media 56.1 14.6 19.9 101.5 112
Extension.Highest 94.1 24.7 45.2 164.2 112
Ext.Media 89.5 24.8 43.0 163.5 112
FT.DELTA.desl 2443.8 4667.7 -7539.9 13109.2 112
FT.DELTA.deslap 1308.0 4593.5 -9944.8 16838.8 112
FT.DELTA.deslml 1816.2 3246.2 -4827.4 10188.4 112
FT.DELTA.vel 2.2 2.3 -3.9 9.9 112
FT.DELTA.vel.ap 1.7 1.5 -2.1 8.3 112
FT.DELTA.vel.ml 1.1 1.9 -3.2 6.6 112
FT.DELTA.area 62.3 117.5 -194.8 422.4 112
FT.DELTA.ampap 2.2 6.5 -11.4 26.9 112
FT.DELTA.ampml 2.2 4.9 -8.1 16.2 112
ST.DELTA.desl 2416.2 4644.4 -17605.3 17867.0 112
ST.DELTA.deslap 1847.5 3560.3 -11430.4 9401.1 112
ST.DELTA.deslml 1229.3 3383.1 -13743.5 13776.1 112
ST.DELTA.vel 2.6 2.9 -7.2 10.7 112
ST.DELTA.vel.ap 1.4 1.8 -3.8 6.2 112
ST.DELTA.vel.ml 1.8 2.4 -7.4 7.7 112
ST.DELTA.area 60.7 118.0 -486.7 491.6 112
ST.DELTA.ampap 1.9 5.2 -20.2 14.1 112
ST.DELTA.ampml 1.6 5.3 -16.8 11.8 112

The response variables are the equilibrium variables and we want to investigate how body composition affect them.

Principal component analysis on equilibrium variables

We perform a principal component analysis on the equilibrium variables to reduce the number of response variables and then project the body composition variables on the principal components.

pca_equilibrium <- prcomp(dat %>% dplyr::select(all_of(equilibrium)), center = TRUE, scale. = TRUE)
summary(pca_equilibrium)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.7220 1.7936 1.4201 1.20489 0.99077 0.85130 0.73392
## Proportion of Variance 0.4116 0.1787 0.1120 0.08065 0.05453 0.04026 0.02992
## Cumulative Proportion  0.4116 0.5904 0.7024 0.78304 0.83757 0.87784 0.90776
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.70693 0.56573 0.51099 0.48618 0.46198 0.28316 0.19844
## Proportion of Variance 0.02776 0.01778 0.01451 0.01313 0.01186 0.00445 0.00219
## Cumulative Proportion  0.93552 0.95330 0.96781 0.98094 0.99280 0.99725 0.99944
##                           PC15    PC16    PC17    PC18
## Standard deviation     0.06251 0.05909 0.03735 0.03549
## Proportion of Variance 0.00022 0.00019 0.00008 0.00007
## Cumulative Proportion  0.99966 0.99985 0.99993 1.00000

The first two principal components explain about 60% of the total variance. We can therefore use these two components as response variables in our analysis.

# Add the principal component scores to the original dataset
dat_pca <- dat %>%
  mutate(PC1 = pca_equilibrium$x[, 1],
         PC2 = pca_equilibrium$x[, 2])

Let’s plot the correlation circle to visualize the most important equilibrium variable into the PC according to the COS2 values and then project the body composition variables.

# Calculate coordinates for the correlation circle
correlation_circle <- as.data.frame(pca_equilibrium$rotation[, 1:2])
correlation_circle$Variable <- rownames(correlation_circle)
correlation_circle <- correlation_circle %>%
  mutate(COS2 = correlation_circle$PC1^2 + correlation_circle$PC2^2)

# Calculate the correlation between PC and body composition variables
body_composition_vars <- dat %>%
  dplyr::select(all_of(nutrition))
pc_scores <- dat_pca %>%
  dplyr::select(PC1, PC2)
correlations <- map_dfr(names(body_composition_vars), function(var) {
  cor_PC1 <- cor(body_composition_vars[[var]], pc_scores$PC1)
  cor_PC2 <- cor(body_composition_vars[[var]], pc_scores$PC2)
  data.frame(PC1 = cor_PC1, PC2 = cor_PC2,Variable = var,COS2 = cor_PC1^2 + cor_PC2^2)
})


# Add the correlations to the correlation circle
correlation_circle <- rbind(correlation_circle, correlations)

# Plot of the correlation circle with body composition variables projected in red to differentiate them from the others. Red-to-green gradient is applied for COS2 values of the body composition and equilibrium variables.
ggplot(correlation_circle, aes(x = PC1, y = PC2, label = Variable)) +
  geom_segment(aes(xend = 0, yend = 0), arrow = arrow(length = unit(0.2, "cm")), color = ifelse(correlation_circle$Variable %in% nutrition, "red", "blue")) +
#  geom_text(aes(color = ifelse(Variable %in% nutrition, COS2, NA)), hjust = 1.2, vjust = 1.2,size=1.5) +
  scale_color_gradient(low = "green", high = "red", na.value = "black") +
  xlim(-0.5, 0.2) + ylim(-0.5, 0.5) +
  labs(title = "Correlation Circle with Body Composition Variables",
       x = "PC1",
       y = "PC2",
       color = "COS2") +
  theme_minimal() +
  theme(legend.position = "right")
## Ignoring unknown labels:
## • colour : "COS2"

Orthogonality between body composition variables and equilibrium variables is a sign of independence. Some body composition variables show high correlation with the first two principal components of equilibrium.

Regression analysis

Here we regress all equilibrium variable jointly against all body composition variables and report the coefficients along with their significance. The regression is the vector of equilibrium variables as response and all body composition variables as predictors.

For each equilibrium variable we calculate the median probability model using BayesVarSel and then we fit the regression for the equilibrium variable using only the selected body composition variables that form the median probability model.

The table below shows the body composition variables that significantly affect each equilibrium variable, along with their coefficients and p-values. These relations are in the so-called Median Probability Model, the model in which a regressor is more probable to affect the response than not to affect it.

# Combine all results into a single data frame
final_results <- bind_rows(results_list)
# Filter results for significant variables (p-value < 0.5) excluding the intercept
significant_results <- final_results %>%
  filter(Variable != "(Intercept)" & `Pr(>|t|)` < 0.5) %>%
  arrange(Equilibrium_Variable, `Pr(>|t|)`)
# Display the significant results
kableExtra::kable(significant_results, digits = 4, caption = "Significant Body Composition Variables Affecting Equilibrium Variables in some Median Probability Models")
Significant Body Composition Variables Affecting Equilibrium Variables in some Median Probability Models
Estimate Std. Error t value Pr(>|t|) Variable Equilibrium_Variable
Calf.C…1 0.3713 0.1591 2.3337 0.0214 Calf.C FT.DELTA.ampml
Zsp 0.0135 0.0069 1.9428 0.0546 Zsp FT.DELTA.ampml
Calf.C…3 0.2296 0.0658 3.4913 0.0007 Calf.C FT.DELTA.vel

This analysis helps to identify which aspects of body composition are most relevant for maintaining equilibrium.

Estimation of the Median Probability models

Here we calculate the regression models for the corresponding equilibrium variable using only those covariates in the median probability model. We report coefficients and model diagnostics. These models globally explain very little variability of the response.

The first is FT.DELTA.ampml regressed on Zsp and Calf.C
# Fit the regression model for FT.DELTA.ampml using only Zsp and Calf.C
model_FT_DELTA_ampml <- lm(FT.DELTA.ampml ~ Zsp + Calf.C, data = dat)
summary(model_FT_DELTA_ampml)
## 
## Call:
## lm(formula = FT.DELTA.ampml ~ Zsp + Calf.C, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5993  -2.6501  -0.1213   2.9498  12.4158 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -18.003950   5.006479  -3.596 0.000487 ***
## Zsp           0.013488   0.006943   1.943 0.054625 .  
## Calf.C        0.371324   0.159115   2.334 0.021446 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.592 on 109 degrees of freedom
## Multiple R-squared:  0.1438, Adjusted R-squared:  0.1281 
## F-statistic: 9.154 on 2 and 109 DF,  p-value: 0.0002113

The model shows that both Zsp and Calf.C significantly affect FT.DELTA.ampml.

The second is FT.DELTA.vel regressed on Zsp only
# Fit the regression model for FT.DELTA.vel using only Zsp
model_FT_DELTA_vel <- lm(FT.DELTA.vel ~ Zsp, data = dat)
summary(model_FT_DELTA_vel)
## 
## Call:
## lm(formula = FT.DELTA.vel ~ Zsp, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1895 -1.4107 -0.2076  1.2214  7.0685 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.198623   1.535675  -0.129    0.897
## Zsp          0.004753   0.002990   1.589    0.115
## 
## Residual standard error: 2.283 on 110 degrees of freedom
## Multiple R-squared:  0.02245,    Adjusted R-squared:  0.01356 
## F-statistic: 2.526 on 1 and 110 DF,  p-value: 0.1148

The model shows that Zsp significantly affects FT.DELTA.vel.