Paired t-tests for LCHF study

Repeated measures low carb high fat diet

Packages

Load the various packages.

# packages for analysis of repeated measures
library(nlme)
library(car)
# document definition and construction
library(knitr)
opts_knit$set(eval.after = "fig.cap")
library(rmarkdown)
#reading and writing files of various types (in this case, csv and text files)
library(readr)
# data manipulation
library(tidyr)
library(dplyr)
library(tibble)
# pipelining commands
library(magrittr)
# plots
library(ggplot2)

Prepare data

obs <- read_csv("mep_observations.csv")

Parsed with column specification:
cols(
  Subject = col_integer(),
  TestNumber = col_integer(),
  MEPHR = col_double(),
  Weight = col_double(),
  MaxRelVO2 = col_double(),
  FatBIA = col_double(),
  FatSF = col_double(),
  Avgkcals = col_integer(),
  AvgFat = col_integer(),
  AvgCHO = col_integer(),
  AvgPro = col_integer()
)

# these models work with factors
obs$Subject <- as.factor(obs$Subject)
obs$TestNumber <- as.factor(obs$TestNumber)
obs <- as_tibble(obs)
obs$Week <- obs$TestNumber
levels(obs$Week) <- c("Baseline", "Week1", "Week2", "Week3")
obs <- obs[,c("Subject", "TestNumber", "Week",   # move Week next to TestNumber
              "MEPHR", "Weight", "MaxRelVO2",    # for easy validation
              "FatBIA", "FatSF", 
              "Avgkcals", "AvgFat", "AvgCHO", "AvgPro")]
obs

---

Average KCal

Linear mixed-effects model

lme.model <- lme(Avgkcals ~ Week, 
                 random = ~1|Subject/Week, 
                 data=obs)
summary(lme.model)

Linear mixed-effects model fit by REML
 Data: obs 
       AIC      BIC    logLik
  266.3145 271.7226 -126.1573

Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    498.0948

 Formula: ~1 | Week %in% Subject
        (Intercept) Residual
StdDev:    376.1294 172.2395

Fixed effects: Avgkcals ~ Week 
             Value Std.Error DF   t-value p-value
(Intercept) 2148.4  289.5646 12  7.419416  0.0000
WeekWeek1   -171.0  261.6408 12 -0.653568  0.5257
WeekWeek2   -228.0  261.6408 12 -0.871424  0.4006
WeekWeek3     88.6  261.6408 12  0.338632  0.7407
 Correlation: 
          (Intr) WekWk1 WekWk2
WeekWeek1 -0.452              
WeekWeek2 -0.452  0.500       
WeekWeek3 -0.452  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.74587052 -0.20976152 -0.03953116  0.23436058  0.76175514 

Number of Observations: 20
Number of Groups: 
          Subject Week %in% Subject 
                5                20 

# print F- and p-values
anova(lme.model)

            numDF denDF  F-value p-value
(Intercept)     1    12 73.71017  <.0001
Week            3    12  0.63289  0.6078

Boxplot

aboxplot <- ggplot(data=obs, aes(x=Week,y=Avgkcals)) + 
            geom_boxplot() +
            labs(x="", y=expression("Mean kcal/day"))
aboxplot

Pairwise t-test

pairwise.t.test(x = obs$Avgkcals,
                g = obs$Week,
                p.adjust.method = "bonferroni",
                pool.sd = FALSE,
                paired = TRUE )

    Pairwise comparisons using paired t tests 

data:  obs$Avgkcals and obs$Week 

      Baseline Week1 Week2
Week1 1        -     -    
Week2 1        1     -    
Week3 1        1     1    

P value adjustment method: bonferroni 

Multivariate approach using lm() / mlm()

1. Convert data frame to matrix. Each row of the matrix is a subject. Each column is the test number.

obs.matrix <- with(obs, # data from obs data frame
                    cbind(
                      Avgkcals[Week=="Baseline"], # subject observations baseline
                      Avgkcals[Week=="Week1"], # subject observations test 1
                      Avgkcals[Week=="Week2"], # etc
                      Avgkcals[Week=="Week3"]
                      ) # cbind column bind, rows are flipped to columns
                  )
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$Week)
# display the matrix
obs.matrix

       Baseline Week1 Week2 Week3
112121     3330  2723  1604  2618
258996     1489  1512  1795  1468
456121     1748  1726  1754  1759
522210     2555  2253  2959  3180
563751     1620  1673  1490  2160

2. Create a multivariate linear model from the matrix.

mlm.model <- lm(obs.matrix ~ 1)
mlm.model

