Regeneration.Diet.IGFs

Load Packages

#load all packages necessary to run statistical and graphic models
library(ggplot2)
library(nlme)
library(multcomp)
library(emmeans)
library(MuMIn)
library(bestNormalize)
library(lmtest)
library(car)
library(Rmisc)
library(dplyr)
library(viridis)
library(hrbrthemes)
library(plotly)

#Set position dodge location for categorized boxplot graphs and data point locations
pd <- position_dodge(0.02)
dodge <- position_dodge(width = 0.6)

Part 1: Pre-experimental body size comparisons

Make sure there are no differences in body size between treatments before starting the experiment

#Read in data set that includes the SVL, mass, and treatment assignment of each individual immediately before beginning the experiment
diet=read.csv("WL.analysis.csv")

#subset data to only include the original body measurements
pre.diet=subset(diet, Week == "-4")

#Run linear model to determine if there was a statistical difference in SVL at the time of beginning the experiment
SVL=lm(SVL~Treatment, data=pre.diet)
#print linear model results
anova(SVL)

## Analysis of Variance Table
## 
## Response: SVL
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Treatment  3   6.984  2.3279  0.6542 0.5829
## Residuals 72 256.214  3.5585

#plot the relationship of SVL across treatments
svl.pl=ggplot(pre.diet, aes(x=as.factor(Treatment), y=SVL, fill=Treatment)) +
    #violin plot displays the density of data at each point displaying the mix and max of each treatment as well
    geom_violin(trim=F, position= dodge, scale="width") +
scale_fill_manual(values = c("slategrey", "lemonchiffon2", "thistle4", "mistyrose3")) +  #insert boxplot overlay displaying the quantiles and means of each group
  geom_boxplot(width=0.1, position= dodge, outlier.shape = NA) +
  geom_jitter(width=0.1) +
  xlab("Treatment")

svl.pl

ggsave(svl.pl, file="svl.png", width=5, height=4, dpi=600)

Result: There was no statistical difference in body size in SVL at the beginning of the experimental period between treatments.

#Run linear model to determine if there was a statistical difference in SVL at the time of beginning the experiment
Mass=lm(Mass~Treatment, data=pre.diet)
#print linear model results
anova(Mass)

## Analysis of Variance Table
## 
## Response: Mass
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Treatment  3  0.3761 0.12538  0.5547 0.6467
## Residuals 72 16.2743 0.22603

#Plot the relationship of SVL across treatments
mass.pl=ggplot(pre.diet, aes(x=as.factor(Treatment), y=Mass, fill=Treatment)) + 
     #violin plot displays the density of data at each point displaying the mix and max of each treatment as well
    geom_violin(trim=F, position= dodge, scale="width") +
scale_fill_manual(values = c("slategrey", "lemonchiffon2", "thistle4", "mistyrose3")) +  geom_boxplot(width=0.1, position= dodge, outlier.shape = NA) +
  geom_jitter(width=0.1) +
  xlab("Treatment")

ggsave(mass.pl, file="mass.png", width=5, height=4, dpi=600)

Result: There was no statistical difference in body size in mass at the beginning of the experimental period between treatments.

Part 2: Assess weight lost post-diet implementation

Weight Loss post-diet implementation but before beginning the experiment

#Look at data set including the body measurements (mass) over time before and after diet implementation. This data set is designed to verify that the diet did result in weight loss in the DR group, and that weight loss levels off before the experiment begins. We do not want to begin the experiment with some animals actively losing weight at a high rate, and others not.

#Subset data to be only the data collected immediately before the experiment begins (week -1). 
diet.an=subset(diet, Week == -1)
#Run linear model determining if the group on a diet lost more weight than the group on the ad lib food 
diet.lm=lm(Perc_WL~Diet, data=diet.an)
#print the results
anova(diet.lm)

## Analysis of Variance Table
## 
## Response: Perc_WL
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Diet       1  608.22  608.22  13.854 0.0004068 ***
## Residuals 67 2941.39   43.90                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Plot the relationship of percent weight loss between the two groups (diet restricted and ad lib)
diet.pl=ggplot(diet.an, aes(x=as.factor(Diet), y=Perc_Orig, fill=Diet)) + 
     #violin plot displays the density of data at each point displaying the mix and max of each treatment as well
    geom_violin(trim=F, position= dodge, scale="width") +
scale_fill_manual(values = c("slategrey", "lemonchiffon2")) +    #insert boxplot overlay displaying the quantiles and means of each group
  geom_boxplot(width=0.1, position= dodge, outlier.shape = NA) +
geom_jitter(width=0.1) +
  xlab("Diet")

ggsave(diet.pl, file="diet.png", width=5, height=4, dpi=600)

diet.pl

Result: The animals on a restricted diet lost significantly more of their body weight than those from the ad lib diet. This indicates that the diet was effective at calorically restricting the animals. The animals on the ad lib diet lost <5% of their total body weight on average, indicating that the ad lib diet was truly ad lib.

#Plot the relationship of percent weight loss between the two groups (diet restricted and ad lib) over time
diet2.pl=ggplot(diet, aes(x=as.factor(Week), y=Mass, fill=Diet)) + 
     #violin plot displays the density of data at each point displaying the mix and max of each treatment as well
    geom_violin(trim=F, position= dodge, scale="width", width=0.5) +
scale_fill_manual(values = c("slategrey", "lemonchiffon2")) +    #insert boxplot overlay displaying the quantiles and means of each group
  geom_boxplot(width=0.1, position= dodge, outlier.shape = NA) +
geom_jitter(position= dodge) +
  xlab("Week")

ggsave(diet2.pl, file="diet2.png", width=5, height=4, dpi=600)

#plot this relationship as a line graph
diet2.pl=ggplot(diet, aes(x=Week, y=Mass, colour=Diet)) + 
  geom_smooth() +
    scale_color_manual(values = c("slategrey", "lemonchiffon2")) 

diet3.pl=ggplot(diet, aes(x=Week, y=Perc_WL, colour=Diet)) + 
  geom_smooth() +
    scale_color_manual(values = c("slategrey", "lemonchiffon2")) 

diet4.pl=ggplot(diet, aes(x=Week, y=Perc_Orig, colour=Diet)) + 
  geom_smooth() +
    scale_color_manual(values = c("slategrey", "lemonchiffon2")) 

ggsave(diet4.pl, file="diet4.png", width=4.75, height=4, dpi=600)

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

diet2.pl

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

diet3.pl

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

diet4.pl

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

Result: There was weight fluctuation in the early week, but the weight loss had tapered off by the beginning of the experiment and steadied.

Part 3: Unedited Generalized Visualization of Data

