Food Shortened dataset
Katherine Ford and Danilo Pereira
4/17/2021
Shortened DataSet
Write some text using this mark to see what happens! And this one.
# Load dataset
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(readr))
#food_shortened <- read_csv("C:/Dropbox/Laboratorio/Carla Prado/food-shortened.csv")
food_shortened <- read_csv("~/Dropbox/Laboratorio/Carla Prado/food-shortened.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double(),
## `ASSIGNED TO:` = col_character(),
## `Respondent ID` = col_character(),
## muscle_loss = col_character()
## )
## ℹ Use `spec()` for the full column specifications.# cleaning names with default naming
food <- tibble::as_tibble(food_shortened, .name_repair = janitor::make_clean_names)
colnames(food)## [1] "assigned_to"
## [2] "respondent_id"
## [3] "importance_nutr"
## [4] "cancer_type"
## [5] "age"
## [6] "sex"
## [7] "stage"
## [8] "wt_chg_6mo_amt_lbs"
## [9] "muscle_loss"
## [10] "sga_appetite"
## [11] "sga_diarrhea"
## [12] "sga_taste"
## [13] "sga_nausea"
## [14] "sga_dry"
## [15] "sga_full"
## [16] "sga_vomit"
## [17] "sga_mouth"
## [18] "sga_tired"
## [19] "sga_constipation"
## [20] "sga_smell"
## [21] "sga_pain"
## [22] "sga_other"
## [23] "tx_chemo"
## [24] "tx_radiation"
## [25] "tx_surgery"
## [26] "tx_hormone"
## [27] "tx_other"
## [28] "tx_unsure"
## [29] "tx_surveillance"
## [30] "tx_none"
## [31] "marital"
## [32] "education"
## [33] "income"
## [34] "occupation"
## [35] "born_canada"
## [36] "ethnic"
## [37] "gram_weight_g"
## [38] "energy_kcal"
## [39] "sga_symmpt_score"
## [40] "energy_from_fat_percent_kcal"
## [41] "energy_from_carbohydrates_percent_kcal"
## [42] "energy_from_protein_percent_kcal"food$stage <- as.factor(food$stage)Frequentist analysis
Robust ANOVA
ANOVA with trimmed-means (t1way function https://rdrr.io/cran/WRS2/man/t1way.html) and Yuen’s trimmed means test https://rdrr.io/cran/WRS2/man/yuen.html for pairwise test. Set type = Robust (Explanatory measure of effect size). See R package information: https://indrajeetpatil.github.io/ggstatsplot/index.html
suppressPackageStartupMessages(library(ggstatsplot))
set.seed(123)
ggbetweenstats(
data = food,
x = muscle_loss,
y = energy_from_protein_percent_kcal,
type = "Robust",
plot.type = "box",
pairwise.comparisons = TRUE,
pairwise.display = "all",
xlab = "Self-perceived muscle loss",
ylab = "Protein intake (% kcal)",
ggsignif.args = list(textsize = 5, tip_length = 0.03),
centrality.label.args = list(size = 4, nudge_x = 0.4, segment.linetype = 4),
ggstatsplot.layer = FALSE
)
set.seed(123)
ggbetweenstats(
data = food,
x = muscle_loss,
y = energy_from_carbohydrates_percent_kcal,
type = "Robust",
plot.type = "box",
pairwise.comparisons = TRUE,
pairwise.display = "all",
xlab = "Self-perceived muscle loss",
ylab = "Carbohydrates intake (% kcal)",
ggsignif.args = list(textsize = 5, tip_length = 0.03),
centrality.label.args = list(size = 4, nudge_x = 0.4, segment.linetype = 4),
ggstatsplot.layer = FALSE
)
set.seed(123)
ggbetweenstats(
data = food,
x = stage,
y = energy_from_protein_percent_kcal,
type = "Robust",
plot.type = "box",
pairwise.comparisons = TRUE,
pairwise.display = "all",
xlab = "Cancer stage",
ylab = "Protein intake (% kcal)",
ggsignif.args = list(textsize = 5, tip_length = 0.03),
centrality.label.args = list(size = 4, nudge_x = 0.4, segment.linetype = 4),
ggstatsplot.layer = FALSE
)
set.seed(123)
ggbetweenstats(
data = food,
x = stage,
y = energy_from_carbohydrates_percent_kcal,
type = "Robust",
plot.type = "box",
pairwise.comparisons = TRUE,
pairwise.display = "all",
xlab = "Cancer stage",
ylab = "Carbohydrates intake (% kcal)",
ggsignif.args = list(textsize = 5, tip_length = 0.03),
centrality.label.args = list(size = 4, nudge_x = 0.4, segment.linetype = 4),
ggstatsplot.layer = FALSE
)
Bayesian analysis
Robust ANOVA (Muscle Loss)
In this case, we do not use trimmed mean for robust, but a t-distribution (Student) replacing the normal (Gaussian) distribution. We also use robust “priors” (median and MAD) in the model. https://bookdown.org/content/3686/metric-predicted-variable-on-one-or-two-groups.html
Energy from Protein (% kcal) and Muscle Loss
suppressPackageStartupMessages(library(brms))
# Variables priors
median_y = median(food$energy_from_carbohydrates_percent_kcal)
mad_y = mad(food$energy_from_carbohydrates_percent_kcal)
stanvars <-
stanvar(median_y, name = "median_y") +
stanvar(mad_y, name = "mad_y") +
stanvar(1/29, name = "one_over_twentynine")
# Model fit
fit.carb <-
brm(data = food,
family = student,
energy_from_carbohydrates_percent_kcal ~ 1 + muscle_loss,
prior = c(prior(normal(median_y, mad_y * 100), class = Intercept),
prior(normal(0, mad_y), class = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvars,
seed = 123
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Start sampling# Print results
summary(fit.carb)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: energy_from_carbohydrates_percent_kcal ~ 1 + muscle_loss
## Data: food (Number of observations: 32)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 45.24 2.27 40.75 49.68 1.00 3124 2811
## muscle_lossno -0.46 2.78 -5.99 5.00 1.00 3147 2917
## muscle_lossyes 5.09 2.74 -0.36 10.45 1.00 3235 3029
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.46 0.82 4.03 7.25 1.00 3261 2469
## nu 32.82 29.60 3.66 112.66 1.00 3136 2350
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot parameters distributions
pairs(fit.carb, off_diag_args = list(size = 1/3, alpha = 1/3))
suppressPackageStartupMessages(library(bayesplot))
posterior <- as.array(fit.carb)
dimnames(posterior)## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2" "chain:3" "chain:4"
##
## $parameters
## [1] "b_Intercept" "b_muscle_lossno" "b_muscle_lossyes" "sigma"
## [5] "nu" "lp__"# Plot posteriors (intervals)
mcmc_intervals(posterior, pars = c("b_muscle_lossno", "b_muscle_lossyes", "sigma"),
prob = 0.95)## Warning: `prob_outer` (0.9) is less than `prob` (0.95)
## ... Swapping the values of `prob_outer` and `prob`
# Plot posteriors (areas)
mcmc_areas(
posterior,
pars = c("b_muscle_lossno", "b_muscle_lossyes", "sigma"),
prob = 0.95, # 95% intervals
prob_outer = 0.99, # 99%
point_est = "mean"
)
Energy from Protein (% kcal) and Muscle Loss
# Variables priors
median_y = median(food$energy_from_protein_percent_kcal)
mad_y = mad(food$energy_from_protein_percent_kcal)
stanvars <-
stanvar(median_y, name = "median_y") +
stanvar(mad_y, name = "mad_y") +
stanvar(1/29, name = "one_over_twentynine")
# Model fit
fit.prot <-
brm(data = food,
family = student,
energy_from_protein_percent_kcal ~ 1 + muscle_loss,
prior = c(prior(normal(median_y, mad_y * 100), class = Intercept),
prior(normal(0, mad_y), class = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvars,
seed = 123
)## Compiling Stan program...## Start sampling# Print results
summary(fit.prot)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: energy_from_protein_percent_kcal ~ 1 + muscle_loss
## Data: food (Number of observations: 32)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 16.40 1.29 13.86 18.90 1.00 3550 2553
## muscle_lossno 0.05 1.61 -3.12 3.14 1.00 3337 2720
## muscle_lossyes -0.00 1.62 -3.10 3.31 1.00 3451 2805
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 3.22 0.44 2.47 4.20 1.00 4048 2731
## nu 39.95 30.26 5.70 121.77 1.00 4554 2861
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot parameters distributions
pairs(fit.prot, off_diag_args = list(size = 1/3, alpha = 1/3))
posterior <- as.array(fit.prot)
dimnames(posterior)## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2" "chain:3" "chain:4"
##
## $parameters
## [1] "b_Intercept" "b_muscle_lossno" "b_muscle_lossyes" "sigma"
## [5] "nu" "lp__"# Plot posteriors (intervals)
mcmc_intervals(posterior, pars = c("b_muscle_lossno", "b_muscle_lossyes", "sigma"),
prob = 0.95)## Warning: `prob_outer` (0.9) is less than `prob` (0.95)
## ... Swapping the values of `prob_outer` and `prob`
# Plot posteriors (areas)
mcmc_areas(
posterior,
pars = c("b_muscle_lossno", "b_muscle_lossyes", "sigma"),
prob = 0.95, # 95% intervals
prob_outer = 0.99, # 99%
point_est = "mean"
)
Robust ANOVA (Cancer Stage)
Energy from Protein (% kcal) and Cancer Stage
# Variables priors
median_y = median(food$energy_from_carbohydrates_percent_kcal)
mad_y = mad(food$energy_from_carbohydrates_percent_kcal)
stanvars <-
stanvar(median_y, name = "median_y") +
stanvar(mad_y, name = "mad_y") +
stanvar(1/29, name = "one_over_twentynine")
# Model fit
fit.carbStage <-
brm(data = food,
family = student,
energy_from_carbohydrates_percent_kcal ~ 1 + stage,
prior = c(prior(normal(median_y, mad_y * 100), class = Intercept),
prior(normal(0, mad_y), class = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvars,
seed = 123
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling# Print results
summary(fit.carbStage)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: energy_from_carbohydrates_percent_kcal ~ 1 + stage
## Data: food (Number of observations: 28)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 47.45 1.51 44.48 50.42 1.00 3709 2391
## stage2 -0.39 2.83 -5.89 5.21 1.00 3805 2478
## stage3 -5.62 3.17 -11.56 0.60 1.00 4045 2964
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.52 0.96 3.70 7.56 1.00 2501 1937
## nu 30.07 28.23 2.66 108.66 1.00 2764 1961
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot parameters distributions
pairs(fit.carbStage, off_diag_args = list(size = 1/3, alpha = 1/3))
posterior <- as.array(fit.carbStage)
dimnames(posterior)## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2" "chain:3" "chain:4"
##
## $parameters
## [1] "b_Intercept" "b_stage2" "b_stage3" "sigma" "nu"
## [6] "lp__"# Plot posteriors (intervals)
mcmc_intervals(posterior, pars = c("b_stage2", "b_stage3", "sigma"),
prob = 0.95)## Warning: `prob_outer` (0.9) is less than `prob` (0.95)
## ... Swapping the values of `prob_outer` and `prob`
# Plot posteriors (areas)
mcmc_areas(
posterior,
pars = c("b_stage2", "b_stage3", "sigma"),
prob = 0.95, # 95% intervals
prob_outer = 0.99, # 99%
point_est = "mean"
)
Energy from Protein (% kcal) and Cancer Stage
# Variables priors
median_y = median(food$energy_from_protein_percent_kcal)
mad_y = mad(food$energy_from_protein_percent_kcal)
stanvars <-
stanvar(median_y, name = "median_y") +
stanvar(mad_y, name = "mad_y") +
stanvar(1/29, name = "one_over_twentynine")
# Model fit
fit.protCancer <-
brm(data = food,
family = student,
energy_from_protein_percent_kcal ~ 1 + stage,
prior = c(prior(normal(median_y, mad_y * 100), class = Intercept),
prior(normal(0, mad_y), class = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvars,
seed = 123
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Start sampling# Print results
summary(fit.protCancer)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: energy_from_protein_percent_kcal ~ 1 + stage
## Data: food (Number of observations: 28)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 16.06 0.79 14.45 17.63 1.00 4518 2859
## stage2 1.39 1.55 -1.67 4.47 1.00 4269 3030
## stage3 2.99 1.77 -0.58 6.44 1.00 3974 3222
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 3.08 0.48 2.27 4.12 1.00 3597 2622
## nu 37.34 30.40 4.72 116.27 1.00 4103 2658
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Plot parameters distributions
pairs(fit.protCancer, off_diag_args = list(size = 1/3, alpha = 1/3))
posterior <- as.array(fit.protCancer)
dimnames(posterior)## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2" "chain:3" "chain:4"
##
## $parameters
## [1] "b_Intercept" "b_stage2" "b_stage3" "sigma" "nu"
## [6] "lp__"# Plot posteriors (intervals)
mcmc_intervals(posterior, pars = c("b_stage2", "b_stage3", "sigma"),
prob = 0.95)## Warning: `prob_outer` (0.9) is less than `prob` (0.95)
## ... Swapping the values of `prob_outer` and `prob`
# Plot posteriors (areas)
mcmc_areas(
posterior,
pars = c("b_stage2", "b_stage3", "sigma"),
prob = 0.95, # 95% intervals
prob_outer = 0.99, # 99%
point_est = "mean"
)
Age effect
suppressPackageStartupMessages(library(rstanarm))
# interaction effects with AGE
fit <- stan_glmer(energy_from_protein_percent_kcal ~ age + (1|muscle_loss), data = food)##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 7.5e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.229081 seconds (Warm-up)
## Chain 1: 0.275587 seconds (Sampling)
## Chain 1: 0.504668 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3.1e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.246781 seconds (Warm-up)
## Chain 2: 0.16653 seconds (Sampling)
## Chain 2: 0.413311 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 2.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.252729 seconds (Warm-up)
## Chain 3: 0.153166 seconds (Sampling)
## Chain 3: 0.405895 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.221091 seconds (Warm-up)
## Chain 4: 0.152203 seconds (Sampling)
## Chain 4: 0.373294 seconds (Total)
## Chain 4:## Warning: There were 7 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: Examine the pairs() plot to diagnose sampling problemsppc_intervals(
y = food$energy_from_protein_percent_kcal,
yrep = posterior_predict(fit),
x = food$age,
prob = 0.5
) +
labs(
x = "Age (years)",
y = "Energy from Protein (% kcal) "
) +
panel_bg(fill = "gray95", color = NA) +
grid_lines(color = "white")
fit2 <- stan_glmer(energy_from_carbohydrates_percent_kcal ~ age + (1|muscle_loss), data = food)##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4.2e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.291721 seconds (Warm-up)
## Chain 1: 0.232787 seconds (Sampling)
## Chain 1: 0.524508 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 3.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.334885 seconds (Warm-up)
## Chain 2: 0.226931 seconds (Sampling)
## Chain 2: 0.561816 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.306782 seconds (Warm-up)
## Chain 3: 0.258303 seconds (Sampling)
## Chain 3: 0.565085 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.330459 seconds (Warm-up)
## Chain 4: 0.241563 seconds (Sampling)
## Chain 4: 0.572022 seconds (Total)
## Chain 4:ppc_intervals(
y = food$energy_from_carbohydrates_percent_kcal,
yrep = posterior_predict(fit2),
x = food$age,
prob = 0.5
) +
labs(
x = "Age (years)",
y = "Energy from Carbohydrates (% kcal) "
) +
panel_bg(fill = "gray95", color = NA) +
grid_lines(color = "white")