Objective
Why is this important?
Set up
Q1. Explore energy intake patterns
- Q1.A. Calculate the mean total energy intake and its standard deviation (SD)
- Q1.B. Calculate the percentage of energy intake from morning meals, snacks, lunch and dinner:
Q2. Home vs. out-of-home eating habits
- Compare the percentage of energy intake from food prepared at home versus food eaten out of home.
Q3. Energy intake in men vs women
Q4. Fiber intake and health marker

The dataset (“24h_recall”) is derived from weighed 24-hour food intake data recording of students. The data is available in different formats (eg. xls, csv) that allow you to import it in different statistical software packages such as R, Stata, SPSS or SPlus. For this exercise, and the solutions, R will be used?

Datasets are often described using ‘dictionaries’, which give an easy to excess and understandable description of all variables present in the dataset.

Objective

The goal of this exercise is to provide hands-on experience in analyzing dietary data and applying energy adjustment methods.

Why is this important?

When analyzing dietary data, individual differences in energy intake can bias relationships between nutrients and health outcomes. For example, people who consume more energy overall may also eat more fruit and fiber. To accurately assess the association between fruit intake and fiber intake, it’s crucial to adjust for total energy intake.

Set up

Get required packages

Before starting your analysis in R, you need to load the necessary packages. If you haven’t installed these packages before, run install.packages() first.

Tip: Install each package only once. After installation, use library() to load them whenever needed.

Install packages

# Run this only once to install the required packages
# install.packages("tidyverse")  # For data science workflows
# install.packages("pastecs")    # For descriptive statistics
# install.packages("modelr")     # For regression models
# install.packages("broom")      # For tidying model outputs
# install.packages("ggplot2")    # For data visualization
# install.packages("readr")      # For reading data files
# install.packages("purrr")      # For functional programming

Load packages

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# Tidyverse: A suite of data science packages including ggplot2, dplyr, tidyr, readr, purrr, tibble, and more. 
# It simplifies tasks like data cleaning, transformation, and visualization.

library(pastecs)

## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

# Pastecs: Provides tools for descriptive statistics, e.g., 'stat.desc()' for means, medians, SDs, and skewness.

library(modelr) 
# Modelr: Works seamlessly with tidyverse to streamline regression modeling workflows.
# Functions like 'add_predictions()' and 'add_residuals()' help integrate model outputs with data.

library(broom)

## 
## Attaching package: 'broom'
## 
## The following object is masked from 'package:modelr':
## 
##     bootstrap

# Broom: Converts statistical model outputs into tidy data frames, making them easier to analyze and visualize.
# Example: 'tidy()' summarizes model estimates, 'augment()' adds predictions, 'glance()' gives overall model stats.

library(ggplot2) 
# GGplot2: A data visualization package to create high-quality, customizable graphs using a layered grammar of graphics.

library(readr) 
# Readr: Provides fast and friendly functions for importing and reading data, e.g., 'read_csv()' for CSV files.

library(purrr) 
# Purrr: Simplifies functional programming tasks, like applying functions over lists and vectors with 'map()'.
# It works well with tidyverse pipelines for iterative operations.

Download database

First, download the file “24h_recall.csv” from Ufora, and save it to a location on your computer where you can easily find it.

Option 1: Set Your Working Directory

To load the file into R, you need to set the folder where you saved it as your working directory. For example, if you saved the file in:

You can set the working directory with the following command:

setwd("/Users/jeroenberden/Library/CloudStorage/OneDrive-UGent/PhD/module5_files")

Option 2: Import Manually

Alternatively, you can import the file directly into RStudio without setting the working directory:

Go to the Environment tab.
Click Import Dataset → From Text (readr).
Browse to the file location and follow the prompts.

Import your dataset

We will use the first option described above. First, set your working directory to the folder where you saved the file. Then, use read_csv() to load the data into R.

setwd("/Users/jeroenberden/Library/CloudStorage/OneDrive-UGent/PhD/module5_files")
energy <- read_csv("24h_recall.csv")

## Rows: 68 Columns: 30
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): gender
## dbl (29): id, quantity, fruit_weight, vegetable_weight, protein_gram, fat_gr...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

energy <- energy %>%
  mutate_all(~replace(., is.na(.), 0))
# 'read_csv()' is a function from the 'readr' package, part of the 'tidyverse'.
# It reads CSV (comma-separated values) files into R and creates a clean "tibble" data frame.
# The dataset "database energy adjustment.csv" is now stored in the 'energy' object 
# for further analysis and manipulation.

