Code from: Genetic and environmental factors shape rates of plasticity: the temporal dynamics of opsin gene expression in aquatic environments
César Bertinetti
March 2025
The datasets needed to reproduce this analysis are deposited on Bertinetti et al. 2025. Data from: Genetic and environmental factors shape rates of plasticity: the temporal dynamics of opsin gene expression in aquatic environments [Dataset]. Zenodo.
Following packages are required:
R.version$platform## [1] "aarch64-apple-darwin20"R.version$version.string## [1] "R version 4.4.1 (2024-06-14)"installed.packages()[names(sessionInfo()$otherPkgs), "Version"]## minpack.lm nls2 proto ggtext ragg colorspace
## "1.2-4" "0.3-4" "1.0.0" "0.1.2" "1.3.3" "2.1-1"
## ggdist ggpubr vegan lattice permute ggfortify
## "3.3.2" "0.6.0" "2.6-8" "0.22-6" "0.9-7" "0.4.17"
## lubridate forcats purrr tibble tidyverse broom
## "1.9.3" "1.0.0" "1.0.2" "3.2.1" "2.0.0" "1.0.7"
## readxl devtools usethis effectsize partR2 rsq
## "1.4.3" "2.4.5" "3.0.0" "0.8.9" "0.9.2" "2.6"
## boot lme4 Matrix performance sjstats visreg
## "1.3-31" "1.1-35.5" "1.7-0" "0.12.3" "0.19.0" "2.7.0"
## pracma gridExtra ggplot2 nlme agricolae corrplot
## "2.4.4" "2.3" "3.5.1" "3.1-166" "1.3-7" "0.95"
## Hmisc car carData zoo MESS MASS
## "5.2-1" "3.1-3" "3.0-5" "1.8-12" "0.5.12" "7.3-61"
## RColorBrewer plotrix tidyr dplyr matrixStats stringr
## "1.1-3" "3.8-4" "1.3.1" "1.1.4" "1.4.1" "1.5.1"
## readr
## "2.1.5"Read Opsin Gene Expression and Plot Distribution
The “Opsin gene expression.csv” contains the raw opsin gene expression values for all individuals
## READ MASTERTABLE GENERATED FROM COMBINE qPCR
opsin <- read.table("Opsin gene expression.csv", header=T, sep =",", dec=".",stringsAsFactors = T) [,-1]
# Combine controls - this part can be left out for further resolutoin before 20 dph
for (i in 1:nrow(opsin)) {
if(opsin$Light[i] == "WR" & opsin$Age[i] < 21){
opsin$Light[i] <- "WW"
}else if (opsin$Light[i] == "RW" & opsin$Age[i] < 21){
opsin$Light[i] <- "RR"
} else {
opsin$Light[i] <- opsin$Light[i]
}
}
## SPLIT TO CREATE "BROAD GROUPS" THAT HELP WHEN PlOTTING (ALL TREATMENTS)
opsin <- opsin %>%
unite("onto", c("Stage","Age"), sep="", remove=F) %>%
unite("combi", c("Light", "onto","Spec"), sep="-", remove=F)
opsin <- opsin[,-5] # removes onto column
## ADD DAYS AFTER CHANGE AS A VARIABLE
for (i in 1:nrow(opsin)) {
if(opsin$Stage[i] == "E" & opsin$Age[i] >= 20){
opsin$day[i] <- as.numeric(opsin$Age[i] - 20)
}else if (opsin$Stage[i] == "L"){
opsin$day[i] <- as.numeric(opsin$Age[i] - 35)
} else {
opsin$day[i] <- NA
}
}
opsin1 <- opsin[,c(3,8:18)] ### EXTRACT IMPORTANT VARIABLES BEFORE AGGREGGATING
## FIRST MAKES SENSE TO AGGREGATE REPLICATES
## THIS IS A SIMPLE FUNCTION TO GET M + SE + SAMPLE SIZE
momred = function(x) {
c(m = mean(x), n = length(x), se = std.error(x))
}
## HERE ONLY FOR SC NOW
mean_exp <- aggregate(cbind(SC_sws1,SC_sws2b,SC_sws2a,PSI_SC) ~ combi, data = opsin1, momred)
mean_exp <- do.call(data.frame, mean_exp)
opsin_mean <- unique(merge( opsin[,c(3:7,19)],mean_exp, by="combi")) ## MERGE TO RETTACH METADATA
opsin_mean <- opsin_mean %>% ## SPLIT TO CREATE "BROAD GROUPS" THAT HELP WHEN PlOTTING (ALL TREATMENTS)
unite("brgroups", c("Light", "Stage","Spec"), sep="-", remove=F)## SUBSETTIN FOR ONLY EARLY TREATMENT
sp <- "XIL" # "NIC"
st <- "E"
#opsin_early <- subset(opsin_mean, Stage == st )
opsin_early <- subset(opsin_mean, Spec == sp & Stage == st & Age < 51)
#genes <- c("sws1","sws2b","sws2a","rh2b","rh2a","lws")
genes <- c("SC_sws1","SC_sws2b","SC_sws2a") ## SAME BUT PROP OF SC OR DC
P_plots <- list()
for (i in 1:length(genes)) {
b <- as.name(paste(as.character(genes[i]),".m", sep="")) ## creates column name for variable "mean"
c <- as.name(paste(as.character(genes[i]),".se", sep=""))
if (i==7) {
pl <- ggplot(opsin_early, aes(x=Age, y=eval(as.name(b)), group=brgroups, color=brgroups)) +
geom_line()+
theme_classic()+
geom_pointrange(aes(ymin=eval(as.name(b))-eval(as.name(c)), ymax=eval(as.name(b))+eval(as.name(c))))+
ylab(as.name(b))+
scale_x_continuous(breaks=seq(0,50,5)) +
geom_vline(xintercept=c(20,35), linetype=2) +
theme(legend.position = "inside",legend.position.inside=c(0.8,0.8)) +
ylim(0,1)
P_plots[[i]] <- ggplotGrob (pl) ## save to be able to plot outside of loop
} else {
pl <- ggplot(opsin_early, aes(x=Age, y=eval(as.name(b)), group=brgroups, color=brgroups)) +
geom_line()+
theme_classic()+
geom_pointrange(aes(ymin=eval(as.name(b))-eval(as.name(c)), ymax=eval(as.name(b))+eval(as.name(c))))+
ylab(as.name(b))+
scale_x_continuous(breaks=seq(0,50,5)) +
geom_vline(xintercept=c(20,35), linetype=2) +
ylim(0,1)+
theme(legend.position = "none") ## don't want legend in all plots
P_plots[[i]] <- ggplotGrob (pl)
}
# P_plots[[i]] <- ggplotGrob (pl)
#print(P_plots[[i]])
}grid.arrange(grobs=P_plots, ncol=3, top = textGrob(sp)) 
Read and visualize single cone opsin predicted sensitivity index
P_plots <- list()
opsin_m <- subset(opsin_mean, Spec == "XIL" & Age < 51)
xil <- ggplot(opsin_m, aes(x=Age, y=PSI_SC.m, group=Light, color=Light)) +
geom_line()+
theme_classic()+
#geom_pointrange(aes(ymin=eval(as.name(b))-eval(as.name(c)), ymax=eval(as.name(b))+eval(as.name(c))))+
geom_pointrange(aes(ymin=PSI_SC.m-PSI_SC.se, ymax=PSI_SC.m+PSI_SC.se))+
labs(x = " Age (days post-hatch)", y = "Single Cone Predicted Sensitivity Index (nm)", title="Xiloá (Clear Lake)")+
scale_x_continuous(breaks=seq(0,50,5)) +
geom_vline(xintercept=c(20,35), linetype=2) +
theme(legend.position = "none", text= element_text(size=12)) +
ylim(350,500)
P_plots[[1]] <- ggplotGrob (xil) ## save to be able to plot outside of loop
opsin_m <- subset(opsin_mean, Spec == "NIC")
nic <- ggplot(opsin_m, aes(x=Age, y=PSI_SC.m, group=Light, color=Light)) +
geom_line()+
theme_classic()+
#geom_pointrange(aes(ymin=eval(as.name(b))-eval(as.name(c)), ymax=eval(as.name(b))+eval(as.name(c))))+
geom_pointrange(aes(ymin=PSI_SC.m-PSI_SC.se, ymax=PSI_SC.m+PSI_SC.se))+
labs(x = " Age (dph)", y = "Single Cone Predicted Sensitivity Index (nm)", title="Nicaragua (Turbid Lake)")+
scale_x_continuous(breaks=seq(0,50,5)) +
geom_vline(xintercept=c(20,35), linetype=2) +
theme(text= element_text(size=12)) +
ylim(350,500)
P_plots[[2]] <- ggplotGrob (nic)grid.arrange(grobs=P_plots, ncol=2) ## put all panels together
Model mean cone opsin predicted sensitivity index PSI_SC
The “LightData.zip” contains the raw irradiance measurements, the normalized irradiance and the photic factors
df <- read.csv("PlasticRate-PSI_SC.csv", sep = ",", dec = ".")
### Rename variables
df$y <- df$PSI_SC.m
df$t <- df$day
df$lake <-df$Spec
df$age <-df$Stage
df$light <-df$Light
### model_key <- nls(y ~ yf + (y0 - yf) * exp(-lambda * t),
### model_used <- nls(y ~ yf + (Dy0yf) * exp(-lambda * t),
df <- df %>% filter(lake %in% c("NIC", "XIL"))
df <- df %>% filter(light %in% c("RW","WR"))
#df$lake_indicator <- ifelse(df$lake == "nicaragua", 0, 1)
df$lake_indicator <- ifelse(df$lake == "NIC", 0, 1)
df$age_indicator <- ifelse(df$Stage == "E", 0, 1)
df$change_indicator <- ifelse(df$light == "WR", 0, 1)full_model <- nlsLM(y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator + yf_change * change_indicator +
yf_chXlake * change_indicator * lake_indicator + yf_chXage * change_indicator * age_indicator +
yf_lakeXage * lake_indicator * age_indicator + yf_lakeXageXch * change_indicator * lake_indicator * age_indicator) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator + Dy0yf_change * change_indicator +
Dy0yf_chXlake * change_indicator * lake_indicator + Dy0yf_chXage * change_indicator * age_indicator +
Dy0yf_lakeXage * lake_indicator * age_indicator + Dy0yf_lakeXageXch * change_indicator * lake_indicator * age_indicator) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator + lambda_change * change_indicator +
lambda_CHxLa * change_indicator * lake_indicator + lambda_chXage * change_indicator * age_indicator +
lambda_lakeXage * lake_indicator * age_indicator + lambda_lakeXageXch * change_indicator * lake_indicator * age_indicator) * t),
data = df,
start = list(Dy0yf = -75, Dy0yf_lake = -20, Dy0yf_age = -34.88, Dy0yf_change = 100, Dy0yf_chXlake = 30,
Dy0yf_chXage = 0, Dy0yf_lakeXage = 0, Dy0yf_lakeXageXch = 0, yf = 465.72,
yf_lake = 1.22, yf_age = 1.22, yf_change = 1.22, yf_chXlake = 1.22, yf_chXage = 1.22,
yf_lakeXage = 1.22, yf_lakeXageXch = 1.22, lambda = 0.0015, lambda_lake = 0.00015,
lambda_age = 0.00015, lambda_change = 0.00015, lambda_CHxLa = 0.00015, lambda_chXage = 0.00015,
lambda_lakeXage = 0.00015, lambda_lakeXageXch = 0.00015),
control=nls.lm.control(maxiter = 200, factor = 10, maxfev = 10000, nprint=F))
summary(full_model)##
## Formula: y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator +
## yf_change * change_indicator + yf_chXlake * change_indicator *
## lake_indicator + yf_chXage * change_indicator * age_indicator +
## yf_lakeXage * lake_indicator * age_indicator + yf_lakeXageXch *
## change_indicator * lake_indicator * age_indicator) + (Dy0yf +
## Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator +
## Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator *
## lake_indicator + Dy0yf_chXage * change_indicator * age_indicator +
## Dy0yf_lakeXage * lake_indicator * age_indicator + Dy0yf_lakeXageXch *
## change_indicator * lake_indicator * age_indicator) * exp(-(lambda +
## lambda_lake * lake_indicator + lambda_age * age_indicator +
## lambda_change * change_indicator + lambda_CHxLa * change_indicator *
## lake_indicator + lambda_chXage * change_indicator * age_indicator +
## lambda_lakeXage * lake_indicator * age_indicator + lambda_lakeXageXch *
## change_indicator * lake_indicator * age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -79.892288 6.684818 -11.951 < 2e-16 ***
## Dy0yf_lake -12.876218 9.438977 -1.364 0.175551
## Dy0yf_age 16.128288 9.563703 1.686 0.094804 .
## Dy0yf_change 143.571384 8.836102 16.248 < 2e-16 ***
## Dy0yf_chXlake 16.033804 12.551052 1.277 0.204359
## Dy0yf_chXage -79.805248 11.173888 -7.142 1.45e-10 ***
## Dy0yf_lakeXage 2.172256 13.270943 0.164 0.870306
## Dy0yf_lakeXageXch 65.897342 16.784612 3.926 0.000158 ***
## yf 463.808182 2.470875 187.710 < 2e-16 ***
## yf_lake -3.378042 3.193849 -1.058 0.292729
## yf_age -0.052515 3.069879 -0.017 0.986385
## yf_change -62.162665 4.847874 -12.823 < 2e-16 ***
## yf_chXlake -13.504005 6.087481 -2.218 0.028772 *
## yf_chXage 50.786387 5.726370 8.869 2.76e-14 ***
## yf_lakeXage 1.092867 4.355230 0.251 0.802376
## yf_lakeXageXch -46.548948 7.548951 -6.166 1.45e-08 ***
## lambda 0.711666 0.132261 5.381 4.82e-07 ***
## lambda_lake 0.076997 0.181223 0.425 0.671833
## lambda_age 1.046064 0.618977 1.690 0.094115 .
## lambda_change -0.461015 0.145466 -3.169 0.002023 **
## lambda_CHxLa 0.006309 0.203237 0.031 0.975296
## lambda_chXage -1.911756 1.518829 -1.259 0.211039
## lambda_lakeXage -1.325600 0.637952 -2.078 0.040255 *
## lambda_lakeXageXch 2.273226 1.530048 1.486 0.140468
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.61 on 101 degrees of freedom
##
## Number of iterations to convergence: 63
## Achieved convergence tolerance: 1.49e-08## NO yf_lakeXageXch NO Dy0yf_lakeXageXch
afm_v3 <- nlsLM(y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator + yf_change * change_indicator + yf_chXlake * change_indicator* lake_indicator+ yf_chXage * change_indicator* age_indicator + yf_lakeXage * lake_indicator* age_indicator ) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator+ Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator* lake_indicator + Dy0yf_chXage * change_indicator* age_indicator + Dy0yf_lakeXage * lake_indicator* age_indicator) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator+ lambda_change * change_indicator + lambda_CHxLa * change_indicator* lake_indicator + lambda_chXage * change_indicator* age_indicator + lambda_lakeXage * lake_indicator* age_indicator + lambda_lakeXageXch * change_indicator * lake_indicator* age_indicator) * t),
data = df,
start = list(Dy0yf = -75, Dy0yf_lake = -20, Dy0yf_age = -30, Dy0yf_change = 120, Dy0yf_chXlake = 60, Dy0yf_chXage=0, Dy0yf_lakeXage=0, yf =460, yf_lake = 0, yf_age = 0 ,yf_change = -50, yf_chXlake = -50, yf_chXage = 30, yf_lakeXage=10,
lambda = 0.8, lambda_lake = 0, lambda_age = -0.7, lambda_change = -0.4, lambda_CHxLa = -.2, lambda_chXage=1, lambda_lakeXage=0, lambda_lakeXageXch=0),
control=nls.lm.control(maxiter = 100))
summary(afm_v3)##
## Formula: y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator +
## yf_change * change_indicator + yf_chXlake * change_indicator *
## lake_indicator + yf_chXage * change_indicator * age_indicator +
## yf_lakeXage * lake_indicator * age_indicator) + (Dy0yf +
## Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator +
## Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator *
## lake_indicator + Dy0yf_chXage * change_indicator * age_indicator +
## Dy0yf_lakeXage * lake_indicator * age_indicator) * exp(-(lambda +
## lambda_lake * lake_indicator + lambda_age * age_indicator +
## lambda_change * change_indicator + lambda_CHxLa * change_indicator *
## lake_indicator + lambda_chXage * change_indicator * age_indicator +
## lambda_lakeXage * lake_indicator * age_indicator + lambda_lakeXageXch *
## change_indicator * lake_indicator * age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -77.124444 6.234228 -12.371 < 2e-16 ***
## Dy0yf_lake -18.510531 8.095180 -2.287 0.02426 *
## Dy0yf_age 10.286864 8.150115 1.262 0.20974
## Dy0yf_change 139.195009 8.037554 17.318 < 2e-16 ***
## Dy0yf_chXlake 25.433710 9.697962 2.623 0.01005 *
## Dy0yf_chXage -23.692679 9.182137 -2.580 0.01128 *
## Dy0yf_lakeXage 13.376405 9.091356 1.471 0.14425
## yf 463.863921 2.499707 185.567 < 2e-16 ***
## yf_lake -3.467858 3.209343 -1.081 0.28242
## yf_age -0.112985 3.092562 -0.037 0.97093
## yf_change -62.796399 5.113593 -12.280 < 2e-16 ***
## yf_chXlake -12.518753 6.205283 -2.017 0.04625 *
## yf_chXage 3.753314 4.843852 0.775 0.44020
## yf_lakeXage 1.262930 4.380464 0.288 0.77369
## lambda 0.691199 0.130868 5.282 7.17e-07 ***
## lambda_lake 0.116531 0.177007 0.658 0.51179
## lambda_age 1.109737 0.616393 1.800 0.07473 .
## lambda_change -0.454872 0.145495 -3.126 0.00230 **
## lambda_CHxLa -0.001979 0.203239 -0.010 0.99225
## lambda_chXage -1.346733 0.620569 -2.170 0.03229 *
## lambda_lakeXage -1.424007 0.631082 -2.256 0.02615 *
## lambda_lakeXageXch 1.714524 0.647081 2.650 0.00933 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.618 on 103 degrees of freedom
##
## Number of iterations to convergence: 17
## Achieved convergence tolerance: 1.49e-08# noDy0yf_lakeXageXch # noDy0yf_lakeXage #noyf_lakeXage
afm_v4 <- nlsLM(y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator + yf_change * change_indicator + yf_chXlake * change_indicator* lake_indicator+ yf_chXage * change_indicator* age_indicator) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator+ Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator* lake_indicator + Dy0yf_chXage * change_indicator* age_indicator) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator+ lambda_change * change_indicator + lambda_CHxLa * change_indicator* lake_indicator + lambda_chXage * change_indicator* age_indicator + lambda_lakeXage * lake_indicator* age_indicator + lambda_lakeXageXch * change_indicator * lake_indicator* age_indicator) * t),
data = df,
start = list(Dy0yf = -75, Dy0yf_lake = -20, Dy0yf_age = -30, Dy0yf_change = 120, Dy0yf_chXlake = 60, Dy0yf_chXage=0, yf =460, yf_lake = 0, yf_age = 0 ,yf_change = -50, yf_chXlake = -50, yf_chXage = 30,
lambda = 0.8, lambda_lake = 0, lambda_age = -0.7, lambda_change = -0.4, lambda_CHxLa = -.2, lambda_chXage=1, lambda_lakeXage=0, lambda_lakeXageXch=0),
control=nls.lm.control(maxiter = 100))
summary(afm_v4)##
## Formula: y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator +
## yf_change * change_indicator + yf_chXlake * change_indicator *
## lake_indicator + yf_chXage * change_indicator * age_indicator) +
## (Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator +
## Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator *
## lake_indicator + Dy0yf_chXage * change_indicator * age_indicator) *
## exp(-(lambda + lambda_lake * lake_indicator + lambda_age *
## age_indicator + lambda_change * change_indicator +
## lambda_CHxLa * change_indicator * lake_indicator +
## lambda_chXage * change_indicator * age_indicator +
## lambda_lakeXage * lake_indicator * age_indicator +
## lambda_lakeXageXch * change_indicator * lake_indicator *
## age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -80.36585 5.82236 -13.803 < 2e-16 ***
## Dy0yf_lake -11.82028 6.68953 -1.767 0.080137 .
## Dy0yf_age 17.25884 6.68571 2.581 0.011218 *
## Dy0yf_change 139.70527 8.18401 17.071 < 2e-16 ***
## Dy0yf_chXlake 25.74500 9.72490 2.647 0.009363 **
## Dy0yf_chXage -25.73991 8.93712 -2.880 0.004819 **
## yf 463.46081 2.03899 227.299 < 2e-16 ***
## yf_lake -2.78953 2.17925 -1.280 0.203351
## yf_age 0.48885 2.18526 0.224 0.823425
## yf_change -63.62189 5.59746 -11.366 < 2e-16 ***
## yf_chXlake -11.41145 6.23361 -1.831 0.069992 .
## yf_chXage 3.75307 4.22982 0.887 0.376952
## lambda 0.72550 0.12673 5.725 9.92e-08 ***
## lambda_lake 0.05238 0.16376 0.320 0.749705
## lambda_age 1.00074 0.59442 1.684 0.095235 .
## lambda_change -0.51516 0.13933 -3.697 0.000348 ***
## lambda_CHxLa 0.12142 0.18433 0.659 0.511506
## lambda_chXage -1.20820 0.59719 -2.023 0.045599 *
## lambda_lakeXage -1.25962 0.59451 -2.119 0.036470 *
## lambda_lakeXageXch 1.46618 0.60303 2.431 0.016737 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.66 on 105 degrees of freedom
##
## Number of iterations to convergence: 17
## Achieved convergence tolerance: 1.49e-08# noDy0yf_lakeXageXch # noDy0yf_lakeXage #noyf_lakeXage
# no yf_chXage # yf_chXlake
afm_v5 <- nlsLM(y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator + yf_change * change_indicator) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator+ Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator* lake_indicator + Dy0yf_chXage * change_indicator* age_indicator) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator+ lambda_change * change_indicator + lambda_CHxLa * change_indicator* lake_indicator + lambda_chXage * change_indicator* age_indicator + lambda_lakeXage * lake_indicator* age_indicator + lambda_lakeXageXch * change_indicator * lake_indicator* age_indicator) * t),
data = df,
start = list(Dy0yf = -75, Dy0yf_lake = -20, Dy0yf_age = -30, Dy0yf_change = 120, Dy0yf_chXlake = 60, Dy0yf_chXage=0, yf =460, yf_lake = 0, yf_age = 0 ,yf_change = -50,
lambda = 0.8, lambda_lake = 0, lambda_age = -0.7, lambda_change = -0.4, lambda_CHxLa = -.2, lambda_chXage=1, lambda_lakeXage=0, lambda_lakeXageXch=0),
control=nls.lm.control(maxiter = 100))
summary(afm_v5)##
## Formula: y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator +
## yf_change * change_indicator) + (Dy0yf + Dy0yf_lake * lake_indicator +
## Dy0yf_age * age_indicator + Dy0yf_change * change_indicator +
## Dy0yf_chXlake * change_indicator * lake_indicator + Dy0yf_chXage *
## change_indicator * age_indicator) * exp(-(lambda + lambda_lake *
## lake_indicator + lambda_age * age_indicator + lambda_change *
## change_indicator + lambda_CHxLa * change_indicator * lake_indicator +
## lambda_chXage * change_indicator * age_indicator + lambda_lakeXage *
## lake_indicator * age_indicator + lambda_lakeXageXch * change_indicator *
## lake_indicator * age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -80.27467 5.81328 -13.809 < 2e-16 ***
## Dy0yf_lake -11.34395 6.70277 -1.692 0.0935 .
## Dy0yf_age 16.63664 6.65021 2.502 0.0139 *
## Dy0yf_change 146.01258 7.20914 20.254 < 2e-16 ***
## Dy0yf_chXlake 16.10642 8.25727 1.951 0.0537 .
## Dy0yf_chXage -20.28437 8.04164 -2.522 0.0131 *
## yf 463.34057 1.86664 248.222 < 2e-16 ***
## yf_lake -3.29869 2.10915 -1.564 0.1208
## yf_age 1.17095 1.84790 0.634 0.5277
## yf_change -72.18381 2.28667 -31.567 < 2e-16 ***
## lambda 0.72858 0.12560 5.801 6.78e-08 ***
## lambda_lake 0.06425 0.16662 0.386 0.7005
## lambda_age 0.94339 0.55531 1.699 0.0923 .
## lambda_change -0.58221 0.13169 -4.421 2.37e-05 ***
## lambda_CHxLa 0.20740 0.18210 1.139 0.2573
## lambda_chXage -1.08839 0.55585 -1.958 0.0528 .
## lambda_lakeXage -1.21818 0.55883 -2.180 0.0315 *
## lambda_lakeXageXch 1.32200 0.56060 2.358 0.0202 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.683 on 107 degrees of freedom
##
## Number of iterations to convergence: 12
## Achieved convergence tolerance: 1.49e-08# noDy0yf_lakeXageXch # noDy0yf_lakeXage #noyf_lakeXage
# no yf_chXage # yf_chXlake # noyf_lakeXageXch
# no yf_lake yf_age
afm_v6 <- nlsLM(y ~ (yf + yf_change * change_indicator) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator+ Dy0yf_change * change_indicator + Dy0yf_chXlake * change_indicator* lake_indicator + Dy0yf_chXage * change_indicator* age_indicator) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator+ lambda_change * change_indicator + lambda_CHxLa * change_indicator* lake_indicator + lambda_chXage * change_indicator* age_indicator + lambda_lakeXage * lake_indicator* age_indicator + lambda_lakeXageXch * change_indicator * lake_indicator* age_indicator) * t),
data = df,
start = list(Dy0yf = -75, Dy0yf_lake = -20, Dy0yf_age = -30, Dy0yf_change = 120, Dy0yf_chXlake = 60, Dy0yf_chXage=0, yf =460,yf_change = -50,
lambda = 0.8, lambda_lake = 0, lambda_age = -0.7, lambda_change = -0.4, lambda_CHxLa = -.2, lambda_chXage=1, lambda_lakeXage=0, lambda_lakeXageXch=0),
control=nls.lm.control(maxiter = 100))
summary(afm_v6)##
## Formula: y ~ (yf + yf_change * change_indicator) + (Dy0yf + Dy0yf_lake *
## lake_indicator + Dy0yf_age * age_indicator + Dy0yf_change *
## change_indicator + Dy0yf_chXlake * change_indicator * lake_indicator +
## Dy0yf_chXage * change_indicator * age_indicator) * exp(-(lambda +
## lambda_lake * lake_indicator + lambda_age * age_indicator +
## lambda_change * change_indicator + lambda_CHxLa * change_indicator *
## lake_indicator + lambda_chXage * change_indicator * age_indicator +
## lambda_lakeXage * lake_indicator * age_indicator + lambda_lakeXageXch *
## change_indicator * lake_indicator * age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -79.39547 5.74015 -13.832 < 2e-16 ***
## Dy0yf_lake -14.26760 6.51851 -2.189 0.03075 *
## Dy0yf_age 17.79705 6.51154 2.733 0.00732 **
## Dy0yf_change 147.14222 7.17344 20.512 < 2e-16 ***
## Dy0yf_chXlake 16.16290 8.27980 1.952 0.05349 .
## Dy0yf_chXage -19.78922 8.04010 -2.461 0.01541 *
## yf 462.39398 1.07764 429.082 < 2e-16 ***
## yf_change -73.87797 2.06762 -35.731 < 2e-16 ***
## lambda 0.75186 0.12450 6.039 2.19e-08 ***
## lambda_lake -0.01484 0.15213 -0.098 0.92246
## lambda_age 1.13677 0.70988 1.601 0.11220
## lambda_change -0.61876 0.12555 -4.929 2.98e-06 ***
## lambda_CHxLa 0.31005 0.16668 1.860 0.06556 .
## lambda_chXage -1.26882 0.71009 -1.787 0.07674 .
## lambda_lakeXage -1.37863 0.71589 -1.926 0.05674 .
## lambda_lakeXageXch 1.45336 0.71963 2.020 0.04588 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.721 on 109 degrees of freedom
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08full_additive_model <- nlsLM(y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator ) +
(Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator ) *
exp(-(lambda + lambda_lake * lake_indicator + lambda_age * age_indicator) * t),
data = df,
start = list(Dy0yf = -95, Dy0yf_lake = -34, Dy0yf_age = -34,
yf =460, yf_lake = 1, yf_age = 1 ,
lambda = 0.001, lambda_lake = 0, lambda_age = 0.001),
control=nls.lm.control(maxiter = 100))
summary(full_additive_model)##
## Formula: y ~ (yf + yf_lake * lake_indicator + yf_age * age_indicator) +
## (Dy0yf + Dy0yf_lake * lake_indicator + Dy0yf_age * age_indicator) *
## exp(-(lambda + lambda_lake * lake_indicator + lambda_age *
## age_indicator) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -0.1760 1.0565 -0.167 0.8680
## Dy0yf_lake 22.5229 13.8528 1.626 0.1067
## Dy0yf_age -27.5546 16.3365 -1.687 0.0944 .
## yf 439.3958 8.5595 51.334 < 2e-16 ***
## yf_lake -34.5946 7.4860 -4.621 9.96e-06 ***
## yf_age 19.4932 10.2489 1.902 0.0597 .
## lambda -0.3609 0.4124 -0.875 0.3833
## lambda_lake 0.3748 0.4092 0.916 0.3616
## lambda_age 2.3493 5.3020 0.443 0.6585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.04 on 116 degrees of freedom
##
## Number of iterations to convergence: 75
## Achieved convergence tolerance: 1.49e-08null_model <- nlsLM(y ~ (yf ) +
(Dy0yf ) *
exp(-(lambda ) * t),
data = df,
start = list(Dy0yf = -75,
yf =460,
lambda = 0.8))
summary(null_model)##
## Formula: y ~ (yf) + (Dy0yf) * exp(-(lambda) * t)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Dy0yf -13.795 11.263 -1.225 0.223
## yf 435.236 2.979 146.087 <2e-16 ***
## lambda 2.596 10.926 0.238 0.813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.73 on 122 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.49e-08MODEL COMPARISON
# Function to compute AIC, BIC, and Degrees of Freedom for nlsLM models
compute_AIC_BIC <- function(model) {
logLik_val <- logLik(model) # Log-likelihood
k <- length(coef(model)) # Number of estimated parameters (degrees of freedom)
AIC_val <- AIC(model)
BIC_val <- BIC(model)
return(data.frame(Model = deparse(substitute(model)), AIC = AIC_val, BIC = BIC_val, DF = k))
}# Function to compute AIC, BIC, and Degrees of Freedom for nlsLM models
compute_AIC_BIC <- function(model) {
logLik_val <- logLik(model) # Log-likelihood
k <- length(coef(model)) # Number of estimated parameters (degrees of freedom)
AIC_val <- AIC(model)
BIC_val <- BIC(model)
return(data.frame(Model = deparse(substitute(model)), AIC = AIC_val, BIC = BIC_val, DF = k))
}
# Compare all models
aic_bic_results <- rbind(
compute_AIC_BIC(full_model),
compute_AIC_BIC(afm_v6),
compute_AIC_BIC(afm_v3),
compute_AIC_BIC(afm_v4),
compute_AIC_BIC(afm_v5),
compute_AIC_BIC(full_additive_model),
compute_AIC_BIC(null_model)
)
# Sort by AIC (best model has the lowest AIC)
aic_bic_results <- aic_bic_results[order(aic_bic_results$AIC), ]
# Display results
print(aic_bic_results)## Model AIC BIC DF
## 2 afm_v6 847.9268 896.0081 16
## 5 afm_v5 848.2039 901.9419 18
## 4 afm_v4 848.9810 908.3756 20
## 3 afm_v3 849.0016 914.0528 22
## 1 full_model 850.2318 920.9396 24
## 6 full_additive_model 1198.8407 1227.1239 9
## 7 null_model 1216.0468 1227.3601 3### BEST MODEL
bmod <- afm_v6
# Compute residual sum of squares (RSS)
rss <- sum(residuals(bmod)^2)
# Compute total sum of squares (TSS)
tss <- sum((df$y - mean(df$y))^2)
# Compute pseudo R²
pseudo_R2 <- 1 - (rss/tss)
pseudo_R2 # 0.96## [1] 0.9577936# Compute RMSE
rmse <- sqrt(mean(residuals(bmod)^2))
rmse # 6.27## [1] 6.276224DIAGNOSTICS PLOT FOR MODEL
par(mfrow=c(2,2))
# Residuals vs. Fitted
plot(fitted(bmod), residuals(bmod),
main="Residuals vs. Fitted", xlab="Fitted values", ylab="Residuals")
# QQ-plot for normality
qqnorm(residuals(bmod))
qqline(residuals(bmod))
# Histogram of residuals
hist(residuals(bmod), main="Histogram of Residuals", breaks=20)
# Residuals vs. time (if `t` is a predictor)
plot(df$t, residuals(bmod), main="Residuals over Time",
xlab="Time", ylab="Residuals")
EXTRACT MODEL ESTIMATES AND CALCULATE STANDARD ERRORS
# Extract estimates and standard errors
estimates <- summary(bmod)$coefficients[, "Estimate"]
std_errors <- summary(bmod)$coefficients[, "Std. Error"]
# Create a data frame for plotting
params <- data.frame(
Parameter = names(estimates),
Estimate = estimates,
Std_Error = std_errors
)
# Variance-covariance matrix
vcov_matrix <- vcov(bmod)
# Calculate the combined standard error for the sum of two parameters (e.g., lambda_lake & lambda_age)
Nic_late_se <- sqrt(vcov_matrix["lambda", "lambda"] +
vcov_matrix["lambda_age", "lambda_age"] +
2 * vcov_matrix["lambda", "lambda_age"])
Xiloa_early_se <- sqrt(vcov_matrix["lambda", "lambda"] +
vcov_matrix["lambda_lake", "lambda_lake"] +
2 * vcov_matrix["lambda", "lambda_lake"])
xiloa_late_se <- sqrt(
vcov_matrix["lambda", "lambda"] +
vcov_matrix["lambda_lake", "lambda_lake"] +
vcov_matrix["lambda_age", "lambda_age"] +
vcov_matrix["lambda_lakeXage", "lambda_lakeXage"] +
2 * vcov_matrix["lambda_lake", "lambda_age"] +
2 * vcov_matrix["lambda", "lambda_age"] +
2 * vcov_matrix["lambda_lake", "lambda"] +
2 * vcov_matrix["lambda", "lambda_lakeXage"] +
2 * vcov_matrix["lambda_lake", "lambda_lakeXage"] +
2 * vcov_matrix["lambda_age", "lambda_lakeXage"]
)PLOT MODEL ESTIMATES
## READ PARAMETERS GENERATED FROM MODEL
plasticity <- read.table("Plasticity_Estimates_Full_Model.csv", header=T, sep =",", dec=".",stringsAsFactors = T) [,-1]
co <- palette.colors()
# Convert categorical variables to factors with defined levels
plasticity$lake <- factor(plasticity$lake, levels = c("Nic", "Xil"))
plasticity$age <- factor(plasticity$age, levels = c("E", "L"))
plasticity$change <- factor(plasticity$change, levels = c("WR", "RW"))
plasticity$Pl_cap <- abs(plasticity$Pl_cap)
#### Plot
# Define color mapping for lakes
#lake_colors <- c("Nic" = "red", "Xil" = "blue")
# change_linetypes <- c("RW" = "solid", "WR" = "dashed")
age_alpha <- c("E" = 1, "L" = 0.5)
lake_colors <- c("Nic" = co[7], "Xil" = co[6])
change_linetypes <- c("WR" = "solid", "RW" = "solid")
# Create a new variable to simplify x-axis labels
#plasticity$x_labels <- ifelse(plasticity$age == "E", "Early", "Late")
plasticity$x_labels <- factor(ifelse(plasticity$change == "WR", "White-to-red", "Red-to-white"),
levels = c( "White-to-red", "Red-to-white"))
# Ensure data is sorted correctly for line plotting
#plasticity <- plasticity[order(plasticity$lake, plasticity$x_labels, plasticity$change), ]
# Create subsets for WR and RW
plasticity_WR <- subset(plasticity, change == "WR")
plasticity_RW <- subset(plasticity, change == "RW")
plasticity_E <- subset(plasticity, age == "E")
plasticity_L <- subset(plasticity, age == "L")
plasticity_L[1,2] <- plasticity_L[1,2] + 0.615
# Define position dodge for x-offset
pos_dodge <- position_dodge(width = 0.2)
# PRATE E (Panel A)
ggplot(plasticity_E, aes(x = x_labels, y = Pl_rate, color = lake, linetype = change)) +
geom_line(aes(group = lake), linewidth = 1.2, position = pos_dodge) + # Connect Nic & Xil
geom_point(aes(alpha = age), size = 3, position = pos_dodge) + # Apply jitter-like effect
geom_errorbar(aes(ymin = Pl_rate - Pl_rate_se, ymax = Pl_rate + Pl_rate_se), width = 0.2, position = pos_dodge) +
scale_color_manual(values = lake_colors) +
scale_linetype_manual(values = change_linetypes) +
scale_alpha_manual(values = age_alpha) +
labs(x = "", y = expression("Plasticity Rate (" * lambda * ")"), title = "") +
coord_cartesian(ylim = c(0, 2.5)) +
theme_classic()