Call:
lm(formula = obs.matrix ~ 1)

Coefficients:
             Baseline  Week1  Week2  Week3
(Intercept)  2148      1977   1920   2237 

The coefficients are equal to the means across subjects at each test.

3. Create a set of factors.

rfactor <- as.factor(levels(obs$Week))
rfactor

[1] Baseline Week1    Week2    Week3   
Levels: Baseline Week1 Week2 Week3

4. Define a new anova model object mlm.model.aov.

mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
                       idata = data.frame(rfactor),
                       idesign = ~rfactor, 
                       type = "III")
summary(mlm.model.aov, 
        multivariate = FALSE) # don't show multivariate tests

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

                  SS num Df Error SS den Df       F   Pr(>F)   
(Intercept) 85764253      1  4654134      4 73.7102 0.001011 **
rfactor       324940      3  2053677     12  0.6329 0.607835   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

        Test statistic p-value
rfactor        0.12623 0.37095


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

         GG eps Pr(>F[GG])
rfactor 0.49116     0.5177

           HF eps Pr(>F[HF])
rfactor 0.7081437  0.5633187

---

Average Fat

Linear mixed-effects model

lme.model <- lme(AvgFat ~ Week, 
                 random = ~1|Subject/Week, 
                 data=obs)
summary(lme.model)

Linear mixed-effects model fit by REML
 Data: obs 
       AIC      BIC    logLik
  178.0157 183.4238 -82.00785

Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    21.31392

 Formula: ~1 | Week %in% Subject
        (Intercept) Residual
StdDev:    26.64973 10.91597

Fixed effects: AvgFat ~ Week 
            Value Std.Error DF  t-value p-value
(Intercept)  82.6  16.02280 12 5.155155  0.0002
WeekWeek1    40.4  18.21391 12 2.218084  0.0466
WeekWeek2    38.8  18.21391 12 2.130240  0.0545
WeekWeek3    58.4  18.21391 12 3.206340  0.0075
 Correlation: 
          (Intr) WekWk1 WekWk2
WeekWeek1 -0.568              
WeekWeek2 -0.568  0.500       
WeekWeek3 -0.568  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.69372218 -0.19205783 -0.01038382  0.08039667  0.60666552 

Number of Observations: 20
Number of Groups: 
          Subject Week %in% Subject 
                5                20 

# print F- and p-values
anova(lme.model)

            numDF denDF  F-value p-value
(Intercept)     1    12 103.4498  <.0001
Week            3    12   3.6468  0.0446

Boxplot

aboxplot <- ggplot(data=obs, aes(x=Week,y=AvgFat)) + 
            geom_boxplot() +
            annotate(geom="text", label="*", x=2.01, y=180, size=6, family="serif") +
            labs(x="", y=expression("Mean Fat g/day"))
aboxplot

Pairwise t-test

pairwise.t.test(x = obs$AvgFat,
                g = obs$Week,
                p.adjust.method = "bonferroni",
                pool.sd = FALSE,
                paired = TRUE )

    Pairwise comparisons using paired t tests 

data:  obs$AvgFat and obs$Week 

      Baseline Week1 Week2
Week1 0.031    -     -    
Week2 1.000    1.000 -    
Week3 0.165    1.000 1.000

P value adjustment method: bonferroni 

Multivariate approach using lm() / mlm()

1. Convert data frame to matrix. Each row of the matrix is a subject. Each column is the test number.

obs.matrix <- with(obs, # data from obs data frame
                    cbind(
                      AvgFat[Week=="Baseline"], # subject observations baseline
                      AvgFat[Week=="Week1"], # subject observations test 1
                      AvgFat[Week=="Week2"], # etc
                      AvgFat[Week=="Week3"]
                      ) # cbind column bind, rows are flipped to columns
                  )
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$Week)
# display the matrix
obs.matrix

       Baseline Week1 Week2 Week3
112121      145   172    85   161
258996       49    89   112   103
456121       68   117   117   113
522210       75   138   181   196
563751       76    99   112   132

2. Create a multivariate linear model from the matrix.

mlm.model <- lm(obs.matrix ~ 1)
mlm.model

Call:
lm(formula = obs.matrix ~ 1)

Coefficients:
             Baseline  Week1  Week2  Week3
(Intercept)   82.6     123.0  121.4  141.0

The coefficients are equal to the means across subjects at each test.

3. Create a set of factors.

rfactor <- as.factor(levels(obs$Week))
rfactor

[1] Baseline Week1    Week2    Week3   
Levels: Baseline Week1 Week2 Week3