Q1. Explore energy intake patterns

Q1.A. Calculate the mean total energy intake and its standard deviation (SD)

Before calculating the mean and SD, it’s good practice to explore the data visually for potential outliers. Outliers can be removed for a cleaner analysis. Here, we visualize the distribution using histograms and boxplots and then filter out extreme outliers.

Remove outliers and save your new dataset

Visualise current dataset for energy intake to visually get an idea of potential outliers.

# Histogram: Energy intake distribution
ggplot(energy, aes(x = energy_kcal)) +
  geom_histogram(binwidth = 100, fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Total Energy Intake (kcal/day)",
       x = "Energy Intake (kcal/day)",
       y = "Frequency") +
  theme_minimal()

# Boxplot: Detect outliers in energy intake
ggplot(energy, aes(y = energy_kcal)) +
  geom_boxplot(fill = "lightgreen", color = "black", alpha = 0.7) +
  labs(title = "Boxplot of Total Energy Intake (kcal/day)",
       y = "Energy Intake (kcal/day)") +
  theme_minimal()

Remove extreme outliers and save as a new dataset. There are different ways to do this, but this exercise opts for a energy intake cut-off value of 4000 kcalories per day.

energy_filtered <- filter(energy, energy_kcal < 4000)

Option 1: Calculate Mean and SD Individually

Visualise the filtered data.

# Histogram:
ggplot(energy_filtered, aes(x = energy_kcal)) +
  geom_histogram(binwidth = 100, fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Total Energy Intake (kcal/day) - After Outlier Removal",
       x = "Energy Intake (kcal/day)",
       y = "Frequency") +
  theme_minimal()

# Boxplot
ggplot(energy_filtered, aes(y = energy_kcal)) +
  geom_boxplot(fill = "lightgreen", color = "black", alpha = 0.7) +
  labs(title = "Boxplot of Total Energy Intake (kcal/day) - After Outlier Removal",
       y = "Energy Intake (kcal/day)") +
  theme_minimal()

Calculate mean and SD individually.

mean(energy_filtered$energy_kcal)

## [1] 2030.047

sd(energy_filtered$energy_kcal)

## [1] 686.8057

Option 2: Calculate Summary Statistics Together

Here, we calculate all relevant summary statistics (mean, SD, min, max, and count) and save them in a single object called sum_energy.

sum_energy <- energy_filtered %>%
  summarise(
    count = n(),                              # Number of observations
    mean_energy = mean(energy_kcal, na.rm = TRUE),  # Mean energy intake
    sd_energy = sd(energy_kcal, na.rm = TRUE),      # Standard deviation
    min_energy = min(energy_kcal, na.rm = TRUE),    # Minimum value
    max_energy = max(energy_kcal, na.rm = TRUE)     # Maximum value
  )

sum_energy

## # A tibble: 1 × 5
##   count mean_energy sd_energy min_energy max_energy
##   <int>       <dbl>     <dbl>      <dbl>      <dbl>
## 1    62       2030.      687.       739.      3982.

Results

After the removal of 6 extreme values (using an arbitrary cut-off of > 4000 kcal for this exercise), the dataset now contains 62 observations.

Descriptive Statistics:

Mean energy intake: 2030 kcal
Standard deviation (SD): 687
Minimum value: 739 kcal
Maximum value: 3982 kcal
Sample size (n): 62

Q1.B. Calculate the percentage of energy intake from morning meals, snacks, lunch and dinner:

Understand Your Data

Before we begin analyzing, it’s important to check the basic summary statistics of our data. This helps us understand the spread and distribution of the values. We’ll calculate the minimum, maximum, mean, and other key measures to get a feel for our data. Checking for any unusual values or outliers is always a good practice.

First, we create a smaller dataset that contains the energy intake for each time period: morning meals, lunch, dinner, and snacks. Then, we’ll calculate the descriptive statistics for these variables.

Create a new dataset with energy intake by time of day.

energy_stat_time <- cbind(
  energy_filtered$energy_morning, 
  energy_filtered$energy_lunch, 
  energy_filtered$energy_dinner, 
  energy_filtered$energy_snack
)

energy_stat_time <- as.data.frame(energy_stat_time)
stat.desc(energy_stat_time)