Outlier analysis and analysis of reduction in regeneration data points

#Length Loss Analysis

ll=read.csv("LengthLossAnalysis.csv")
ll_sub=subset(ll, Week>2)
reg <- read.csv("R.analysis.currated.csv", stringsAsFactors = FALSE)


pl.ll=ggplot(reg, aes(x=Week, y=Regeneration, colour=Animal_ID)) +
  geom_line(aes(group = Animal_ID)) +
    facet_wrap(~Treatment, scales="free", ncol=2) +
  geom_text(aes(label=Animal_ID),hjust=0, vjust=0, size=2)+
  theme(legend.position = "none")
  
pl.ll

ggsave(pl.ll, file="length_loss.png", width=8, height=8, dpi=600)

pl.ll.sub=ggplot(ll_sub, aes(x=Week, y=Regeneration, colour=Animal_ID)) +
  geom_line(aes(group = Animal_ID)) +
    facet_wrap(~Treatment, scales="free", ncol=2) +
    geom_text(aes(label=Animal_ID),hjust=0, vjust=0, size=2)+
  theme(legend.position = "none")
pl.ll.sub

ggsave(pl.ll.sub, file="length_loss_sub.png", width=8, height=8, dpi=600)

Result: There are clear inconsistencies following the 4 week timepoint. Beginning at week 5, there were points in time where there were decreases in regeneration length (meaning tail length was lost). Therefore, we decided to zoom in to the first four weeks of time. Additionally, this decision is supported by the Beatty, Mote, Schwartz Regeneration paper indicating the majority of investment into tail regeneration occurs during the first 4 weeks of tail regeneration. There were no additional inconsistencies in the tail regeneration patterns within the first 4 weeks.

Overall visualization of data

#upload raw data set in long format with all 8 weeks of regeneration, reproduction, and weight data following the acclimation period.
reg=read.csv("R.analysis.currated.csv")

#The primary data set used in analysis will be weeks 1-4. We chose this because our previous study (Beatty, Mote and Schwartz 2020) shows that nearly all investment into regeneration occurs in the first four weeks of regeneration. Additionally, the study was performed late in the season, as reproduction tapered off, and there were issues with measurements, COVID, and feeding during the end of the experiment.
early_reg=subset(reg, Week < 5)

pl.reg.rate=ggplot(early_reg, aes(x=Week, y=Rate_Reg, colour=Animal_ID)) +
  geom_line(aes(group = Animal_ID)) +
    facet_wrap(~Treatment, scales="free", ncol=2) +
  geom_text(aes(label=Animal_ID),hjust=0, vjust=0, size=2)+
  theme(legend.position = "none")
  
pl.reg.rate

## Warning: Removed 121 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 121 rows containing missing values or values outside the scale range
## (`geom_text()`).

pl.perc.reg=ggplot(early_reg, aes(x=Week, y=Perc.Regeneration, colour=Animal_ID)) +
  geom_line(aes(group = Animal_ID)) +
    facet_wrap(~Treatment, scales="free", ncol=2) +
  geom_text(aes(label=Animal_ID),hjust=0, vjust=0, size=2)+
  theme(legend.position = "none")
  
pl.perc.reg

## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 54 rows containing missing values or values outside the scale range
## (`geom_text()`).

Result: Analyses were previously completed using the rate of tail regeneration rather than the regenerative length. This showed that there was a large spike in regenerative investment between weeks 2 and 3, which then decreased significantly again by week 4 (also consistant with our previous findings). This was the case for nearly all individials regardless of treatment assignment. The validation of regeneration using this metric was a sanity check. For analyses, percent regneration is used, as it accounts for differences in total amount needed to regenerate, and possible variations in body size of the animals (volume of tissue needed to regenerate varies between animal).

outlierTest(lm(Rate_Reg~ Treatment*Week, random=~1|Animal_ID , data=early_reg, na.action=na.omit))

## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'random' will be disregarded

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 169 2.990036          0.0031587      0.62857

outlierTest(lm(Perc.Regeneration~ Treatment*Week, random=~1|Animal_ID , data=early_reg, na.action=na.omit))

## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'random' will be disregarded

##     rstudent unadjusted p-value Bonferroni p
## 274 4.241648         3.0990e-05    0.0082432
## 275 4.217107         3.4316e-05    0.0091280
## 277 4.051736         6.7389e-05    0.0179250

Result: When using percent regeneration, BC333 and BC335 at week 4 are statistical outliers. These are the two highest individuals in the ad lib group. Additionally, in the IGF2 group individual BC337 is labeled as an outlier at week 4. Using the corrected p-values, the rate of regeneration does not have any statistical outliers. These samples did not have any questionable data, and followed reliable trends indicating that they are statistical outliers, but biologically accurate. The individual of concern (BC327) was not indicated as an outlier in either test. Therefore, no data points were excluded as data points. Again, regeneration within the first 4 weeks appears to be consistent across individuals and a reliable measure for analysis.

Preliminary reproduction analysis:

Plot of total reproduction by time to justify sperm storage and reference our other regeneration paper. This is another way to test the seasonal effects in our dataset as well. NOTE: we did not reliably collect eggs during the acclimation period, so we cannot produce a plot that includes that data.

ggplot(reg, aes(x = factor(Week), y = Eggs, group=Treatment, color=Treatment)) + 
  geom_line(stat = "summary", fun = "mean", position="dodge", size=2)

repro.pl=ggplot(reg, aes(x = factor(Week), y = Eggs, group=Treatment, color=Treatment)) + 
  geom_line(stat = "summary", fun = "sum", position="dodge", size=2) +
  scale_color_manual(values = c("slategrey", "lemonchiffon2", "thistle4", "mistyrose3")) +
      theme_ipsum() 

repro.pl

ggsave(repro.pl, file="reproduction.png", width=4, height=4, dpi=600)

There is alot of variation in regeneration over the time period. Whether you are looking at the average number of eggs, or the sum of egg production between treatments, there is a spike at week 2 in all treatments except for IGF2. Then they all decreased to zero for week 3, with the exception of IGF2. The patterns remain inconsistent through week 8. Should we include egg production in the model? I dont think we should because the total is so low, personally.

RQ1- Does restricted diet decrease regeneration rate?

#Analysis 2. Does a restricted diet decrease rate of tail regeneration. >rate~treatment*week + rand (ID) using the adlib and vehicle groups only, as both were given saline. >repro~treatment: using the sum of all individuals across all 8 weeks for reproduction. Do this ONLY if reproduction was important in the preliminary analysis (looking for seasonal effects).