4. Define a new anova model object mlm.model.aov.

mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
                       idata = data.frame(rfactor),
                       idesign = ~rfactor, 
                       type = "III")
summary(mlm.model.aov, 
        multivariate = FALSE) # don't show multivariate tests

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

                SS num Df Error SS den Df        F    Pr(>F)    
(Intercept) 273780      1  10586.0      4 103.4498 0.0005263 ***
rfactor       9074      3   9952.4     12   3.6468 0.0445511 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

        Test statistic p-value
rfactor       0.063552 0.21124


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

         GG eps Pr(>F[GG])
rfactor 0.43137     0.1096

         HF eps Pr(>F[HF])
rfactor 0.55072 0.09017872

---

Average Carbohydrate

Linear mixed-effects model

lme.model <- lme(AvgCHO ~ Week, 
                 random = ~1|Subject/Week, 
                 data=obs)
summary(lme.model)

Linear mixed-effects model fit by REML
 Data: obs 
       AIC      BIC    logLik
  183.5791 188.9872 -84.78954

Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    34.10852

 Formula: ~1 | Week %in% Subject
        (Intercept) Residual
StdDev:    29.22064 12.98414

Fixed effects: AvgCHO ~ Week 
             Value Std.Error DF   t-value p-value
(Intercept)  249.0  20.90849 12 11.909038   0e+00
WeekWeek1   -116.8  20.22309 12 -5.775576   1e-04
WeekWeek2   -136.6  20.22309 12 -6.754655   0e+00
WeekWeek3   -103.6  20.22309 12 -5.122857   3e-04
 Correlation: 
          (Intr) WekWk1 WekWk2
WeekWeek1 -0.484              
WeekWeek2 -0.484  0.500       
WeekWeek3 -0.484  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.59192367 -0.37207554  0.06267628  0.19456751  0.61974464 

Number of Observations: 20
Number of Groups: 
          Subject Week %in% Subject 
                5                20 

# print F- and p-values
anova(lme.model)

            numDF denDF  F-value p-value
(Intercept)     1    12 89.92273  <.0001
Week            3    12 18.21229   1e-04

Boxplot

aboxplot <- ggplot(data=obs, aes(x=Week,y=AvgCHO)) + 
            geom_boxplot() +
            annotate(geom="text", label="*", x=2.01, y=190, size=6, family="serif") +
            labs(x="", y=expression("Mean  Carbohydrates g/day"))
aboxplot

Pairwise t-test

pairwise.t.test(x = obs$AvgCHO,
                g = obs$Week,
                p.adjust.method = "bonferroni",
                pool.sd = FALSE,
                paired = TRUE )

    Pairwise comparisons using paired t tests 

data:  obs$AvgCHO and obs$Week 

      Baseline Week1 Week2
Week1 0.044    -     -    
Week2 0.052    1.000 -    
Week3 0.024    1.000 0.732

P value adjustment method: bonferroni 

Multivariate approach using lm() / mlm()

1. Convert data frame to matrix. Each row of the matrix is a subject. Each column is the test number.

obs.matrix <- with(obs, # data from obs data frame
                    cbind(
                      AvgCHO[Week=="Baseline"], # subject observations baseline
                      AvgCHO[Week=="Week1"], # subject observations test 1
                      AvgCHO[Week=="Week2"], # etc
                      AvgCHO[Week=="Week3"]
                      ) # cbind column bind, rows are flipped to columns
                  )
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$Week)
# display the matrix
obs.matrix

       Baseline Week1 Week2 Week3
112121      319   176    91   176
258996      195   122   110    91
456121      196   106   110   130
522210      339   143   161   197
563751      196   114    90   133

2. Create a multivariate linear model from the matrix.

mlm.model <- lm(obs.matrix ~ 1)
mlm.model

Call:
lm(formula = obs.matrix ~ 1)

Coefficients:
             Baseline  Week1  Week2  Week3
(Intercept)  249.0     132.2  112.4  145.4

The coefficients are equal to the means across subjects at each test.

3. Create a set of factors.

rfactor <- as.factor(levels(obs$Week))
rfactor

[1] Baseline Week1    Week2    Week3   
Levels: Baseline Week1 Week2 Week3

4. Define a new anova model object mlm.model.aov.

mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
                       idata = data.frame(rfactor),
                       idesign = ~rfactor, 
                       type = "III")
summary(mlm.model.aov, 
        multivariate = FALSE) # don't show multivariate tests

HF eps > 1 treated as 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

                SS num Df Error SS den Df      F    Pr(>F)    