##                        V1           V2           V3           V4
## nbr.val         62.000000    62.000000    62.000000    62.000000
## nbr.null        43.000000    42.000000    42.000000    43.000000
## nbr.na           0.000000     0.000000     0.000000     0.000000
## min              0.000000     0.000000     0.000000     0.000000
## max           1222.132476  1129.831360  1188.460577  1325.499718
## range         1222.132476  1129.831360  1188.460577  1325.499718
## sum          10104.561574 11753.289660 10980.760507  7813.429185
## median           0.000000     0.000000     0.000000     0.000000
## mean           162.976800   189.569188   177.109040   126.023051
## SE.mean         38.369326    39.527638    37.903395    34.007897
## CI.mean.0.95    76.724210    79.040397    75.792522    68.002993
## var          91276.722273 96870.918151 89073.375051 71705.297258
## std.dev        302.120377   311.240933   298.451629   267.778448
## coef.var         1.853763     1.641833     1.685129     2.124837

Generating New Variables: Share of Energy Intake

Next, we want to calculate the percentage of total energy intake that comes from each time period (morning meals, lunch, dinner, and snacks). To do this, we divide the energy intake for each time period (e.g., energy_morning) by the total energy intake for the day (energy_kcal).

energy_share <- energy_filtered %>%
  mutate(
    share_morning = energy_morning / energy_kcal,
    share_lunch = energy_lunch / energy_kcal,
    share_dinner = energy_dinner / energy_kcal,
    share_snack = energy_snack / energy_kcal,
    
    total_share = share_morning + share_lunch + share_dinner + share_snack
  ) %>%
  filter(total_share != 0)

Calculate the Mean Percentage for Each Time Period

Once we have the percentages, we can calculate the average (mean) percentage of total energy intake for each time period (morning meals, lunch, dinner, snacks). This helps us understand how energy is distributed across the day.

means <- energy_share %>%
  summarise(
    mean_morning = mean(share_morning, na.rm = TRUE)*100,
    mean_lunch = mean(share_lunch, na.rm = TRUE)*100,
    mean_dinner = mean(share_dinner, na.rm = TRUE)*100,
    mean_snack = mean(share_snack, na.rm = TRUE)*100
  )

print(means)

## # A tibble: 1 × 4
##   mean_morning mean_lunch mean_dinner mean_snack
##          <dbl>      <dbl>       <dbl>      <dbl>
## 1         23.3       30.1        28.2       18.4

Visualize the data

To get a better understanding of how energy intake is distributed across the day, we’ll visualize the data using a boxplot. This will show the spread of the data for each time period, and we’ll add points for the average (mean) percentage.

energy_long <- energy_share %>%
  select(share_morning, share_lunch, share_dinner, share_snack) %>%
  pivot_longer(
    cols = everything(),
    names_to = "time_period",
    values_to = "share"
  )

# Create a boxplot
ggplot(energy_long, aes(x = time_period, y = share, fill = time_period)) +
  geom_boxplot(outlier.colour = "red", outlier.size = 1.5) +
  stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "blue") +  # Add mean points
  scale_y_continuous(labels = scales::percent) +
  labs(
    title = "Distribution of Energy Intake by Time Period",
    x = "Time Period",
    y = "Percentage of Total Energy Intake",
    fill = "Time Period"
  ) +
  theme_minimal()

Result

Finally, after calculating the percentages and visualizing the data, you will see the distribution of energy intake across the time periods. Based on the analysis, the typical distribution of energy intake might look something like this:

Morning meals: 23% of total energy intake
Lunch: 30% of total energy intake
Dinner: 28% of total energy intake
Snacks: 18% of total energy intake

Q2. Home vs. out-of-home eating habits

Compare the percentage of energy intake from food prepared at home versus food eaten out of home.

We will compare the percentage of energy intake from food prepared at home versus food eaten out of home. By combining all out-of-home energy sources into one variable, we can easily analyze the data.

Combine Out-of-home eating habits

We first create a new variable called energy_allout that sums up energy intake from UGent, other places, and unspecified sources. This will represent all food consumed out of home.

energy_out_of_home <- energy_filtered %>%
  mutate(
    energy_allout = energy_ugent + energy_out + energy_no
  )

Descriptive Statistics for Home vs. Out-of-Home Energy Intake

To understand the distribution of energy intake for both home-prepared food (energy_home) and out-of-home food (energy_allout), we calculate descriptive statistics using the stat.desc() function.

stat.desc(energy_out_of_home[c("energy_home", "energy_allout")])

