AR-WethPLot2_22_beta
Ralf Schwamborn
12 de abril de 2017
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
# AR-WethPlot2 v.2.2 #
### Automatic Routines for Powell-Wetherall plot methods
# (Copyright: R. Schwamborn, 2017)
# AR-WethPlot2 is a script (i.e., a sequence of commands) for the "R" language and environment. It is an automatic routine ("A.R.") for the analysis of length-frequency-distributions (LFDs) with the Wetherall Plot ("Wethplot") methods in "R". The data for regression are selected by an automatic "gamma" post-hoc selection routine. After calculation of mean lengths, it selects a specific data range from L = LFD mode+10% to L = 93%Lmax (="gamma" selection), fits a regression model on these points and calculates the Wetherall estimates (Linf and Z/K). AR-Wethplot allows the application of the
# original Powell-Wetherall plot ("P-W plot", Wetherall, 1986) and of the
# "modified" version of the method ("mP-W plot", Pauly, 1986),
# with and without gamma selection.
# Available (free) at: http://www.ufpe.br/ecocost/AR_WETHPLOT
# Functions: mP_W.fun, P_W.fun, mP_W.ffun, P_W.ffun
# explanations: see below
#### Test Data
### 1. Generate test data
# Lengths and Catch data
L <- seq(4.5,17.5)
Catch <- c(2, 4, 16, 28,83,116,102,75,48,28,12,5,2, 1)
plot(L,Catch)
### FUNCTIONS ###
#### I. FUNCTIONs of Wetherall methods without using selection of points
## (uses ALL data points for regression)
######## FUNCTION "mP_W.fun"
#### I.a, FUNCTION: mP-W plot as Function
## "modified" version of the method ("mP-W plot", Pauly, 1986),
# without gamma selction (uses ALL Data of the plot)
mP_W.fun <- function(Len, Cat) {
N_c <- length(Len) # N of length classes
int_Len <- (max(Len)-min(Len))/ (N_c-1) # interval (class width)
L_c <- Len - (int_Len/2) # cutoff lengths (L - 0.5)
# L_c <- Len - 0.5
# plot(Len, Cat)
### Calculate Mean lengths
#Catch*L
Catch_L <- Cat * Len
#Sum of "Catch*L" above each Lc
#write loop
Sum_of_Catch_times_L <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch_times_L[i] <- sum(Catch_L[i:N_c]);
}
#Sum of Catch above Lc
#write loop
Sum_of_Catch <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch[i] <- sum(Cat[i:N_c]);
}
# Calculate mean length for each Lc
mean_length <- Sum_of_Catch_times_L / Sum_of_Catch
#Pauly's modified P-W_plot method
# mP-W plot
mean_length # mean length
L_c # cutoff length
plot((mean_length - L_c) ~ L_c)
mP_W_mod1 <- lm((mean_length-L_c) ~ L_c)
a.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[1:1]
b.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[2,1]
# Linf = a / -b and Z/K = (1 + b)/(- b)
L_inf_mP_W <- a.mP_W_mod1 / (-1 * b.mP_W_mod1)
Z_K_mP_W = (1 + b.mP_W_mod1)/(-1 * b.mP_W_mod1)
; c("LinF:", L_inf_mP_W, "Z/K:", Z_K_mP_W, "mP-W method") # output
}
mP_W.fun(L, Catch)
## [1] "LinF:" "16.495993981392" "Z/K:"
## [4] "1.51721432168017" "mP-W method"######## FUNCTION "P_W.fun"
#### I.b, FUNCTION P-W plot as Function
## original, "unmodifed" version of the method ("P-W plot", Wetherall, 1986),
# without gamma selction (uses ALL Data of the plot)
#### FUNCTION original P-W method as Function
P_W.fun <- function(Len, Cat) {
N_c <- length(Len) # N of length classes
int_Len <- (max(Len)-min(Len))/ (N_c-1) # interval (class width)
L_c <- Len - (int_Len/2) # cutoff lengths (e.g. L - 0.5)
# L_c <- Len - 0.5
# plot(Len, Cat)
### Calculate Mean lengths
#Catch*L
Catch_L <- Cat * Len
#Sum of "Catch*L" above each Lc
#write loop
Sum_of_Catch_times_L <- rep (0, N_c)
for(i in 1:N_c){
Sum_of_Catch_times_L[i] <- sum(Catch_L[i:N_c]);
}
#Sum of Catch above Lc
#write loop
Sum_of_Catch <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch[i] <- sum(Cat[i:N_c]);
}
# Calculate mean length for each Lc
mean_length <- Sum_of_Catch_times_L / Sum_of_Catch
plot(mean_length ~ L_c)
# linear model
P_W_mod1 <- lm(mean_length ~ L_c)
a.P_W_mod1 <- summary(P_W_mod1)$coefficients[1:1]
b.P_W_mod1 <- summary(P_W_mod1)$coefficients[2,1]
#calculate Linf and Z/K
# Linf = a / (1 - b), Z/K = b / (1 - b)
L_inf_P_W <- a.P_W_mod1 / (1 - b.P_W_mod1)
Z_K_P_W = (b.P_W_mod1)/(1 - b.P_W_mod1)
; c("LinF:", L_inf_P_W, "Z/K:", Z_K_P_W, "P-W method") # output
} #### EOF # end of Function "P_W.fun"
P_W.fun(L, Catch)
## [1] "LinF:" "16.495993981392" "Z/K:"
## [4] "1.51721432168017" "P-W method"#################################
#### II. FUNCTIONs of Wetherall methods using "gamma" selection of points
### selection of points for regression
## from mode+10% to 93%Lmax
#
######## FUNCTION "mP_W.ffun"
#### II.a, FUNCTION mP-W plot as Function "mP_W.ffun"
## original, "unmodifed" version of the method ("P-W plot", Wetherall, 1986),
# without gamma selction (uses ALL Data of the plot)
mP_W.ffun <- function(Len, Cat) {
N_c <- length(Len) # N of length classes
int_Len <- (max(Len)-min(Len))/ (N_c-1) # interval (class width)
L_c <- Len - (int_Len/2) # cutoff lengths (L - 0.5)
# L_c <- Len - 0.5
# plot(Len, Cat)
### Calculate Mean lengths
#Catch*L
Catch_L <- Cat * Len
#Sum of "Catch*L" above each Lc
#write loop
Sum_of_Catch_times_L <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch_times_L[i] <- sum(Catch_L[i:N_c]);
}
#Sum of Catch above Lc
#write loop
Sum_of_Catch <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch[i] <- sum(Cat[i:N_c]);
}
# Calculate mean length for each Lc
mean_length <- Sum_of_Catch_times_L / Sum_of_Catch
#Pauly's modified P-W_plot method
# mP-W plot
# GAMMA selection: select data from Mode+20% to 90%Lmax
mean_length # mean length
L_c # cutoff length
df_1 <- data.frame(mean_length, L_c, Len, Cat)
d_mode <- df_1$Len[df_1$Cat == max(df_1$Cat)] # mode of original LFD data
df_gamma <- subset(df_1, Len > (1.1 * d_mode))
df_gamma <- subset(df_gamma, Len < (0.93 * max(Len)))
mean_length_G <- df_gamma$mean_length # Gamma selected mean lengths
L_c_G <- df_gamma$L_c # Gamma selected cutoff lengths
y_p <- (mean_length - L_c)
y_p_G <- (mean_length_G - L_c_G)
plot((mean_length - L_c) ~ L_c)
# plot((mean_length - L_c) ~ L_c, xlim = c(min(L_c), max(L_c)), ylim = c(min(y_p), max(y_p)))
points((mean_length_G - L_c_G) ~ L_c_G, pch = 16, col = "blue")
mP_W_mod1 <- lm((mean_length_G - L_c_G) ~ L_c_G)
abline(mP_W_mod1, col = "blue")
a.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[1:1]
b.mP_W_mod1 <- summary(mP_W_mod1)$coefficients[2,1]
# Linf = a / -b and Z/K = (1 + b)/(- b)
L_inf_mP_W <- a.mP_W_mod1 / (-1 * b.mP_W_mod1)
Z_K_mP_W = (1 + b.mP_W_mod1)/(-1 * b.mP_W_mod1)
; c("LinF:", L_inf_mP_W, "Z/K:", Z_K_mP_W, "mP-W method") # output
} #### EOF # end of Function "mP_W.ffun"
mP_W.ffun(L, Catch)
## [1] "LinF:" "21.1794203958001" "Z/K:"
## [4] "5.59171283957967" "mP-W method"######## FUNCTION "P_W.ffun"
#### II.b. Original "unmodified" P-W method (Wetherall, 1986)
## FUNCTION P-W method as Function with blue colour plot
# and GAMMA selection mode+10% to 93%Lmax
P_W.ffun <- function(Len, Cat) {
N_c <- length(Len) # N of length classes
int_Len <- (max(Len)-min(Len))/ (N_c-1) # interval (class width)
L_c <- Len - (int_Len/2) # cutoff lengths (e.g. L - 0.5)
# L_c <- Len - 0.5
# plot(Len, Cat)
### Calculate Mean lengths
#Catch*L
Catch_L <- Cat * Len
#Sum of "Catch*L" above each Lc
#write loop
Sum_of_Catch_times_L <- rep (0, N_c)
for(i in 1:N_c){
Sum_of_Catch_times_L[i] <- sum(Catch_L[i:N_c]);
}
#Sum of Catch above Lc
#write loop
Sum_of_Catch <- rep(0, N_c)
for(i in 1:N_c){
Sum_of_Catch[i] <- sum(Cat[i:N_c]);
}
# Calculate mean length for each Lc
mean_length <- Sum_of_Catch_times_L / Sum_of_Catch
#plot(mean_length ~ L_c)
# linear model
# GAMMA selection: select data from Mode+20% to 90%Lmax
mean_length # mean length
L_c # cutoff length
df_1 <- data.frame(mean_length, L_c, Len, Cat)
d_mode <- df_1$Len[df_1$Cat == max(df_1$Cat)] # mode of original LFD data
df_gamma <- subset(df_1, Len > (1.1 * d_mode))
df_gamma <- subset(df_gamma, Len < (0.93 * max(Len)))
mean_length_G <- df_gamma$mean_length # Gamma selected mean lengths
L_c_G <- df_gamma$L_c # Gamma selected cutoff lengths
y_p <- (mean_length)
y_p_G <- (mean_length_G)
plot((mean_length) ~ L_c, xlim = c(min(L_c), max(L_c)),
ylim = c(min(y_p), max(y_p)))
points((mean_length_G) ~ L_c_G, pch = 16, col = "blue")
P_W_mod1 <- lm((mean_length_G) ~ L_c_G)
abline(P_W_mod1, col = "blue")
a.P_W_mod1 <- summary(P_W_mod1)$coefficients[1:1]
b.P_W_mod1 <- summary(P_W_mod1)$coefficients[2,1]
#calculate Linf and Z/K
# Linf = a / (1 - b), Z/K = b / (1 - b)
L_inf_P_W <- a.P_W_mod1 / (1 - b.P_W_mod1)
Z_K_P_W = (b.P_W_mod1)/(1 - b.P_W_mod1)
; c("LinF:", L_inf_P_W, "Z/K:", Z_K_P_W, "P-W method") # output
} #### EOF # end of Function "P_W.ffun"
P_W.ffun(L, Catch)
## [1] "LinF:" "21.1794203958001" "Z/K:"
## [4] "5.59171283957969" "P-W method" ##########################################################
##########################################################
##########################################################
##########################################################
## test and compare with TropfishR (to run, remove "#" from lines below)
# library(TropFishR)
#select month 1 to12
# data(synLFQ5)
# powell_wetherall(synLFQ5, catch_columns = 1:12)
# Cat_data_1y <- transform(Cat_data_1y, SUM=rowSums(Cat_data_1y))
#View(Cat_data_1y)
#plot(synLFQ5$midLengths, Cat_data_1y$SUM) # Z/K = 2.135 for a straight line
# mP_W.fun(L, Catch)
# mP_W.fun(synLFQ5$midLengths, Cat_data_1y$SUM) # without GAMMA selection
# mP_W.ffun(synLFQ5$midLengths, Cat_data_1y$SUM) # with GAMMA selection, Z/K: 2.24
#P_W.fun(L, Catch)
#P_W.fun(synLFQ5$midLengths, Cat_data_1y$SUM) # without GAMMA selection
#P_W.ffun(synLFQ5$midLengths, Cat_data_1y$SUM) # with GAMMA selection, Z/K: 2.24
##########################################################