(Intercept) 510401      1    22704      4 89.923 0.0006901 ***
rfactor      55863      3    12269     12 18.212 9.196e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

        Test statistic p-value
rfactor        0.37407 0.76403


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

         GG eps Pr(>F[GG])   
rfactor 0.62838   0.001396 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          HF eps   Pr(>F[HF])
rfactor 1.170419 9.196289e-05

---

Average Protein

Linear mixed-effects model

lme.model <- lme(AvgPro ~ Week, 
                 random = ~1|Subject/Week, 
                 data=obs)
summary(lme.model)

Linear mixed-effects model fit by REML
 Data: obs 
       AIC     BIC    logLik
  169.7359 175.144 -77.86793

Random effects:
 Formula: ~1 | Subject
        (Intercept)
StdDev:    21.94177

 Formula: ~1 | Week %in% Subject
        (Intercept) Residual
StdDev:    19.00878 8.424649

Fixed effects: AvgPro ~ Week 
            Value Std.Error DF   t-value p-value
(Intercept) 101.8  13.51850 12  7.530419  0.0000
WeekWeek1   -18.4  13.15003 12 -1.399236  0.1871
WeekWeek2    -8.4  13.15003 12 -0.638782  0.5350
WeekWeek3     0.0  13.15003 12  0.000000  1.0000
 Correlation: 
          (Intr) WekWk1 WekWk2
WeekWeek1 -0.486              
WeekWeek2 -0.486  0.500       
WeekWeek3 -0.486  0.500  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.68323582 -0.21701377  0.05656356  0.26974876  0.43819577 

Number of Observations: 20
Number of Groups: 
          Subject Week %in% Subject 
                5                20 

# print F- and p-values
anova(lme.model)

            numDF denDF  F-value p-value
(Intercept)     1    12 76.70675  <.0001
Week            3    12  0.88502  0.4765

Boxplot

aboxplot <- ggplot(data=obs, aes(x=Week,y=AvgPro)) + 
            geom_boxplot() +
            labs(x="", y=expression("Mean Protein g/day"))
aboxplot

Pairwise t-test

pairwise.t.test(x = obs$AvgPro,
                g = obs$Week,
                p.adjust.method = "bonferroni",
                pool.sd = FALSE,
                paired = TRUE )

    Pairwise comparisons using paired t tests 

data:  obs$AvgPro and obs$Week 

      Baseline Week1 Week2
Week1 0.53     -     -    
Week2 1.00     1.00  -    
Week3 1.00     1.00  1.00 

P value adjustment method: bonferroni 

Multivariate approach using lm() / mlm()

1. Convert data frame to matrix. Each row of the matrix is a subject. Each column is the test number.

obs.matrix <- with(obs, # data from obs data frame
                    cbind(
                      AvgPro[Week=="Baseline"], # subject observations baseline
                      AvgPro[Week=="Week1"], # subject observations test 1
                      AvgPro[Week=="Week2"], # etc
                      AvgPro[Week=="Week3"]
                      ) # cbind column bind, rows are flipped to columns
                  )
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$Week)
# display the matrix
obs.matrix

       Baseline Week1 Week2 Week3
112121      138   120    75   129
258996       74    54    96    59
456121       85    72    79    67
522210      132    86   138   146
563751       80    85    79   108

2. Create a multivariate linear model from the matrix.

mlm.model <- lm(obs.matrix ~ 1)
mlm.model

Call:
lm(formula = obs.matrix ~ 1)

Coefficients:
             Baseline  Week1  Week2  Week3
(Intercept)  101.8      83.4   93.4  101.8

The coefficients are equal to the means across subjects at each test.

3. Create a set of factors.

rfactor <- as.factor(levels(obs$Week))
rfactor

[1] Baseline Week1    Week2    Week3   
Levels: Baseline Week1 Week2 Week3

4. Define a new anova model object mlm.model.aov.

mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
                       idata = data.frame(rfactor),
                       idesign = ~rfactor, 
                       type = "III")
summary(mlm.model.aov, 
        multivariate = FALSE) # don't show multivariate tests

HF eps > 1 treated as 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

                SS num Df Error SS den Df      F    Pr(>F)    
(Intercept) 180880      1   9432.3      4 76.707 0.0009368 ***
rfactor       1148      3   5187.7     12  0.885 0.4764684    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

        Test statistic p-value
rfactor        0.40954  0.7994


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

        GG eps Pr(>F[GG])
rfactor 0.6553     0.4484

          HF eps Pr(>F[HF])
rfactor 1.283041  0.4764684