##               energy_home energy_allout
## nbr.val      6.200000e+01  6.200000e+01
## nbr.null     0.000000e+00  1.700000e+01
## nbr.na       0.000000e+00  0.000000e+00
## min          1.142985e+02  0.000000e+00
## max          3.982048e+03  2.392214e+03
## range        3.867750e+03  2.392214e+03
## sum          9.392621e+04  3.193671e+04
## median       1.300452e+03  2.884729e+02
## mean         1.514939e+03  5.151083e+02
## SE.mean      1.116469e+02  6.836276e+01
## CI.mean.0.95 2.232517e+02  1.366998e+02
## var          7.728316e+05  2.897550e+05
## std.dev      8.791084e+02  5.382889e+02
## coef.var     5.802930e-01  1.045002e+00

Visualise energy intake (home vs out of home)

Next, we reshape the data into a long format for visualization and create a boxplot to compare energy intake between the two locations.

energy_long <- energy_out_of_home %>%
  select(energy_home, energy_allout) %>%
  pivot_longer(cols = everything(), names_to = "location", values_to = "energy_intake")

ggplot(energy_long, aes(x = location, y = energy_intake, fill = location)) +
  geom_boxplot() +
  stat_summary(fun = mean, geom = "point", shape = 20, size = 3, color = "blue") +  # Add mean points
  labs(
    title = "Comparison of Energy Intake: Home vs. Out of Home",
    x = "Location",
    y = "Energy Intake (kcal)",
    fill = "Location"
  ) +
  theme_minimal()

Calculate the percentage of energy intake

We now calculate the percentage of total energy intake (energy_kcal) coming from:

Home-prepared food (share_home)
Out-of-home food (share_allout)

energy_out_of_home <- energy_out_of_home %>%
  mutate(
    share_home = energy_home / energy_kcal,     # Percentage of energy intake from home
    share_allout = energy_allout / energy_kcal  # Percentage of energy intake from out-of-home sources
  )

Calculate the mean percentages

To summarize the results, we calculate the mean percentage of energy intake for home-prepared and out-of-home food.

energy_out_of_home_shares <- energy_out_of_home %>%
  summarise(
    mean_home = mean(share_home, na.rm = TRUE) * 100,     # Convert to percentage
    mean_allout = mean(share_allout, na.rm = TRUE) * 100  # Convert to percentage
  )

print(energy_out_of_home_shares)

## # A tibble: 1 × 2
##   mean_home mean_allout
##       <dbl>       <dbl>
## 1      72.5        27.5

Visualise the mean percentages

To make the results more intuitive, we create a bar plot to show the average percentage of energy intake from home and out-of-home sources.

energy_shares_long <- pivot_longer(energy_out_of_home_shares, cols = everything(), names_to = "source", values_to = "percentage")

# Bar plot
ggplot(energy_shares_long, aes(x = source, y = percentage, fill = source)) +
  geom_bar(stat = "identity", width = 0.6) +
  labs(
    title = "Percentage of Energy Intake: Home vs. Out of Home",
    x = "Source of Energy Intake",
    y = "Mean Percentage (%)",
    fill = "Source"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c("mean_home" = "skyblue", "mean_allout" = "salmon")) +
  geom_text(aes(label = round(percentage, 1)), vjust = -0.5)

Statistical Test for Normality

The Shapiro-Wilk test is used to statistically test if the data is normally distributed:

Null hypothesis: The data is normally distributed.
Alternative hypothesis: The data is not normally distributed.
If the p-value < 0.05, we reject the null hypothesis, indicating the data is not normally distributed.

shapiro.test(energy_out_of_home$energy_kcal)

## 
##  Shapiro-Wilk normality test
## 
## data:  energy_out_of_home$energy_kcal
## W = 0.94405, p-value = 0.006966

shapiro.test(energy_out_of_home$share_home)

## 
##  Shapiro-Wilk normality test
## 
## data:  energy_out_of_home$share_home
## W = 0.88293, p-value = 2.526e-05

shapiro.test(energy_out_of_home$share_allout)

## 
##  Shapiro-Wilk normality test
## 
## data:  energy_out_of_home$share_allout
## W = 0.88293, p-value = 2.526e-05

The distribution of energy_kcal, share_home and share_allout are normally distributed.

Paired t-test

t.test(energy_out_of_home$share_home, energy_out_of_home$share_allout, paired= TRUE)