diet.reg=subset(early_reg, Treatment == "DR" | Treatment == "AdLib")
anova(lme(Perc.Regeneration~ Treatment*Week, random=~1|Animal_ID , data=diet.reg, na.action=na.omit))

##                numDF denDF  F-value p-value
## (Intercept)        1    99 218.8295  <.0001
## Treatment          1    32   2.9581  0.0951
## Week               1    99 376.6708  <.0001
## Treatment:Week     1    99   5.3101  0.0233

anova(lme(Eggs~ Treatment, random=~1|Animal_ID , data=diet.reg, na.action=na.omit))

##             numDF denDF  F-value p-value
## (Intercept)     1   101 8.090978  0.0054
## Treatment       1    32 0.172737  0.6805

#get week by week comparisons
w1=subset(diet.reg, Week == "1")
w2=subset(diet.reg, Week == "2")
w3=subset(diet.reg, Week == "3")
w4=subset(diet.reg, Week == "4")

anova(lm(Perc.Regeneration~ Treatment, data=w2, na.action=na.omit))

## Analysis of Variance Table
## 
## Response: Perc.Regeneration
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Treatment  1 0.00366 0.003665  0.1075 0.7452
## Residuals 31 1.05638 0.034077

anova(lm(Perc.Regeneration~ Treatment, data=w3, na.action=na.omit))

## Analysis of Variance Table
## 
## Response: Perc.Regeneration
##           Df Sum Sq Mean Sq F value Pr(>F)
## Treatment  1  17.90  17.903  1.0569 0.3116
## Residuals 32 542.05  16.939

anova(lm(Perc.Regeneration~ Treatment, data=w4, na.action=na.omit))

## Analysis of Variance Table
## 
## Response: Perc.Regeneration
##           Df  Sum Sq Mean Sq F value  Pr(>F)  
## Treatment  1  147.96  147.96  3.9216 0.05632 .
## Residuals 32 1207.37   37.73                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

diet.reg.av=Rmisc::summarySE(data=diet.reg, measurevar="Perc.Regeneration", groupvars=c("Week", "Treatment"), na.rm=T, conf.interval=0.95)



diet.reg.pl=ggplot(diet.reg, aes(x=Week, y=Perc.Regeneration, color=Treatment)) + 
  geom_smooth(linewidth=2) +
      scale_color_manual(values = c("slategrey", "lemonchiffon2")) +
        theme(legend.position="top") +
    ylim(0,16) 

diet.reg.pl

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

ggsave(diet.reg.pl, file="diet.reg.pl.png", height=4, width=4, dpi = 300)

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

RQ2-Does IGF injection increase regeneration rate?

#Analysis 3. Does Increasing levels of IGF1 or IGF2 will increase the rate of tail regeneration >rate~treat*week + rand(ID) using the saline, IGF1, and IGF2 groups. Consider making density plots. Then complete emmeans pairwise analyses between treatments at each week.

inj.reg=subset(early_reg, Treatment != "AdLib")
anova(lme(Perc.Regeneration~ Treatment*Week, random=~1|Animal_ID , data=inj.reg, na.action=na.omit))

##                numDF denDF  F-value p-value
## (Intercept)        1   141 303.4932  <.0001
## Treatment          2    48   1.0038  0.3741
## Week               1   141 522.9761  <.0001
## Treatment:Week     2   141   0.5019  0.6065

inj.reg.average=Rmisc::summarySE(data=inj.reg, measurevar="Perc.Regeneration", groupvars=c("Week","Treatment"), na.rm=T, conf.interval=0.95)

inj.reg.pl=ggplot(inj.reg, aes(x=Week, y=Perc.Regeneration, color=Treatment)) +  
    geom_smooth(linewidth=2) +
     # geom_jitter(data=inj.reg, alpha=1, size=0.5) +
      scale_color_manual(values = c("lemonchiffon2", "thistle4", "mistyrose3")) +
        theme(legend.position="top") +
    ylim(0,16) 

inj.reg.pl

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

ggsave(inj.reg.pl, file="inj.reg.pl.png", height=4, width=4, dpi = 300)

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

Additional Plot Types

ggplot(diet.reg,aes(x=as.factor(Week), y=Perc.Regeneration, fill=Treatment)) + 
  geom_boxplot()  + 
      scale_fill_manual(values = c("slategrey", "lemonchiffon2")) +
          theme(legend.position="top") +
  facet_wrap(.~Week, nrow=1, scales="free") +
  ylim(0,35)

## Warning: Removed 25 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

ggplot(inj.reg,aes(x=as.factor(Week), y=Perc.Regeneration, fill=Treatment)) + 
  geom_boxplot()  + 
      scale_fill_manual(values = c("lemonchiffon2", "thistle4", "mistyrose3")) +
          theme(legend.position="top") +
  facet_wrap(.~Week, nrow=1, scales="free") +
  ylim(0,30)

## Warning: Removed 45 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

Print alternative view of results

ggplot(diet.reg.av, aes(x=Week, y=Perc.Regeneration, color=Treatment, fill=Treatment)) +
    geom_area(position = "identity", alpha=0.02, size=2) +
    geom_jitter(data=diet.reg, alpha=0.9, size=0.5) +
      scale_color_manual(values = c("slategrey", "lemonchiffon2")) +
        scale_fill_manual(values = c("slategrey", "lemonchiffon2")) +
      ggtitle("Effect of Diet Restriction") +
    theme_ipsum() +
        theme(legend.position="top") 

## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(inj.reg.average, aes(x=Week, y=Perc.Regeneration, color=Treatment, fill=Treatment)) +
       geom_area(position = "identity", alpha=0.02, size=2 ) +
    geom_jitter(data=inj.reg, alpha=0.9, size=0.5) +
      scale_color_manual(values = c("lemonchiffon2", "thistle4", "mistyrose3")) +
        scale_fill_manual(values = c("lemonchiffon2", "thistle4", "mistyrose3")) +
    ggtitle("Effect of Supplemental IGF Injection") +
    theme_ipsum() +
        theme(legend.position="top",
              plot.title = element_text(size=14))

## Warning: Removed 45 rows containing missing values or values outside the scale range
## (`geom_point()`).

Response to Reviews

# Is there a relationship between change in Body Size and Regeneration rate
library(dplyr)
library(ggplot2)
library(lme4)

## Loading required package: Matrix

## Registered S3 method overwritten by 'lme4':
##   method           from
##   na.action.merMod car

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:nlme':
## 
##     lmList

data <- read.csv("R.analysis.currated.csv", stringsAsFactors = FALSE)
colnames(data) <- trimws(colnames(data))

