-
Notifications
You must be signed in to change notification settings - Fork 0
/
nlr.R
119 lines (94 loc) · 5.31 KB
/
nlr.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
# =============================================================================
# Nonlinear Regression
# Supporting Information for Early-Capistrán et al. (2020), PeerJ
# [email protected] - April 2020
# =============================================================================
# =============================================================================
# Install and load libraries
# =============================================================================
# Check if required libraries are installed and install if necessary ..........
packages <- c("easynls", "nlstools", "dplyr", "ggplot2", "ggthemes", "car",
"DescTools")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
install.packages(setdiff(packages, rownames(installed.packages())))
}
# Load libraries ..............................................................
library(easynls)
library(nlstools)
library(car)
library(DescTools)
library(ggplot2)
library(ggthemes)
library(dplyr)
# .............................................................................
# NOTE:
# If you did not open this script directly from ".../quantifying_lek_data_code",
# please define the working directory to ".../quantifying_lek_data_code":
#
# setwd(".../quantifying_lek_data_code/")
# .............................................................................
source("utils/resi_utils.R")
# =============================================================================
# Load and prepare data
# =============================================================================
# Load dataset from csv file ..................................................
st_data = read.csv("data/standardised_dataset.csv", header = TRUE)
# Standardised CPUE data can be found in ~/datasets/standardised_dataset.csv
# Select columns with X (serialYear) and Y (cpue) variables ...................
st_data <- select(st_data, yearSerial, stCpue)
# Remove rows with missing values .............................................
st_data = st_data[complete.cases(st_data),]
# =============================================================================
# Fitting the exponential model "y~a*exp(b*x)" to standardised CPUE values
# Starting values from LabFit 7.2.49 (http://labfit.net/)
# =============================================================================
# Define the model with nlstools ..............................................
exp.model <- as.formula(stCpue~a*exp(b*yearSerial))
# Define starting values and run model ........................................
cpue.exp <- nls(exp.model, start = list(a=0.189111E+02, b=-.264181E+00),
data = st_data, model=TRUE)
# Plot nonlinear regression ...................................................
plotfit(cpue.exp, smooth = TRUE, main="Exponential decay model",
xlab="Year (serialised)", ylab="Mean standardised CPUE (turtles/night)")
# View results ................................................................
overview(cpue.exp)
# easynls provides an R^2 value ...............................................
nlsfit(st_data, model=6, start=c(0.189111E+02 ,-.264181E+00))
# =============================================================================
# Residual analysis
# =============================================================================
# Residual analysis function ..................................................
# Returns normality plot, plot of residual vs. fitted values, autocorrelation
# plot, residual mean, Shapiro-Wilk normality test, Levene Test for
# homogeneity of variance, and Run's test for randomness
resi.analysis(cpue.exp)
# Test residual mean = 0 ......................................................
resi.t.test.nls(cpue.exp)
# Evaluate residual autocorrelation with autocorrelation plot and correlation
# test (residuals vs. lagged residuals) .......................................
# Create a data frame with residuals and lagged residuals
resi <- residuals(cpue.exp)
resDf <- as.data.frame(resi) %>% select(resi)
resDf$resi_lag <- c(resDf$resi[-1], NA) # Append NA to vector of lagged
# residuals to match vector lengths
# Generate autocorrelation plot of residuals with linear trend line
ggplot(aes(x=resi, y=resi_lag), data=resDf) + geom_point(size=2, na.rm=TRUE) +
geom_smooth(method = "lm", se=FALSE, na.rm=TRUE) +
labs(y = "Lagged residuals", x = "Residuals")
# Run a linear regression of residuals vs. lagged residuals
lag_lm <- lm(resDf$resi ~ resDf$resi_lag)
summary(lag_lm)
# Run Pearson correlation test for residuals vs. lagged residuals
cor.test(resDf$resi, resDf$resi_lag, method = "pearson")
# =============================================================================
# Custom NLR plot with ggplot
# =============================================================================
# Plot with ggplot ............................................................
custom_plot <- ggplot(st_data, aes(x = yearSerial, y = stCpue), ylim=20) +
geom_point(size=2.5, colour = "steelblue", na.rm=TRUE) +
geom_smooth(method="nls", formula = y ~ a*exp(b*x), na.rm=TRUE,
method.args = list(start = c(a=0.189111E+02, b=-.264181E+00)),
se = F, colour = "grey35") +
labs(y = "Mean standardised CPUE (turtles/night)", x = "Year") +
theme_hc()
custom_plot