## 
##  Paired t-test
## 
## data:  energy_out_of_home$share_home and energy_out_of_home$share_allout
## t = 6.6247, df = 61, p-value = 1.019e-08
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.3137455 0.5850360
## sample estimates:
## mean difference 
##       0.4493908

Result

•   Home-prepared food accounts for about 72.5% of total energy intake.
•   Out-of-home food accounts for about 27.5% of total energy intake.

Q3. Energy intake in men vs women

Visualize energy intake in men and women separately.

Make boxplots.

ggplot(energy_filtered, aes(x = gender, y = energy_kcal, fill = gender)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7) +
  labs(
    title = "Total Energy Intake by Gender",
    x = "Gender",
    y = "Total Energy Intake (kcal)"
  ) +
  theme_minimal()

T-test

shapiro_results <- energy_filtered %>%
  group_by(gender) %>%
  summarize(shapiro_p_value = shapiro.test(energy_kcal)$p.value)

print(shapiro_results)

## # A tibble: 2 × 2
##   gender shapiro_p_value
##   <chr>            <dbl>
## 1 female         0.490  
## 2 male           0.00512

As energy intake is normally distributed for women, but not for men, we will both perform a t-test and a non-parametric test to compare results.

t.test(energy_kcal ~ gender, data = energy_filtered)

## 
##  Welch Two Sample t-test
## 
## data:  energy_kcal by gender
## t = -0.42899, df = 49.989, p-value = 0.6698
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -428.1023  277.4176
## sample estimates:
## mean in group female   mean in group male 
##             1992.376             2067.718

wilcox.test(energy_kcal ~ gender, data = energy_filtered)

## 
##  Wilcoxon rank sum exact test
## 
## data:  energy_kcal by gender
## W = 499, p-value = 0.8014
## alternative hypothesis: true location shift is not equal to 0

Result

There is no signficant difference between men and women for energy intake, both the t-test and the non-paramteric alternative lead to the same conclusions.

Q4. Fiber intake and health marker

Q4.A. Calculate the mean fiber intake and its standard deviation (SD)

Visualise first.

# Histogram
ggplot(energy_filtered, aes(x = fibre_gram)) +
  geom_histogram(binwidth = 100, fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Total Fibre Intake (g/day)",  # Update title to reflect fibre intake in grams
       x = "Fibre Intake (g/day)",  # Update x-axis label
       y = "Frequency") +
  theme_minimal()

# Boxplot
ggplot(energy_filtered, aes(y = fibre_gram)) +
  geom_boxplot(fill = "lightgreen", color = "black", alpha = 0.7) +
  labs(title = "Boxplot of Total Fibre Intake (g/day)",  # Update title to reflect fibre intake in grams
       y = "Fibre Intake (g/day)") +  # Update y-axis label
  theme_minimal()

mean(energy_filtered$fibre_gram, na.rm = TRUE)

## [1] 20.74819

sd(energy_filtered$fibre_gram, na.rm = TRUE)

## [1] 13.02773

Result

Summary of fibre intake in grams per day:

mean = 21
SD = 13

Q4.B. Calculate the mean health marker and its standard deviation (SD) for men and women separately

# Histogram
ggplot(energy_filtered, aes(x = health_marker)) +
  geom_histogram(binwidth = 100, fill = "skyblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Health Marker",  # Updated title to reflect health marker
       x = "Health Marker",  # Updated x-axis label
       y = "Frequency") +
  theme_minimal()

# Boxplot
ggplot(energy_filtered, aes(y = health_marker)) +
  geom_boxplot(fill = "lightgreen", color = "black", alpha = 0.7) +
  labs(title = "Boxplot of Health Marker",  # Updated title to reflect health marker
       y = "Health Marker") +  # Updated y-axis label
  theme_minimal()

mean(energy_filtered$health_marker, na.rm = TRUE)

## [1] 2262.216

sd(energy_filtered$health_marker, na.rm = TRUE)

## [1] 915.6574

Result

Summary of the health marker:

mean = 2262
SD = 916

Q4.C. Assess the association between fiber intake and the health marker

Scatter plot for the relationship fiber intake and health marker

ggplot(data = energy_filtered, mapping = aes(y = health_marker, x = fibre_gram)) +
  geom_point() +  # Scatter plot
  geom_smooth(se = TRUE) +  # Add smoothing line with confidence interval
  labs(title = "Health Marker vs. Fibre Intake",
       x = "Fibre Intake (grams)",
       y = "Health Marker") +
  theme_minimal()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