analysis_data <- data %>%
  group_by(Animal_ID) %>%
  summarize(
    mean_weight_loss = mean(Perc.WeightLoss, na.rm = TRUE),
    mean_regen_rate  = mean(Rate_Reg, na.rm = TRUE),
    Treatment        = first(Treatment)
  )

model_mixed <- lmer(Rate_Reg ~ Perc.WeightLoss + (1 | Animal_ID), data = data)
summary(model_mixed)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Rate_Reg ~ Perc.WeightLoss + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 1935.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1316 -0.7990 -0.2811  0.5622  5.9167 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 0.06057  0.2461  
##  Residual              7.38336  2.7172  
## Number of obs: 398, groups:  Animal_ID, 68
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      2.50071    0.28336   8.825
## Perc.WeightLoss  0.01787    0.01440   1.241
## 
## Correlation of Fixed Effects:
##             (Intr)
## Prc.WghtLss -0.870

Result: Does weight loss predict regeneration rate, while accounting for repeated measurements per animal? There is no strong evidence that percent weight loss affects tail regeneration rate. A 1% increase in weight loss → only ~0.018 increase in regeneration rate. Regeneration rate varies a lot within individuals over time. Weight loss isn’t explaining much of that variation

summary(lmer(Rate_Reg ~ Perc.WeightLoss * Treatment + (1 | Animal_ID), data = data))

## Linear mixed model fit by REML ['lmerMod']
## Formula: Rate_Reg ~ Perc.WeightLoss * Treatment + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 1937.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3983 -0.7633 -0.2416  0.5594  6.0244 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 0.08015  0.2831  
##  Residual              7.26570  2.6955  
## Number of obs: 398, groups:  Animal_ID, 68
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                       2.978100   0.374034   7.962
## Perc.WeightLoss                   0.034070   0.023723   1.436
## TreatmentDR                      -0.414679   0.975090  -0.425
## TreatmentDR.IGF1                 -0.842106   0.754979  -1.115
## TreatmentDR.IGF2                 -2.119445   0.896293  -2.365
## Perc.WeightLoss:TreatmentDR      -0.028237   0.045678  -0.618
## Perc.WeightLoss:TreatmentDR.IGF1 -0.006972   0.044263  -0.158
## Perc.WeightLoss:TreatmentDR.IGF2  0.052077   0.044466   1.171
## 
## Correlation of Fixed Effects:
##               (Intr) Prc.WL TrtmDR TDR.IGF1 TDR.IGF2 Pr.WL:TDR P.WL:TDR.IGF1
## Prc.WghtLss   -0.675                                                        
## TreatmentDR   -0.384  0.259                                                 
## TrtmDR.IGF1   -0.495  0.335  0.190                                          
## TrtmDR.IGF2   -0.417  0.282  0.160  0.207                                   
## Prc.WgL:TDR    0.351 -0.519 -0.881 -0.174   -0.146                          
## P.WL:TDR.IGF1  0.362 -0.536 -0.139 -0.841   -0.151    0.278                 
## P.WL:TDR.IGF2  0.360 -0.534 -0.138 -0.179   -0.876    0.277     0.286

Result: When you account for treatment, No strong evidence that weight loss predicts regeneration rate. No evidence for a tradeoff (energy vs regeneration) that depends on treatment.Treatments affect regeneration directly, not via weight loss. DR + IGF2 reduces regeneration rate significantly.

#Create a "dose" varibale and look to see if it interacts as a continuous variable.

library(dplyr)
library(tidyr)

## 
## Attaching package: 'tidyr'

## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack

data <- read.csv("R.analysis.currated.csv", stringsAsFactors = FALSE)

data <- data %>%
  group_by(Animal_ID) %>%
  fill(Dose, SVL, .direction = "downup") %>%
  ungroup()


data$Animal_ID <- as.factor(data$Animal_ID)
data$Treatment <- as.factor(data$Treatment)
data$Week <- as.numeric(data$Week)

model_correct <- lmer(
  Rate_Reg ~ Week * Treatment + Week:Dose + (1 | Animal_ID),
  data = data
)

## boundary (singular) fit: see help('isSingular')

summary(model_correct)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Rate_Reg ~ Week * Treatment + Week:Dose + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 1877.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3200 -0.6992 -0.1925  0.5281  5.9955 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 0.000    0.000   
##  Residual              7.578    2.753   
## Number of obs: 385, groups:  Animal_ID, 66
## 
## Fixed effects:
##                        Estimate Std. Error t value
## (Intercept)            3.896711   0.703390   5.540
## Week                  -0.121489   0.139831  -0.869
## TreatmentDR           -0.584570   1.036264  -0.564
## TreatmentDR.IGF1      -0.922248   1.026528  -0.898
## TreatmentDR.IGF2      -0.771645   1.009424  -0.764
## Week:TreatmentDR      -0.019441   0.211571  -0.092
## Week:TreatmentDR.IGF1  0.035191   0.273946   0.128
## Week:TreatmentDR.IGF2  0.016777   0.276434   0.061
## Week:Dose              0.002969   0.086749   0.034
## 
## Correlation of Fixed Effects:
##             (Intr) Week   TrtmDR TDR.IGF1 TDR.IGF2 Wk:TDR W:TDR.IGF1 W:TDR.IGF2
## Week        -0.922                                                             
## TreatmentDR -0.679  0.626                                                      
## TrtmDR.IGF1 -0.685  0.632  0.465                                               
## TrtmDR.IGF2 -0.697  0.642  0.473  0.477                                        
## Wk:TrtmntDR  0.609 -0.661 -0.921 -0.417   -0.425                               
## Wk:TDR.IGF1  0.471 -0.510 -0.351 -0.677   -0.330    0.418                      
## Wk:TDR.IGF2  0.466 -0.506 -0.349 -0.317   -0.665    0.418  0.730               
## Week:Dose    0.000  0.000  0.046 -0.003    0.004   -0.120 -0.676     -0.698    
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

library(ggplot2)

ggplot(data, aes(x = Week, y = Rate_Reg, color = Treatment)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun = mean, geom = "point") +
  theme_minimal()

## Warning: Removed 242 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 242 rows containing non-finite outside the scale range
## (`stat_summary()`).

Result: There is no evidence that treatment, dose, or time affects regeneration rate in this dataset.Regeneration rate is statistically improper as it likely normalizes time effects, that is why week has a base 0 value all the time.

