IGF2 Attenuation of the Weight Loss–Regeneration Tradeoff
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)
library(broom.mixed)
library(patchwork)
theme_set(theme_classic())Load Data
dat <- read.csv("DataSheet - Body Measurements.csv")
weeks <- 1:8
long <- map_dfr(weeks, function(w){
tibble(
AnimalID = dat$AnimalID,
Treatment = dat$Treatment,
Week = factor(w),
PercentWL = dat[[paste0("PercWL_Week",w)]],
PercentRegeneration = dat[[paste0("RegenTail_Week",w)]]
)
}) %>%
filter(Treatment %in% c("Control","IGF2")) %>%
drop_na(PercentWL, PercentRegeneration)
long$Treatment <- factor(long$Treatment,
levels=c("Control","IGF2"))Mixed Model
mod <- lmer(
PercentRegeneration ~ PercentWL*Treatment*Week +
(1|AnimalID),
data=long
)
anova(mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## PercentWL 52.09 52.087 1 97.276 10.7810 0.001426 **
## Treatment 10.62 10.625 1 68.615 2.1992 0.142666
## Week 1832.47 305.412 6 97.335 63.2148 < 2.2e-16 ***
## PercentWL:Treatment 15.98 15.975 1 97.276 3.3066 0.072080 .
## PercentWL:Week 71.75 11.958 6 95.037 2.4751 0.028720 *
## Treatment:Week 81.50 13.584 6 97.335 2.8117 0.014464 *
## PercentWL:Treatment:Week 83.13 13.855 6 95.037 2.8677 0.013021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PercentRegeneration ~ PercentWL * Treatment * Week + (1 | AnimalID)
## Data: long
##
## REML criterion at convergence: 749.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0635 -0.4183 -0.0320 0.3533 3.1645
##
## Random effects:
## Groups Name Variance Std.Dev.
## AnimalID (Intercept) 27.722 5.265
## Residual 4.831 2.198
## Number of obs: 149, groups: AnimalID, 29
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.22241 3.00714 78.34909 1.404 0.164232
## PercentWL -0.02704 0.09504 102.29445 -0.284 0.776629
## TreatmentIGF2 2.80330 3.80894 75.07437 0.736 0.464037
## Week3 14.33756 1.96114 99.45166 7.311 6.82e-11
## Week4 22.44322 2.60271 95.33007 8.623 1.41e-13
## Week5 28.62670 2.77815 95.67149 10.304 < 2e-16
## Week6 35.34768 3.61411 96.82273 9.780 4.06e-16
## Week7 39.31685 4.28082 100.77821 9.184 5.71e-15
## Week8 36.55037 3.20867 95.03857 11.391 < 2e-16
## PercentWL:TreatmentIGF2 -0.07325 0.11970 107.85588 -0.612 0.541835
## PercentWL:Week3 -0.16749 0.06583 96.22178 -2.544 0.012539
## PercentWL:Week4 -0.21723 0.08298 96.36906 -2.618 0.010273
## PercentWL:Week5 -0.27472 0.09341 97.73982 -2.941 0.004083
## PercentWL:Week6 -0.38664 0.11114 94.96971 -3.479 0.000762
## PercentWL:Week7 -0.41289 0.11403 95.06457 -3.621 0.000473
## PercentWL:Week8 -0.38966 0.10567 95.31193 -3.688 0.000376
## TreatmentIGF2:Week3 -5.59486 2.51283 97.18683 -2.227 0.028289
## TreatmentIGF2:Week4 -6.19710 3.36454 96.21950 -1.842 0.068572
## TreatmentIGF2:Week5 -8.20697 3.56469 97.31334 -2.302 0.023449
## TreatmentIGF2:Week6 -13.68382 4.40997 99.50152 -3.103 0.002494
## TreatmentIGF2:Week7 -14.99155 5.24006 104.03072 -2.861 0.005106
## TreatmentIGF2:Week8 -11.88782 3.93527 96.67851 -3.021 0.003225
## PercentWL:TreatmentIGF2:Week3 0.21604 0.08914 95.26936 2.424 0.017250
## PercentWL:TreatmentIGF2:Week4 0.22577 0.10922 96.26236 2.067 0.041417
## PercentWL:TreatmentIGF2:Week5 0.28490 0.11923 97.85596 2.390 0.018786
## PercentWL:TreatmentIGF2:Week6 0.42156 0.13598 96.70669 3.100 0.002534
## PercentWL:TreatmentIGF2:Week7 0.41547 0.13759 96.20268 3.020 0.003242
## PercentWL:TreatmentIGF2:Week8 0.40170 0.12793 95.67857 3.140 0.002248
##
## (Intercept)
## PercentWL
## TreatmentIGF2
## Week3 ***
## Week4 ***
## Week5 ***
## Week6 ***
## Week7 ***
## Week8 ***
## PercentWL:TreatmentIGF2
## PercentWL:Week3 *
## PercentWL:Week4 *
## PercentWL:Week5 **
## PercentWL:Week6 ***
## PercentWL:Week7 ***
## PercentWL:Week8 ***
## TreatmentIGF2:Week3 *
## TreatmentIGF2:Week4 .
## TreatmentIGF2:Week5 *
## TreatmentIGF2:Week6 **
## TreatmentIGF2:Week7 **
## TreatmentIGF2:Week8 **
## PercentWL:TreatmentIGF2:Week3 *
## PercentWL:TreatmentIGF2:Week4 *
## PercentWL:TreatmentIGF2:Week5 *
## PercentWL:TreatmentIGF2:Week6 **
## PercentWL:TreatmentIGF2:Week7 **
## PercentWL:TreatmentIGF2:Week8 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itEstimated Slopes by Week
slopes <- emtrends(
mod,
pairwise ~ Treatment | Week,
var="PercentWL",
adjust="holm"
)
summary(slopes$emtrends)## Week = 2:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.0270 0.0969 103 -0.219 0.1652
## IGF2 -0.1003 0.0741 116 -0.247 0.0465
##
## Week = 3:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.1945 0.0835 103 -0.360 -0.0290
## IGF2 -0.0517 0.0815 116 -0.213 0.1097
##
## Week = 4:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.2443 0.0989 117 -0.440 -0.0484
## IGF2 -0.0918 0.0940 119 -0.278 0.0944
##
## Week = 5:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.3018 0.1030 121 -0.507 -0.0970
## IGF2 -0.0901 0.0995 119 -0.287 0.1069
##
## Week = 6:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.4137 0.1310 121 -0.672 -0.1553
## IGF2 -0.0654 0.1060 119 -0.275 0.1440
##
## Week = 7:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.4399 0.1350 121 -0.707 -0.1727
## IGF2 -0.0977 0.1030 119 -0.301 0.1057
##
## Week = 8:
## Treatment PercentWL.trend SE df lower.CL upper.CL
## Control -0.4167 0.1210 121 -0.657 -0.1766
## IGF2 -0.0882 0.0953 119 -0.277 0.1005
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Post Hoc Comparisons of Slopes
pairs(slopes$contrasts)## Week = 2:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 3:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 4:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 5:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 6:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 7:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Week = 8:
## contrast estimate SE df z.ratio p.value
## (nothing) nonEst NA NA NA NA
##
## Degrees-of-freedom method: kenward-rogerFDR-adjusted Comparisons
summary(
emtrends(
mod,
pairwise ~ Treatment | Week,
var="PercentWL",
adjust="fdr"
)$contrasts
)## Week = 2:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 0.0733 0.122 109 0.600 0.5495
##
## Week = 3:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.1428 0.117 110 -1.224 0.2236
##
## Week = 4:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.1525 0.136 118 -1.118 0.2660
##
## Week = 5:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.2116 0.144 120 -1.475 0.1429
##
## Week = 6:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.3483 0.168 121 -2.074 0.0402
##
## Week = 7:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.3422 0.170 121 -2.017 0.0459
##
## Week = 8:
## contrast estimate SE df t.ratio p.value
## Control - IGF2 -0.3284 0.154 121 -2.129 0.0353
##
## Degrees-of-freedom method: kenward-rogerSave Tables
write.csv(
as.data.frame(summary(slopes$emtrends)),
"Estimated_Slopes_by_Week.csv",
row.names=FALSE
)
write.csv(
as.data.frame(summary(slopes$contrasts)),
"Slope_Comparisons_Holm.csv",
row.names=FALSE
)Regression Plots
ggplot(long,
aes(PercentWL,
PercentRegeneration,
color=Treatment))+
geom_point(size=2)+
geom_smooth(method="lm",se=TRUE)+
facet_wrap(~Week)+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
Estimated Slopes Figure
slope.df <- as.data.frame(summary(slopes$emtrends))
ggplot(slope.df,
aes(Week, PercentWL.trend,
color=Treatment,
group=Treatment))+
geom_point(size=3)+
geom_line()+
geom_errorbar(aes(ymin=lower.CL,
ymax=upper.CL),
width=.15)+
ylab("Estimated slope of regeneration vs weight loss")
library(emmeans)
library(dplyr)
library(ggplot2)
# Estimated weight-loss slopes for each treatment at each week
slope_results <- emtrends(
mod,
specs = ~ Treatment | Week,
var = "PercentWL"
)
# Compare slopes as IGF2 - Control within each week
slope_difference <- contrast(
slope_results,
method = list(
"IGF2 - Control" = c(-1, 1)
),
by = "Week",
adjust = "holm"
)
# Convert results to a data frame with confidence intervals
difference_df <- as.data.frame(
summary(
slope_difference,
infer = c(TRUE, TRUE)
)
) %>%
mutate(
Week = as.numeric(as.character(Week)),
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.10 ~ "\u2020",
TRUE ~ ""
),
p_label = case_when(
p.value < 0.001 ~ "p < .001",
TRUE ~ paste0("p = ", sprintf("%.3f", p.value))
)
)
difference_df## contrast Week estimate SE df lower.CL upper.CL
## 1 IGF2 - Control 2 -0.07325351 0.1220145 108.5699 -0.315093122 0.1685861
## 2 IGF2 - Control 3 0.14279127 0.1166659 109.8779 -0.088416020 0.3739986
## 3 IGF2 - Control 4 0.15251344 0.1364541 118.1933 -0.117698323 0.4227252
## 4 IGF2 - Control 5 0.21164809 0.1435180 120.3995 -0.072498015 0.4957942
## 5 IGF2 - Control 6 0.34830594 0.1679447 120.6860 0.015806449 0.6808054
## 6 IGF2 - Control 7 0.34221926 0.1696304 120.7564 0.006384319 0.6780542
## 7 IGF2 - Control 8 0.32844599 0.1542723 120.8273 0.023018783 0.6338732
## t.ratio p.value significance p_label
## 1 -0.6003671 0.54951277 p = 0.550
## 2 1.2239332 0.22359542 p = 0.224
## 3 1.1176901 0.26596609 p = 0.266
## 4 1.4747143 0.14289949 p = 0.143
## 5 2.0739328 0.04021235 * p = 0.040
## 6 2.0174408 0.04586663 * p = 0.046
## 7 2.1290011 0.03528402 * p = 0.035# Position p-value labels above the confidence interval
difference_df <- difference_df %>%
mutate(
label_y = upper.CL + 0.06
)
p_difference <- ggplot(
difference_df,
aes(
x = Week,
y = estimate
)
) +
annotate(
"rect",
xmin = -Inf,
xmax = Inf,
ymin = 0,
ymax = Inf,
fill = "#DDE2E8",
alpha = 0.85
) +
annotate(
"rect",
xmin = -Inf,
xmax = Inf,
ymin = -Inf,
ymax = 0,
fill = "#FFF0E6",
alpha = 0.85
) +
geom_hline(
yintercept = 0,
linetype = "dashed",
linewidth = 0.8
) +
geom_line(
group = 1,
linewidth = 1.1,
color = "#102A4C"
) +
geom_errorbar(
aes(
ymin = lower.CL,
ymax = upper.CL
),
width = 0.12,
linewidth = 1.1,
color = "#102A4C"
) +
geom_point(
size = 4.5,
color = "#102A4C"
) +
geom_text(
aes(
y = label_y,
label = significance
),
size = 7,
fontface = "bold",
color = "#102A4C",
na.rm = TRUE
) +
scale_x_continuous(
breaks = sort(unique(difference_df$Week))
) +
labs(
title = "IGF2 Attenuates the Energetic Constraint on Regeneration",
subtitle = expression(
"Difference in weight-loss slopes: " *
Delta * beta == beta[IGF2] - beta[Control]
),
x = "Week",
y = expression(
Delta * " slope relative to control"
),
caption = paste0(
"Positive values indicate that the relationship between weight loss ",
"and regeneration is less negative under IGF2. ",
"Error bars are 95% confidence intervals. ",
"Holm-adjusted comparisons: \u2020p < .10, *p < .05, **p < .01, ***p < .001."
)
) +
theme_classic(base_size = 18) +
theme(
plot.title = element_text(
face = "bold",
size = 24
),
plot.subtitle = element_text(
face = "bold",
size = 17,
margin = margin(b = 12)
),
axis.title = element_text(
face = "bold",
size = 18
),
axis.text = element_text(
size = 15,
color = "black"
),
plot.caption = element_text(
size = 11,
hjust = 0
)
)
p_difference
p_difference <- p_difference +
annotate(
"label",
x = max(difference_df$Week) - 1,
y = max(difference_df$upper.CL, na.rm = TRUE) + 0.18,
label = "Greater Attenuation",
color = "#102A4C",
fill = "white",
fontface = "bold",
size = 6
) +
annotate(
"label",
x = max(difference_df$Week) - 1,
y = min(difference_df$lower.CL, na.rm = TRUE) - 0.12,
label = "No Attenuation",
color = "#E87517",
fill = "white",
fontface = "bold",
size = 6
)
p_difference
ggsave(
"IGF2_Difference_in_Slopes.png",
p_difference,
width = 10,
height = 7,
units = "in",
dpi = 600,
bg = "white"
)
ggsave(
"IGF2_Difference_in_Slopes.pdf",
p_difference,
width = 10,
height = 7,
units = "in"
)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Positive values indicate that the relationship between
## weight loss and regeneration is less negative under IGF2. Error bars are 95%
## confidence intervals. Holm-adjusted comparisons: †p < .10, *p < .05, **p < .01,
## ***p < .001.' in 'mbcsToSbcs': for † (U+2020)# library(patchwork)
#
# combined_figure <- p_slopes / p_difference +
# plot_annotation(
# tag_levels = "A",
# theme = theme(
# plot.tag = element_text(
# face = "bold",
# size = 20
# )
# )
# )
#
# combined_figure
#
# ggsave(
# "IGF2_Slopes_and_Slope_Differences.png",
# combined_figure,
# width = 11,
# height = 13,
# units = "in",
# dpi = 600,
# bg = "white"
# )library(emmeans)
library(dplyr)
library(ggplot2)
# Estimated weight-loss slopes for each treatment at each week
slope_results <- emtrends(
mod,
specs = ~ Treatment | Week,
var = "PercentWL"
)
# Compare slopes as IGF2 - Control within each week
slope_difference <- contrast(
slope_results,
method = list(
"IGF2 - Control" = c(-1, 1)
),
by = "Week",
adjust = "holm"
)
# Convert results to a data frame with confidence intervals
difference_df <- as.data.frame(
summary(
slope_difference,
infer = c(TRUE, TRUE)
)
) %>%
mutate(
Week = as.numeric(as.character(Week)),
significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.10 ~ "\u2020",
TRUE ~ ""
),
p_label = case_when(
p.value < 0.001 ~ "p < .001",
TRUE ~ paste0("p = ", sprintf("%.3f", p.value))
)
)
difference_df## contrast Week estimate SE df lower.CL upper.CL
## 1 IGF2 - Control 2 -0.07325351 0.1220145 108.5699 -0.315093122 0.1685861
## 2 IGF2 - Control 3 0.14279127 0.1166659 109.8779 -0.088416020 0.3739986
## 3 IGF2 - Control 4 0.15251344 0.1364541 118.1933 -0.117698323 0.4227252
## 4 IGF2 - Control 5 0.21164809 0.1435180 120.3995 -0.072498015 0.4957942
## 5 IGF2 - Control 6 0.34830594 0.1679447 120.6860 0.015806449 0.6808054
## 6 IGF2 - Control 7 0.34221926 0.1696304 120.7564 0.006384319 0.6780542
## 7 IGF2 - Control 8 0.32844599 0.1542723 120.8273 0.023018783 0.6338732
## t.ratio p.value significance p_label
## 1 -0.6003671 0.54951277 p = 0.550
## 2 1.2239332 0.22359542 p = 0.224
## 3 1.1176901 0.26596609 p = 0.266
## 4 1.4747143 0.14289949 p = 0.143
## 5 2.0739328 0.04021235 * p = 0.040
## 6 2.0174408 0.04586663 * p = 0.046
## 7 2.1290011 0.03528402 * p = 0.035# Position p-value labels above the confidence interval
difference_df <- difference_df %>%
mutate(
label_y = upper.CL + 0.06
)
p_difference <- ggplot(
difference_df,
aes(
x = Week,
y = estimate
)
) +
annotate(
"rect",
xmin = -Inf,
xmax = Inf,
ymin = 0,
ymax = Inf,
fill = "#DDE2E8",
alpha = 0.85
) +
annotate(
"rect",
xmin = -Inf,
xmax = Inf,
ymin = -Inf,
ymax = 0,
fill = "#FFF0E6",
alpha = 0.85
) +
geom_hline(
yintercept = 0,
linetype = "dashed",
linewidth = 0.8
) +
geom_line(
group = 1,
linewidth = 1.1,
color = "#102A4C"
) +
geom_errorbar(
aes(
ymin = lower.CL,
ymax = upper.CL
),
width = 0.12,
linewidth = 1.1,
color = "#102A4C"
) +
geom_point(
size = 4.5,
color = "#102A4C"
) +
geom_text(
aes(
y = label_y,
label = significance
),
size = 7,
fontface = "bold",
color = "#102A4C",
na.rm = TRUE
) +
scale_x_continuous(
breaks = sort(unique(difference_df$Week))
) +
labs(
title = "IGF2 Attenuates the Energetic Constraint on Regeneration",
subtitle = expression(
"Difference in weight-loss slopes: " *
Delta * beta == beta[IGF2] - beta[Control]
),
x = "Week",
y = expression(
Delta * " slope relative to control"
),
caption = paste0(
"Positive values indicate that the relationship between weight loss ",
"and regeneration is less negative under IGF2. ",
"Error bars are 95% confidence intervals. ",
"Holm-adjusted comparisons: \u2020p < .10, *p < .05, **p < .01, ***p < .001."
)
) +
theme_classic(base_size = 18) +
theme(
plot.title = element_text(
face = "bold",
size = 24
),
plot.subtitle = element_text(
face = "bold",
size = 17,
margin = margin(b = 12)
),
axis.title = element_text(
face = "bold",
size = 18
),
axis.text = element_text(
size = 15,
color = "black"
),
plot.caption = element_text(
size = 11,
hjust = 0
)
)
p_difference
p_difference <- p_difference +
annotate(
"label",
x = max(difference_df$Week) - 1,
y = max(difference_df$upper.CL, na.rm = TRUE) + 0.18,
label = "Greater Attenuation",
color = "#102A4C",
fill = "white",
fontface = "bold",
size = 6
) +
annotate(
"label",
x = max(difference_df$Week) - 1,
y = min(difference_df$lower.CL, na.rm = TRUE) - 0.12,
label = "No Attenuation",
color = "#E87517",
fill = "white",
fontface = "bold",
size = 6
)
p_difference
ggsave(
"IGF2_Difference_in_Slopes.png",
p_difference,
width = 10,
height = 7,
units = "in",
dpi = 600,
bg = "white"
)
ggsave(
"IGF2_Difference_in_Slopes.pdf",
p_difference,
width = 10,
height = 7,
units = "in"
)## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## conversion failure on 'Positive values indicate that the relationship between
## weight loss and regeneration is less negative under IGF2. Error bars are 95%
## confidence intervals. Holm-adjusted comparisons: †p < .10, *p < .05, **p < .01,
## ***p < .001.' in 'mbcsToSbcs': for † (U+2020)# library(patchwork)
#
# combined_figure <- p_slopes / p_difference +
# plot_annotation(
# tag_levels = "A",
# theme = theme(
# plot.tag = element_text(
# face = "bold",
# size = 20
# )
# )
# )
#
# combined_figureDiagnostics
plot(mod)
qqnorm(resid(mod))
qqline(resid(mod))
Optional Random Slope Model
# mod.rs <- lmer(
# PercentRegeneration ~ PercentWL*Treatment*Week +
# (1 + Week|AnimalID),
# data=long,
# REML=FALSE
# )
# anova(mod,mod.rs)