The scatter plot seems to indicate an association between fibre intake and health.

Assess the relationships using the different energy adjustment methods

Use the slides to help you interpret results!

After each model, don’t forget to check:

•   Independence of errors: Residuals vs fitted plot.
•   Homoscedasticity: Residuals vs fitted plot.
•   Normality of errors: Histogram, Q-Q plot, or Shapiro-Wilk test.

Method 0. Unadjusted model

model0 <- lm(health_marker ~ fibre_gram, data = energy_filtered)
summary(model0)

## 
## Call:
## lm(formula = health_marker ~ fibre_gram, data = energy_filtered)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2095.06  -466.10   -45.23   736.62  1553.79 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1970.869    217.321   9.069 7.56e-13 ***
## fibre_gram    14.042      8.891   1.579     0.12    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 904.6 on 60 degrees of freedom
## Multiple R-squared:  0.03991,    Adjusted R-squared:  0.02391 
## F-statistic: 2.494 on 1 and 60 DF,  p-value: 0.1195

confint(model0)

##                   2.5 %     97.5 %
## (Intercept) 1536.161704 2405.57622
## fibre_gram    -3.742308   31.82641

ggplot(data = data.frame(fitted = fitted(model0), residuals = residuals(model0)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model0))

Method 1. Energy partition model

energy_filtered$Remaining_Energy <- energy_filtered$energy_kcal - energy_filtered$fibre_gram * 2
model1 <- lm(health_marker ~ fibre_gram + Remaining_Energy, data = energy_filtered)
summary(model1)

## 
## Call:
## lm(formula = health_marker ~ fibre_gram + Remaining_Energy, data = energy_filtered)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1757.0  -448.7  -107.2   559.0  2435.1 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3188.2479   329.1250   9.687 8.35e-14 ***
## fibre_gram         52.3671    11.4769   4.563 2.61e-05 ***
## Remaining_Energy   -1.0121     0.2241  -4.517 3.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 786.4 on 59 degrees of freedom
## Multiple R-squared:  0.2866, Adjusted R-squared:  0.2624 
## F-statistic: 11.85 on 2 and 59 DF,  p-value: 4.71e-05

ggplot(data = data.frame(fitted = fitted(model1), residuals = residuals(model1)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model1))

Method 2. Standard model

model2 <- lm(health_marker ~ fibre_gram + energy_kcal, data = energy_filtered)
summary(model2)

## 
## Call:
## lm(formula = health_marker ~ fibre_gram + energy_kcal, data = energy_filtered)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1757.0  -448.7  -107.2   559.0  2435.1 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3188.2479   329.1250   9.687 8.35e-14 ***
## fibre_gram    54.3912    11.8120   4.605 2.25e-05 ***
## energy_kcal   -1.0121     0.2241  -4.517 3.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 786.4 on 59 degrees of freedom
## Multiple R-squared:  0.2866, Adjusted R-squared:  0.2624 
## F-statistic: 11.85 on 2 and 59 DF,  p-value: 4.71e-05

ggplot(data = data.frame(fitted = fitted(model2), residuals = residuals(model2)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model2))

Method 3. Nutrient density model

energy_filtered$fiber_per_kcal <- energy_filtered$fibre_gram/energy_filtered$energy_kcal
model3 <- lm(health_marker ~ fiber_per_kcal, data = energy_filtered)
summary(model3)

## 
## Call:
## lm(formula = health_marker ~ fiber_per_kcal, data = energy_filtered)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1764.72  -481.11   -58.62   606.17  1678.37 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      1153.8      316.4   3.647 0.000557 ***
## fiber_per_kcal 111538.6    30005.6   3.717 0.000445 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 832.4 on 60 degrees of freedom
## Multiple R-squared:  0.1872, Adjusted R-squared:  0.1736 
## F-statistic: 13.82 on 1 and 60 DF,  p-value: 0.0004448

ggplot(data = data.frame(fitted = fitted(model3), residuals = residuals(model3)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model3))

Method 4. Residual model

model_res <- lm(fibre_gram ~ energy_kcal, data = energy_filtered)
summary(model_res)