model5 <- lmer(
  Regeneration ~ Week * Treatment + Week:Dose + (1 | Animal_ID),
  data = data
)
summary(model5)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Regeneration ~ Week * Treatment + Week:Dose + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2342.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.69957 -0.61449  0.05928  0.59938  2.96223 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 6.484    2.546   
##  Residual              8.086    2.844   
## Number of obs: 451, groups:  Animal_ID, 67
## 
## Fixed effects:
##                         Estimate Std. Error t value
## (Intercept)           -3.8252866  0.8139124  -4.700
## Week                   3.3667449  0.1219617  27.605
## TreatmentDR            1.0040461  1.1922796   0.842
## TreatmentDR.IGF1       0.8708884  1.1754040   0.741
## TreatmentDR.IGF2       0.5236999  1.1798358   0.444
## Week:TreatmentDR      -0.6768179  0.1857158  -3.644
## Week:TreatmentDR.IGF1 -0.7163933  0.3832053  -1.869
## Week:TreatmentDR.IGF2 -0.5337399  0.3916686  -1.363
## Week:Dose              0.0001053  0.1590532   0.001
## 
## Correlation of Fixed Effects:
##             (Intr) Week   TrtmDR TDR.IGF1 TDR.IGF2 Wk:TDR W:TDR.IGF1 W:TDR.IGF2
## Week        -0.589                                                             
## TreatmentDR -0.683  0.402                                                      
## TrtmDR.IGF1 -0.692  0.408  0.473                                               
## TrtmDR.IGF2 -0.690  0.406  0.471  0.478                                        
## Wk:TrtmntDR  0.387 -0.657 -0.591 -0.267   -0.267                               
## Wk:TDR.IGF1  0.187 -0.318 -0.153 -0.265   -0.130    0.363                      
## Wk:TDR.IGF2  0.183 -0.311 -0.150 -0.125   -0.258    0.360  0.899               
## Week:Dose    0.000  0.000  0.028 -0.003    0.001   -0.173 -0.889     -0.900

ggplot(data, aes(x = Week, y = Regeneration, color = Treatment)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun = mean, geom = "point") +
  theme_minimal()

## Warning: Removed 176 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 176 rows containing non-finite outside the scale range
## (`stat_summary()`).

model6 <- lmer(
  Perc.Regeneration ~ Week * Treatment + Week:Dose + (1 | Animal_ID),
  data = data
)
summary(model6)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Perc.Regeneration ~ Week * Treatment + Week:Dose + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2755.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9316 -0.5501  0.0606  0.5840  3.6280 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 17.74    4.212   
##  Residual              20.08    4.481   
## Number of obs: 452, groups:  Animal_ID, 67
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)            -6.1212     1.3177  -4.645
## Week                    5.4003     0.1924  28.075
## TreatmentDR             1.8247     1.9299   0.945
## TreatmentDR.IGF1        1.8224     1.8974   0.960
## TreatmentDR.IGF2        1.3275     1.9107   0.695
## Week:TreatmentDR       -1.3341     0.2930  -4.554
## Week:TreatmentDR.IGF1  -1.9379     0.6106  -3.174
## Week:TreatmentDR.IGF2  -1.6339     0.6245  -2.617
## Week:Dose               0.1710     0.2542   0.673
## 
## Correlation of Fixed Effects:
##             (Intr) Week   TrtmDR TDR.IGF1 TDR.IGF2 Wk:TDR W:TDR.IGF1 W:TDR.IGF2
## Week        -0.573                                                             
## TreatmentDR -0.683  0.391                                                      
## TrtmDR.IGF1 -0.695  0.398  0.474                                               
## TrtmDR.IGF2 -0.690  0.395  0.471  0.479                                        
## Wk:TrtmntDR  0.376 -0.657 -0.575 -0.261   -0.260                               
## Wk:TDR.IGF1  0.181 -0.315 -0.148 -0.254   -0.126    0.363                      
## Wk:TDR.IGF2  0.177 -0.308 -0.145 -0.121   -0.248    0.360  0.902               
## Week:Dose    0.000  0.000  0.027 -0.002    0.001   -0.175 -0.892     -0.902

Result:

When using regeneration: Regeneration strongly increases over time. No baseline treatment differences (without taking time into account). For the interaction, compared to control, DR slowed regneration, IGF1 has a borderline effect, and IGF2 is non-significant (effect is 0.00 with t= 0.001). So. 1. Regeneration is time drive 2. DR suppresses regeneration 3. IGF treatment does not rescue strongly

A linear mixed-effects model revealed a strong effect of time on regeneration (β = 3.37, t = 27.6), indicating consistent increases in regeneration across weeks. While no baseline differences among treatments were observed, a significant interaction between week and dietary restriction (DR) demonstrated that DR significantly slowed regeneration over time (β = -0.68, t = -3.64). The DR + IGF1 group showed a similar but weaker trend, while DR + IGF2 did not significantly differ from controls. No effect of dose on regeneration trajectory was detected. These results indicate that dietary restriction suppresses regeneration, with limited evidence that IGF treatments rescue this effect.

When using Percent regeneration:

1. Regeneration is time driven
2. ALL treatments decrease regeneration significantly

At this point: When using raw regeneration, DR decreases compared to ad lib, and IGF1/IGF2 are not different from adlib. Which may suggest rescue. But we need a direct DR/IGF1 DR/IGF2 comparison to state that.

library(emmeans)

em <- emtrends(model5, ~ Treatment, var = "Week")

## Registered S3 methods overwritten by 'broom':
##   method        from 
##   nobs.fitdistr MuMIn
##   nobs.multinom MuMIn

pairs(em)

##  contrast          estimate    SE  df t.ratio p.value
##  AdLib - DR          0.6768 0.186 398   3.641  0.0017
##  AdLib - DR.IGF1     0.7164 0.386 373   1.856  0.2491
##  AdLib - DR.IGF2     0.5337 0.395 370   1.353  0.5300
##  DR - DR.IGF1        0.0396 0.363 382   0.109  0.9995
##  DR - DR.IGF2       -0.1431 0.371 378  -0.386  0.9804
##  DR.IGF1 - DR.IGF2  -0.1827 0.174 394  -1.047  0.7218
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

library(ggplot2)
library(dplyr)

ggplot(data, aes(x = Week, y = Regeneration, color = Treatment, group = Treatment)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1.2) +
  stat_summary(fun = mean, geom = "point", size = 2.8) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.15, linewidth = 0.6) +
  theme_classic() +
  labs(
    x = "Week",
    y = "Raw Regeneration",
    color = "Treatment"
  )

## Warning: Removed 176 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 176 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 176 rows containing non-finite outside the scale range
## (`stat_summary()`).

em <- emtrends(model6, ~ Treatment, var = "Week")
pairs(em)