## 
## Call:
## lm(formula = fibre_gram ~ energy_kcal, data = energy_filtered)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.3796  -4.9480   0.6486   4.5226  25.6133 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8.372366   3.430933  -2.440   0.0176 *  
## energy_kcal  0.014345   0.001602   8.953 1.19e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.595 on 60 degrees of freedom
## Multiple R-squared:  0.5719, Adjusted R-squared:  0.5648 
## F-statistic: 80.15 on 1 and 60 DF,  p-value: 1.185e-12

energy_filtered$residual_fiber <- residuals(model_res)
model4 <- lm(health_marker ~ residual_fiber, data = energy_filtered)
summary(model4)

## 
## Call:
## lm(formula = health_marker ~ residual_fiber, data = energy_filtered)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2150.2  -428.9  -123.3   693.1  2092.3 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2262.22     101.11  22.373  < 2e-16 ***
## residual_fiber    54.39      11.96   4.548 2.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 796.2 on 60 degrees of freedom
## Multiple R-squared:  0.2564, Adjusted R-squared:  0.244 
## F-statistic: 20.69 on 1 and 60 DF,  p-value: 2.684e-05

ggplot(data = data.frame(fitted = fitted(model4), residuals = residuals(model4)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model4))

Method 5. All-component model

model5 <- lm(health_marker ~ fibre_gram + vitc_mg + protein_gram + fat_gram + carbohydrates_gram + na_mg, data = energy_filtered)
summary(model5)

## 
## Call:
## lm(formula = health_marker ~ fibre_gram + vitc_mg + protein_gram + 
##     fat_gram + carbohydrates_gram + na_mg, data = energy_filtered)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8354 -3.4568 -0.1075  3.1763 12.3069 
## 
## Coefficients:
##                      Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)         4.278e+03  2.318e+00  1845.53   <2e-16 ***
## fibre_gram          3.496e+01  7.847e-02   445.49   <2e-16 ***
## vitc_mg             8.914e-01  5.183e-03   172.00   <2e-16 ***
## protein_gram        8.378e-01  3.040e-02    27.56   <2e-16 ***
## fat_gram           -1.444e+00  3.367e-02   -42.88   <2e-16 ***
## carbohydrates_gram -1.487e+00  8.046e-03  -184.83   <2e-16 ***
## na_mg              -8.997e-01  8.051e-04 -1117.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.972 on 55 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 3.447e+05 on 6 and 55 DF,  p-value: < 2.2e-16

ggplot(data = data.frame(fitted = fitted(model5), residuals = residuals(model5)), aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs Fitted (Homoscedasticity)",
       x = "Fitted values",
       y = "Residuals") +
  theme_minimal()

hist(residuals(model5))

Result

Estimate of relationship between fiber intake in grams per day and the health marker for each model:

Model 0: 14.04
Model 1: 52.37
Model 2: 54.39
Model 3: 11539
Model 4: 54.40
Model 5: 34.96

Module 5 Exercise Session 24-hour Recall Data

Jeroen Berden

Objective

Why is this important?

Set up

Get required packages

Install packages

Load packages

Download database

Import your dataset

Q1. Explore energy intake patterns

Q1.A. Calculate the mean total energy intake and its standard deviation (SD)

Remove outliers and save your new dataset

Option 1: Calculate Mean and SD Individually

Option 2: Calculate Summary Statistics Together

Results

Q1.B. Calculate the percentage of energy intake from morning meals, snacks, lunch and dinner:

Understand Your Data

Generating New Variables: Share of Energy Intake

Calculate the Mean Percentage for Each Time Period

Visualize the data

Result

Q2. Home vs. out-of-home eating habits

Compare the percentage of energy intake from food prepared at home versus food eaten out of home.

Combine Out-of-home eating habits

Descriptive Statistics for Home vs. Out-of-Home Energy Intake

Visualise energy intake (home vs out of home)

Calculate the percentage of energy intake

Calculate the mean percentages

Visualise the mean percentages

Statistical Test for Normality

Paired t-test

Result

Q3. Energy intake in men vs women

Visualize energy intake in men and women separately.

T-test

Result

Q4. Fiber intake and health marker

Q4.A. Calculate the mean fiber intake and its standard deviation (SD)

Result

Q4.B. Calculate the mean health marker and its standard deviation (SD) for men and women separately

Result

Q4.C. Assess the association between fiber intake and the health marker

Assess the relationships using the different energy adjustment methods

Method 0. Unadjusted model

Method 1. Energy partition model

Method 2. Standard model

Method 3. Nutrient density model

Method 4. Residual model

Method 5. All-component model

Result

Module 5
Exercise Session
24-hour Recall Data