##  contrast          estimate    SE  df t.ratio p.value
##  AdLib - DR           1.334 0.293 398   4.550 <0.0001
##  AdLib - DR.IGF1      1.938 0.615 385   3.151  0.0095
##  AdLib - DR.IGF2      1.634 0.629 382   2.598  0.0478
##  DR - DR.IGF1         0.604 0.577 392   1.046  0.7227
##  DR - DR.IGF2         0.300 0.591 389   0.508  0.9573
##  DR.IGF1 - DR.IGF2   -0.304 0.274 395  -1.108  0.6848
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

DR regenerates significantly slower than Ad lib. IGF groups are not detectably different from Ad lib. IGF groups are not detectably different from DR. DR significantly impairs regeneration, and IGF treatment does not produce a statistically detectable rescue effect.

Pairwise comparisons of regeneration slopes revealed that dietary restriction significantly reduced regeneration relative to ad lib controls (p = 0.0017). However, neither IGF-treated group differed significantly from either ad lib or dietary restriction groups, indicating no statistical evidence that IGF treatment rescues or modifies the effect of dietary restriction on regeneration.

Key takeaway:

Regeneration Percentage = All had stat. sig. negative effect. Raw Regeneration = DR Had negative effect, and IGF1 and IGF2 trended but showed no statistical sign of rescue. IGF not statistically worse than control or better than DR.

#Diet as a continuous variable
library(dplyr)

data <- data %>%
  mutate(
    Diet = case_when(
      Treatment == "AdLib" ~ "AdLib",
      Treatment %in% c("DR", "DR.IGF1", "DR.IGF2") ~ "Restricted"
    )
  )
data$Diet <- factor(data$Diet, levels = c("AdLib", "Restricted"))

model_weightloss <- lmer(
  Perc.Regeneration ~ Week * Perc.WeightLoss + (1 | Animal_ID),
  data = data
)
summary(model_weightloss)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Perc.Regeneration ~ Week * Perc.WeightLoss + (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2882.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1479 -0.5381  0.0486  0.5775  3.3623 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 18.41    4.291   
##  Residual              21.28    4.614   
## Number of obs: 465, groups:  Animal_ID, 69
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          -6.07577    1.13441  -5.356
## Week                  4.85717    0.19194  25.305
## Perc.WeightLoss       0.07376    0.05308   1.390
## Week:Perc.WeightLoss -0.03017    0.01002  -3.011
## 
## Correlation of Fixed Effects:
##             (Intr) Week   Prc.WL
## Week        -0.676              
## Prc.WghtLss -0.792  0.638       
## Wk:Prc.WghL  0.575 -0.852 -0.747

data <- data %>%
  mutate(WL_group = cut(Perc.WeightLoss, breaks = 3))

ggplot(data, aes(x = Week, y = Perc.Regeneration, color = WL_group)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1.2) +
  stat_summary(fun = mean, geom = "point") +
  theme_classic()

## Warning: Removed 175 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 175 rows containing non-finite outside the scale range
## (`stat_summary()`).

Does the degree of restriction (weight loss) affect how regeneration changes over time? As weight loss increases, the rate of regeneration over time decreases.

Low weight loss → faster regeneration High weight loss → slower regeneration

This is a diet dose-dependent suppression effect

A linear mixed-effects model revealed a strong effect of time on percent regeneration (β = 4.86, t = 25.3), indicating consistent increases across weeks. While percent weight loss did not affect baseline regeneration, a significant interaction between week and weight loss (β = −0.030, t = −3.01) demonstrated that increasing weight loss was associated with a reduced rate of regeneration over time. These results indicate that regeneration is sensitive to energetic state, with greater nutritional deficit leading to slower regenerative growth.

data$IGF <- ifelse(data$Treatment %in% c("DR.IGF1", "DR.IGF2"), "IGF", "No_IGF")
data$IGF <- factor(data$IGF)

library(lme4)

model_full <- lmer(
  Perc.Regeneration ~ 
    Week * Perc.WeightLoss +     # energetic effect
    Week * IGF +                 # IGF effect
    Week:Dose +                  # dose effect (only meaningful in IGF)
    (1 | Animal_ID),
  data = data
)

summary(model_full)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Perc.Regeneration ~ Week * Perc.WeightLoss + Week * IGF + Week:Dose +  
##     (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2789.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0001 -0.5795  0.0481  0.5564  3.1733 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 17.85    4.226   
##  Residual              20.88    4.569   
## Number of obs: 452, groups:  Animal_ID, 67
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          -5.45807    1.40959  -3.872
## Week                  4.36686    0.57964   7.534
## Perc.WeightLoss       0.05403    0.05389   1.003
## IGFNo_IGF            -0.70074    1.39854  -0.501
## Week:Perc.WeightLoss -0.02591    0.01015  -2.554
## Week:IGFNo_IGF        0.85035    0.56702   1.500
## Week:Dose             0.03336    0.25399   0.131
## 
## Correlation of Fixed Effects:
##             (Intr) Week   Prc.WL IGFN_I W:P.WL W:IGFN
## Week        -0.222                                   
## Prc.WghtLss -0.717  0.180                            
## IGFNo_IGF   -0.599  0.108  0.153                     
## Wk:Prc.WghL  0.519 -0.198 -0.750 -0.102              
## Wk:IGFN_IGF  0.098 -0.934  0.003 -0.191 -0.059       
## Week:Dose   -0.041 -0.922  0.061  0.022 -0.113  0.933

model_rescue <- lmer(
  Perc.Regeneration ~ 
    Week * Perc.WeightLoss * IGF + 
    Week:Dose +
    (1 | Animal_ID),
  data = data
)

summary(model_rescue)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Perc.Regeneration ~ Week * Perc.WeightLoss * IGF + Week:Dose +  
##     (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2781.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9346 -0.5802  0.0584  0.5439  3.2222 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 17.31    4.16    
##  Residual              20.16    4.49    
## Number of obs: 452, groups:  Animal_ID, 67
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    -2.06192    2.11169  -0.976
## Week                            3.32697    0.62605   5.314
## Perc.WeightLoss                -0.13838    0.10010  -1.382
## IGFNo_IGF                      -5.13420    2.49455  -2.058
## Week:Perc.WeightLoss            0.03625    0.01851   1.958
## Week:IGFNo_IGF                  2.27937    0.65995   3.454
## Perc.WeightLoss:IGFNo_IGF       0.26426    0.11782   2.243
## Week:Dose                       0.01792    0.24989   0.072
## Week:Perc.WeightLoss:IGFNo_IGF -0.08783    0.02188  -4.014
## 
## Correlation of Fixed Effects:
##             (Intr) Week   Prc.WL IGFN_I Wk:P.WL W:IGFN P.WL:I Wek:Ds
## Week        -0.365                                                  
## Prc.WghtLss -0.889  0.357                                           
## IGFNo_IGF   -0.846  0.311  0.752                                    
## Wk:Prc.WghL  0.665 -0.446 -0.772 -0.562                             
## Wk:IGFN_IGF  0.347 -0.939 -0.339 -0.412  0.424                      
## P.WL:IGFN_I  0.754 -0.317 -0.849 -0.834  0.655   0.399              
## Week:Dose   -0.056 -0.828  0.064  0.044 -0.082   0.775 -0.038       
## W:P.WL:IGFN -0.560  0.414  0.650  0.619 -0.842  -0.535 -0.762  0.025

ggplot(data, aes(x = Perc.WeightLoss, y = Perc.Regeneration, color = IGF)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ Week) +
  theme_classic()

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

## Warning: Removed 175 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 175 rows containing missing values or values outside the scale range
## (`geom_point()`).

Model Rescue: Does IGF change how weight loss affects regeneration over time? Does the effect of weight loss differ between IGF and non-IGF animals?

1. Regeneration increases over time
2. In No IGF animals: Weight loss strongly reduces regeneration over time
3. In IGF animals: That negative effect is much weaker (or absent)

Weight loss harms regeneration, but IGF reduces that harm. IGF does not fully restore regeneration to control levels, but it buffers the negative effect of energetic deficit.

A significant three-way interaction between week, percent weight loss, and IGF treatment (β = −0.088, t = −4.01) indicates that IGF modifies the relationship between energetic deficit and regeneration. Specifically, while increased weight loss is associated with reduced regeneration over time in non-IGF animals, this negative effect is attenuated in IGF-treated animals. These results suggest that IGF partially buffers the impact of energetic limitation on regeneration.

IMPORTANT: Same pattern is seen with raw regeneration when using Precent.WeightLoss instead of “Diet”

library(dplyr)

data <- data %>%
  mutate(
    IGF_type = case_when(
      Treatment %in% c("AdLib", "DR") ~ "No_IGF",
      Treatment == "DR.IGF1" ~ "IGF1",
      Treatment == "DR.IGF2" ~ "IGF2"
    )
  )

data$IGF_type <- factor(data$IGF_type, levels = c("No_IGF", "IGF1", "IGF2"))

library(lme4)

model_igf_split <- lmer(
  Perc.Regeneration ~ Week * Perc.WeightLoss * IGF_type + Week:Dose + (1 | Animal_ID),
  data = data
)

summary(model_igf_split)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Perc.Regeneration ~ Week * Perc.WeightLoss * IGF_type + Week:Dose +  
##     (1 | Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2772
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9805 -0.5759  0.0440  0.5106  3.2589 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 17.34    4.164   
##  Residual              19.65    4.433   
## Number of obs: 452, groups:  Animal_ID, 67
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                       -7.20039    1.31779  -5.464
## Week                               5.60541    0.22344  25.087
## Perc.WeightLoss                    0.12611    0.06162   2.047
## IGF_typeIGF1                       4.38031    2.92191   1.499
## IGF_typeIGF2                       3.14519    3.94815   0.797
## Week:Perc.WeightLoss              -0.05149    0.01164  -4.422
## Week:IGF_typeIGF1                 -1.54634    0.71797  -2.154
## Week:IGF_typeIGF2                 -2.72059    0.80250  -3.390
## Perc.WeightLoss:IGF_typeIGF1      -0.20687    0.14575  -1.419
## Perc.WeightLoss:IGF_typeIGF2      -0.18734    0.17827  -1.051
## Week:Dose                          0.01291    0.24820   0.052
## Week:Perc.WeightLoss:IGF_typeIGF1  0.03403    0.02943   1.156
## Week:Perc.WeightLoss:IGF_typeIGF2  0.11873    0.02981   3.983

## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

ggplot(data, aes(x = Perc.WeightLoss, y = Perc.Regeneration, color = IGF_type)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ Week) +
  theme_classic()

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

## Warning: Removed 175 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 175 rows containing missing values or values outside the scale range
## (`geom_point()`).

library(ggeffects)
library(ggplot2)

pred <- ggpredict(
  model_igf_split,
  terms = c("Perc.WeightLoss", "IGF_type", "Week [2,4,6,8]")
)

pred2 <- ggpredict(
  model_igf_split,
  terms = c("Perc.WeightLoss", "IGF_type", "Week [2,4]")
)

ggplot(pred, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  facet_wrap(~ facet) +
  theme_classic() +
  labs(
    x = "Percent weight loss",
    y = "Predicted percent regeneration",
    color = "IGF treatment",
    fill = "IGF treatment"
  )

range(data$Perc.WeightLoss, na.rm = TRUE)

## [1] -55.11265  49.84560

pred <- ggpredict(
  model_igf_split,
  terms = c("Perc.WeightLoss [-55.2:48.85]", "IGF_type", "Week [2,4,6,8]")
)

ggplot(pred, aes(x = x, y = predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.2) +
  facet_wrap(~ facet) +
  theme_classic() +
  labs(
    x = "Percent weight loss",
    y = "Predicted percent regeneration",
    color = "IGF treatment",
    fill = "IGF treatment"
  )

library(emmeans)

em <- emtrends(
  model_igf_split,
  ~ IGF_type | Week,
  var = "Perc.WeightLoss",
  at = list(Week = c(2,4, 6, 8))
)

pairs(em)

## Week = 2:
##  contrast      estimate     SE  df t.ratio p.value
##  No_IGF - IGF1  0.13882 0.1070 400   1.296  0.3985
##  No_IGF - IGF2 -0.05011 0.1370 386  -0.365  0.9293
##  IGF1 - IGF2   -0.18893 0.1610 382  -1.175  0.4689
## 
## Week = 4:
##  contrast      estimate     SE  df t.ratio p.value
##  No_IGF - IGF1  0.07076 0.0923 408   0.767  0.7237
##  No_IGF - IGF2 -0.28757 0.1130 333  -2.542  0.0307
##  IGF1 - IGF2   -0.35833 0.1330 359  -2.685  0.0207
## 
## Week = 6:
##  contrast      estimate     SE  df t.ratio p.value
##  No_IGF - IGF1  0.00271 0.1120 439   0.024  0.9997
##  No_IGF - IGF2 -0.52502 0.1180 388  -4.464 <0.0001
##  IGF1 - IGF2   -0.52773 0.1470 422  -3.584  0.0011
## 
## Week = 8:
##  contrast      estimate     SE  df t.ratio p.value
##  No_IGF - IGF1 -0.06535 0.1530 434  -0.426  0.9047
##  No_IGF - IGF2 -0.76247 0.1480 439  -5.142 <0.0001
##  IGF1 - IGF2   -0.69713 0.1940 438  -3.601  0.0010
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

em_df <- as.data.frame(em)

Result: Does IGF1 vs IGF2 independently buffer the effect of weight loss?

IGF1 does not significantly change the effect of weight loss. No clear evidence of rescue. IGF2 is Strongly significant and positive. IGF2 reduces the negative effect of weight loss, evidence of rescue.

data$IGF <- ifelse(data$Treatment %in% c("DR.IGF1", "DR.IGF2"), "IGF", "No_IGF")
data$IGF <- factor(data$IGF)

library(lme4)

model_rescue <- lmer(
  Regeneration ~ 
    Week * Perc.WeightLoss * IGF + 
    Week:Dose +
    (1 | Animal_ID),
  data = data
)

summary(model_rescue)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Regeneration ~ Week * Perc.WeightLoss * IGF + Week:Dose + (1 |  
##     Animal_ID)
##    Data: data
## 
## REML criterion at convergence: 2374.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.87777 -0.56843  0.05979  0.60292  2.89221 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Animal_ID (Intercept) 6.619    2.573   
##  Residual              8.200    2.864   
## Number of obs: 451, groups:  Animal_ID, 67
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    -1.26923    1.33750  -0.949
## Week                            2.51286    0.39681   6.333
## Perc.WeightLoss                -0.10342    0.06353  -1.628
## IGFNo_IGF                      -2.70142    1.57860  -1.711
## Week:Perc.WeightLoss            0.02648    0.01179   2.247
## Week:IGFNo_IGF                  0.80765    0.41860   1.929
## Perc.WeightLoss:IGFNo_IGF       0.14182    0.07482   1.895
## Week:Dose                      -0.10391    0.15792  -0.658
## Week:Perc.WeightLoss:IGFNo_IGF -0.04254    0.01394  -3.053
## 
## Correlation of Fixed Effects:
##             (Intr) Week   Prc.WL IGFN_I Wk:P.WL W:IGFN P.WL:I Wek:Ds
## Week        -0.371                                                  
## Prc.WghtLss -0.890  0.362                                           
## IGFNo_IGF   -0.847  0.317  0.754                                    
## Wk:Prc.WghL  0.667 -0.450 -0.773 -0.565                             
## Wk:IGFN_IGF  0.352 -0.939 -0.344 -0.419  0.427                      
## P.WL:IGFN_I  0.755 -0.321 -0.848 -0.837  0.655   0.404              
## Week:Dose   -0.054 -0.826  0.062  0.042 -0.081   0.772 -0.036       
## W:P.WL:IGFN -0.562  0.417  0.651  0.622 -0.842  -0.538 -0.764  0.024

ggplot(data, aes(x = Perc.WeightLoss, y = Perc.Regeneration, color = IGF)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ Week) +
  theme_classic()

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

## Warning: Removed 175 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 175 rows containing missing values or values outside the scale range
## (`geom_point()`).

library(ggplot2)

ggplot(em_df, aes(x = factor(Week), y = Perc.WeightLoss.trend, color = IGF_type)) +
  geom_point(position = position_dodge(width = 0.4), size = 3) +
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.4),
    width = 0.2
  ) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_classic() +
  labs(
    x = "Week",
    y = "Effect of weight loss on regeneration (slope)",
    color = "IGF type"
  )

library(car)

library(car)

model_data <- model.frame(model_igf_split)

model_data$resid <- residuals(model_igf_split)

leveneTest(resid ~ IGF_type, data = model_data)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   2  2.8544 0.05864 .
##       449                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_data_igf1 <- model_data %>%
  filter(IGF_type %in% c("No_IGF", "IGF1"))

leveneTest(resid ~ IGF_type, data = model_data_igf1)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   1   4.759 0.02985 *
##       331                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_data_igf2 <- model_data %>%
  filter(IGF_type %in% c("No_IGF", "IGF2"))

leveneTest(resid ~ IGF_type, data = model_data_igf2)

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1   2.377 0.1241
##       343

The authors would have to shown negative covariance between weight change and tail regeneration to show a change in prioritization of energy use.

library(dplyr)

final_data <- data %>%
  group_by(Animal_ID) %>%
  summarize(
    final_regeneration = max(Regeneration, na.rm = TRUE),
    weight_loss = first(Perc.WeightLoss),
    Treatment = first(Treatment)
  )

## Warning: There were 11 warnings in `summarize()`.
## The first warning was:
## ℹ In argument: `final_regeneration = max(Regeneration, na.rm = TRUE)`.
## ℹ In group 15: `Animal_ID = BC234`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 10 remaining warnings.

cor.test(final_data$weight_loss, final_data$final_regeneration)

## 
##  Pearson's product-moment correlation
## 
## data:  final_data$weight_loss and final_data$final_regeneration
## t = 0.41102, df = 65, p-value = 0.6824
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1916378  0.2876057
## sample estimates:
##        cor 
## 0.05091453

lm_tradeoff_igf <- lm(final_regeneration ~ weight_loss * Treatment, data = final_data)
summary(lm_tradeoff_igf)

## 
## Call:
## lm(formula = final_regeneration ~ weight_loss * Treatment, data = final_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.2959  -3.8851   0.0258   4.0045  13.7746 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  17.65763    3.44728   5.122 3.49e-06 ***
## weight_loss                   0.03591    0.21863   0.164    0.870    
## TreatmentDR                  -5.07518    4.34895  -1.167    0.248    
## TreatmentDR.IGF1              1.28540    5.20592   0.247    0.806    
## TreatmentDR.IGF2             -9.61449    6.63895  -1.448    0.153    
## weight_loss:TreatmentDR       0.10406    0.24833   0.419    0.677    
## weight_loss:TreatmentDR.IGF1 -0.35296    0.28919  -1.220    0.227    
## weight_loss:TreatmentDR.IGF2  0.38573    0.33482   1.152    0.254    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.407 on 59 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1881, Adjusted R-squared:  0.09182 
## F-statistic: 1.953 on 7 and 59 DF,  p-value: 0.07706

Result: We tested for a direct covariance between weight loss and regeneration at the individual level, but this relationship was not statistically significant. However, a mixed-effects model revealed a significant interaction between time and weight loss, indicating that individuals experiencing greater energetic deficit exhibited slower regeneration over time. These results suggest that energetic state influences regeneration dynamics, even if a simple endpoint tradeoff is not detectable.

R2R: While a simple covariance between weight loss and final regeneration was not detected, such analyses collapse dynamic processes into a single endpoint and may obscure temporal relationships. Our longitudinal analysis demonstrates that weight loss significantly alters regeneration trajectories, providing evidence that energetic state constrains regenerative processes over time.