Mangochi_Fishery
Paul Tuda
2024-05-17
Setup
set default chunk options
Load packages
library(tidyverse)
library(readr)
library(dplyr)
library(stats)
library(ggplot2)
library(reshape2)
library(vegan)
library(ggforce)
library(forecast)
library(tseries)Load the data only for Mangochi
## # A tibble: 6 × 23
## Year Gear Statistics Chambo Kasawala Kambuzi Mbaba Kampango Mcheni Sanjika
## <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2002 Gill… Catch 1.43e3 4.24e+1 2.31e1 6.46e+2 8.4 e+1 2.13e1 0
## 2 2002 Gill… Effort 1.06e6 1.06e+6 1.06e6 1.06e+6 1.06e+6 1.06e6 1058426
## 3 2002 Gill… CPUE 1.4 e0 0 0 6 e-1 1 e-1 0 0
## 4 2002 Kamb… Catch 1.30e2 4 e-1 3.55e2 1.69e+2 2 e-1 1.9 e1 0
## 5 2002 Kamb… Effort 4.80e3 4.80e+3 4.80e3 4.80e+3 4.80e+3 4.80e3 4804
## 6 2002 Kamb… CPUE 2.71e1 1 e-1 7.38e1 3.52e+1 0 3.9 e0 0
## # ℹ 13 more variables: Ntchila <dbl>, Mpasa <dbl>, Mlamba <dbl>, Utaka <dbl>,
## # Usipa <dbl>, Chisawasawa <dbl>, Nkholokolo <dbl>, Ndunduma <dbl>,
## # `Other Tilapia` <dbl>, Others <dbl>, ...21 <lgl>, ...22 <chr>, ...23 <dbl>## spc_tbl_ [612 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Year : num [1:612] 2002 2002 2002 2002 2002 ...
## $ Gear : chr [1:612] "Gillnet" "Gillnet" "Gillnet" "Kambuzi seine" ...
## $ Statistics : chr [1:612] "Catch" "Effort" "CPUE" "Catch" ...
## $ Chambo : num [1:612] 1.43e+03 1.06e+06 1.40 1.30e+02 4.80e+03 ...
## $ Kasawala : num [1:612] 4.24e+01 1.06e+06 0.00 4.00e-01 4.80e+03 ...
## $ Kambuzi : num [1:612] 2.31e+01 1.06e+06 0.00 3.55e+02 4.80e+03 ...
## $ Mbaba : num [1:612] 6.46e+02 1.06e+06 6.00e-01 1.69e+02 4.80e+03 ...
## $ Kampango : num [1:612] 8.40e+01 1.06e+06 1.00e-01 2.00e-01 4.80e+03 ...
## $ Mcheni : num [1:612] 2.13e+01 1.06e+06 0.00 1.90e+01 4.80e+03 ...
## $ Sanjika : num [1:612] 0 1058426 0 0 4804 ...
## $ Ntchila : num [1:612] 0 1058426 0 0 4804 ...
## $ Mpasa : num [1:612] 8.30 1.06e+06 0.00 0.00 4.80e+03 ...
## $ Mlamba : num [1:612] 9.20e+01 1.06e+06 1.00e-01 1.00 4.80e+03 ...
## $ Utaka : num [1:612] 6.55e+02 1.06e+06 6.00e-01 0.00 4.80e+03 ...
## $ Usipa : num [1:612] 5.45e+02 1.06e+06 5.00e-01 1.00e-01 4.80e+03 ...
## $ Chisawasawa : num [1:612] 1.30e+02 1.06e+06 1.00e-01 0.00 4.80e+03 ...
## $ Nkholokolo : num [1:612] 8.80 1.06e+06 0.00 0.00 4.80e+03 ...
## $ Ndunduma : num [1:612] 1.98e+02 1.06e+06 2.00e-01 0.00 4.80e+03 ...
## $ Other Tilapia: num [1:612] 1.20e+02 1.06e+06 1.00e-01 0.00 4.80e+03 ...
## $ Others : num [1:612] 1.27e+01 1.06e+06 0.00 0.00 4.80e+03 ...
## $ ...21 : logi [1:612] NA NA NA NA NA NA ...
## $ ...22 : chr [1:612] NA NA NA NA ...
## $ ...23 : num [1:612] NA NA NA NA NA NA NA NA NA NA ...
## - attr(*, "spec")=
## .. cols(
## .. Year = col_double(),
## .. Gear = col_character(),
## .. Statistics = col_character(),
## .. Chambo = col_double(),
## .. Kasawala = col_double(),
## .. Kambuzi = col_double(),
## .. Mbaba = col_double(),
## .. Kampango = col_double(),
## .. Mcheni = col_double(),
## .. Sanjika = col_double(),
## .. Ntchila = col_double(),
## .. Mpasa = col_double(),
## .. Mlamba = col_double(),
## .. Utaka = col_double(),
## .. Usipa = col_double(),
## .. Chisawasawa = col_double(),
## .. Nkholokolo = col_double(),
## .. Ndunduma = col_double(),
## .. `Other Tilapia` = col_double(),
## .. Others = col_double(),
## .. ...21 = col_logical(),
## .. ...22 = col_character(),
## .. ...23 = col_double()
## .. )
## - attr(*, "problems")=<externalptr># Load the data frame previously defined
# Assuming `df` is already loaded
# Remove unwanted columns
df <- select(df, -c(...21, ...22, ...23)) # Adjust column names as needed
# Check for missing values in the dataframe
missing_values <- colSums(is.na(df))
print(missing_values)## Year Gear Statistics Chambo Kasawala
## 0 0 0 2 2
## Kambuzi Mbaba Kampango Mcheni Sanjika
## 2 2 2 2 2
## Ntchila Mpasa Mlamba Utaka Usipa
## 2 2 2 2 2
## Chisawasawa Nkholokolo Ndunduma Other Tilapia Others
## 2 2 2 2 2Verify and Clean Gear Types
Before analyzing the data, it’s essential to ensure that the
Gear column is free from typographical errors and
standardized for consistency. This section outlines the steps taken to
clean and standardize the Gear column.
# Display unique gear types before cleaning to identify inconsistencies
original_gear_types <- unique(df$Gear)
print(original_gear_types)## [1] "Gillnet" "Kambuzi seine" "Chambo seine" "Chilimira"
## [5] "Mosquito net" "Kandwindwi" "Nkacha" "Longline"
## [9] "Handline" "FishTrap" "Cast Net" "Cast net"
## [13] "Fish Trap" "Usipa seine" "chambo seine" "Mosquito Net"# Convert gear types to lower case to avoid case-sensitive mismatches
df$Gear <- tolower(df$Gear)
# Replace spaces with underscores to ensure consistent naming conventions
df$Gear <- str_replace_all(df$Gear, " ", "_")
# Recode similar terms to a standard form using the dplyr recode function
df$Gear <- recode(df$Gear,
"cast_net" = "cast_net",
"fishtrap" = "fish_trap",
"mosquito net" = "mosquito_net",
"chambo seine" = "chambo_seine")
# Display unique gear types after cleaning to verify changes
cleaned_gear_types <- unique(df$Gear)
print(cleaned_gear_types)## [1] "gillnet" "kambuzi_seine" "chambo_seine" "chilimira"
## [5] "mosquito_net" "kandwindwi" "nkacha" "longline"
## [9] "handline" "fish_trap" "cast_net" "usipa_seine"Data Transformation
Convert Data from Wide to Long Format
Data needs to be in long format to facilitate easier aggregation and analysis by species. This transformation will make subsequent analyses more straightforward.
# Convert the dataset from wide to long format, focusing on species columns
dat <- df %>%
pivot_longer(cols = Chambo:Others, names_to = "species", values_to = "total_catch")
# Display the first few rows of the transformed data to verify changes
head(dat)## # A tibble: 6 × 5
## Year Gear Statistics species total_catch
## <dbl> <chr> <chr> <chr> <dbl>
## 1 2002 gillnet Catch Chambo 1433
## 2 2002 gillnet Catch Kasawala 42.4
## 3 2002 gillnet Catch Kambuzi 23.1
## 4 2002 gillnet Catch Mbaba 646.
## 5 2002 gillnet Catch Kampango 84
## 6 2002 gillnet Catch Mcheni 21.3Summarize Data by Year
After transforming the data, summarize the total catches by year for each species identified as a ‘Catch’ in the ‘Statistics’ column (Catch, Effort, CPUE).
# Filter data for entries marked as 'Catch' and group by year
total_catch_by_year <- dat %>%
filter(Statistics == "Catch") %>%
group_by(Year) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE))
# Display the summarized data
print(total_catch_by_year)## # A tibble: 21 × 2
## Year total_catch
## <dbl> <dbl>
## 1 2002 10621.
## 2 2003 11217
## 3 2004 13136
## 4 2005 19977.
## 5 2006 19493.
## 6 2007 30006.
## 7 2008 28121.
## 8 2009 28447.
## 9 2010 36060.
## 10 2011 29234.
## # ℹ 11 more rowsPlot Total Catch by Year
# Create a line plot with points showing the total catch by year
ggplot(total_catch_by_year, aes(x = Year, y = total_catch)) +
geom_line(size = 1.2) + # Connect points with lines
geom_point(size = 3, shape = 21, fill = "white") + # Add points with white fill for visibility
geom_smooth(method="lm",se=FALSE,linetype="dashed",color="tomato",linewidth=0.5) +
scale_color_brewer(palette = "Dark2") + # A distinct color palette
labs(title = "Total Catch by Year",
x = "Year",
y = "Total Catch") +
theme_bw() +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12)
)
sort the discrepency in 2016
Modify the dataset to replace 2016 total catch with the calculated mean from 2014 and 2015
total_catch_by_year_mod <- dat %>%
filter(Statistics == "Catch") %>%
mutate(total_catch = if_else(Year == 2016, moving_average_2013_2015, total_catch)) %>%
group_by(Year) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE))
# Display the updated summarized data
print(total_catch_by_year_mod)## # A tibble: 21 × 2
## Year total_catch
## <dbl> <dbl>
## 1 2002 10621.
## 2 2003 11217
## 3 2004 13136
## 4 2005 19977.
## 5 2006 19493.
## 6 2007 30006.
## 7 2008 28121.
## 8 2009 28447.
## 9 2010 36060.
## 10 2011 29234.
## # ℹ 11 more rows### plot new trend
ggplot(total_catch_by_year_mod, aes(x = Year, y = total_catch)) +
geom_line(size = 1.2) + # Connect points with lines
geom_point(size = 3, shape = 21, fill = "white") + # Add points with white fill for visibility
geom_smooth(method="lm",se=FALSE,linetype="dashed",color="tomato",linewidth=0.5) +
scale_color_brewer(palette = "Dark2") + # A distinct color palette
labs(title = "Total Catch by Year",
x = "Year",
y = "Total Catch") +
theme_bw() +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12)
)
# Filter data for entries marked as 'Catch' and group by year
total_cpue_by_year <- dat %>%
filter(Statistics == "CPUE") %>%
group_by(Year) %>%
summarise(total_CPUE = sum(total_catch, na.rm = TRUE))
# Display the summarized data
print(total_cpue_by_year)## # A tibble: 21 × 2
## Year total_CPUE
## <dbl> <dbl>
## 1 2002 506.
## 2 2003 265.
## 3 2004 416.
## 4 2005 2017
## 5 2006 464.
## 6 2007 595.
## 7 2008 1202.
## 8 2009 706.
## 9 2010 652.
## 10 2011 877.
## # ℹ 11 more rowsPlot Total Catch by Year
# Create a line plot with points showing the total catch by year
ggplot(total_cpue_by_year, aes(x = Year, y = total_CPUE)) +
geom_line(size = 1.2) + # Connect points with lines
geom_point(size = 3, shape = 21, fill = "white") + # Add points with white fill for visibility
geom_smooth(method="lm",se=FALSE,linetype="dashed",color="tomato",linewidth=0.5) +
scale_color_brewer(palette = "Dark2") + # A distinct color palette
labs(title = "Total CPUE by Year",
x = "Year",
y = "CPUE") +
theme_bw() +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12)
)
Correct for cpue
Modify the dataset to replace 2016 cpue with the calculated mean from 2019 and 2021
total_cpue_by_year_mod <- dat %>%
filter(Statistics == "CPUE") %>%
mutate(total_catch = if_else(Year == 2020, mean_cpue_2018_2021, total_catch)) %>%
group_by(Year) %>%
summarise(total_cpue = sum(total_catch, na.rm = TRUE))
# Display the updated summarized data
print(total_cpue_by_year_mod)## # A tibble: 21 × 2
## Year total_cpue
## <dbl> <dbl>
## 1 2002 506.
## 2 2003 265.
## 3 2004 416.
## 4 2005 2017
## 5 2006 464.
## 6 2007 595.
## 7 2008 1202.
## 8 2009 706.
## 9 2010 652.
## 10 2011 877.
## # ℹ 11 more rows# Create a line plot with points showing the total catch by year
ggplot(total_cpue_by_year_mod, aes(x = Year, y = total_cpue)) +
geom_line(size = 1.2) + # Connect points with lines
geom_point(size = 3, shape = 21, fill = "white") + # Add points with white fill for visibility
geom_smooth(method="lm",se=FALSE,linetype="dashed",color="tomato",linewidth=0.5) +
scale_color_brewer(palette = "Dark2") + # A distinct color palette
labs(title = "Total CPUE by Year",
x = "Year",
y = "CPUE") +
theme_bw() +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12)
)
### plot trend of catch by gear by year
# Summarize data by year and gear
total_catch_by_year_gear <- dat %>%
filter(Statistics == "Catch") %>%
group_by(Year, Gear) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE))
# Sort the data frame by Year and Gear to ensure correct line connections
total_catch_by_year_gear <- total_catch_by_year_gear %>%
arrange(Gear, Year)
# Plot the updated dataset
ggplot(total_catch_by_year_gear, aes(x = Year, y = total_catch, color = Gear, group = Gear)) +
geom_line(size = 0.5) + # Connect points with lines
geom_point(size = 1.5, shape = 21, fill = "white") + # Add points with white fill for visibility
scale_color_brewer(palette = "Dark2") + # A distinct color palette
labs(title = "Total Catch by Gear and Year",
x = "Year",
y = "Total Catch",
color = "Gear Type") +
theme_bw() +
theme(
text = element_text(size = 12),
plot.title = element_text(face = "bold", size = 16),
plot.subtitle = element_text(size = 12)
)
### Prepare gear_species_wide by transforming
df
gear_species_data <- dat %>%
filter(Statistics == "Catch") %>%
group_by(Gear, species) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE), .groups = 'drop')
gear_species_wide <- gear_species_data %>%
pivot_wider(names_from = species, values_from = total_catch, values_fill = list(total_catch = 0))
# Check the structure of gear_species_wide to ensure it is created correctly
str(gear_species_wide)## tibble [12 × 18] (S3: tbl_df/tbl/data.frame)
## $ Gear : chr [1:12] "cast_net" "chambo_seine" "chilimira" "fish_trap" ...
## $ Chambo : num [1:12] 0 840.4 3150.1 65.9 20926.1 ...
## $ Chisawasawa : num [1:12] 0 0 541 0 2390 ...
## $ Kambuzi : num [1:12] 0 170 3396 186 4017 ...
## $ Kampango : num [1:12] 0 3.6 154.9 0 1747.1 ...
## $ Kasawala : num [1:12] 0 45.6 42.7 91.4 917.4 ...
## $ Mbaba : num [1:12] 0 43.9 370.8 1.9 24426.2 ...
## $ Mcheni : num [1:12] 6.8 7.4 9895.6 0.5 3171.8 ...
## $ Mlamba : num [1:12] 0.3 15.3 166.5 3 3446.5 ...
## $ Mpasa : num [1:12] 0 0 9.8 0 38.7 0 0 0 0 0 ...
## $ Ndunduma : num [1:12] 0 0 3170 0 6584 ...
## $ Nkholokolo : num [1:12] 0 0 225 0 392 ...
## $ Ntchila : num [1:12] 0 0 0 0 24.5 0 0 0 0 0 ...
## $ Other Tilapia: num [1:12] 0 8.9 55.5 19.7 1539.9 ...
## $ Others : num [1:12] 0 0.1 10 8.1 431.7 ...
## $ Sanjika : num [1:12] 0 0 0 0 50.7 0 0 0 0 0.7 ...
## $ Usipa : num [1:12] 227 389 253461 0 1720 ...
## $ Utaka : num [1:12] 0.8 76.8 20333.6 0.2 21733.5 ...Standardize and perform PCA
Scree plot
The scree plot is used to determine the number of factors to retain in an exploratory factor analysis (FA) or principal components to keep in a principal component analysis (PCA).
var_explained <- pca_result$sdev^2 / sum(pca_result$sdev^2) * 100
df_var_explained <- data.frame(PC = 1:length(var_explained), Variance = var_explained)
ggplot(df_var_explained, aes(x = PC, y = Variance)) +
geom_bar(stat = "identity", fill = "steelblue") +
labs(title = "Scree Plot", x = "Principal Component", y = "Percentage of Variance Explained") +
theme_minimal()
NMDS and Dendrogram
Calculate Bray-Curtis dissimilarity and plot NMDS.
## Run 0 stress 0.1011328
## Run 1 stress 0.1550643
## Run 2 stress 0.1677644
## Run 3 stress 0.1011328
## ... Procrustes: rmse 3.431823e-05 max resid 7.210148e-05
## ... Similar to previous best
## Run 4 stress 0.1784278
## Run 5 stress 0.1011328
## ... New best solution
## ... Procrustes: rmse 8.242118e-06 max resid 1.705362e-05
## ... Similar to previous best
## Run 6 stress 0.08652564
## ... New best solution
## ... Procrustes: rmse 0.09079027 max resid 0.2635656
## Run 7 stress 0.1723964
## Run 8 stress 0.1254047
## Run 9 stress 0.08652562
## ... New best solution
## ... Procrustes: rmse 4.610883e-05 max resid 0.0001078343
## ... Similar to previous best
## Run 10 stress 0.1011328
## Run 11 stress 0.1884851
## Run 12 stress 0.08652563
## ... Procrustes: rmse 7.289252e-05 max resid 0.0001781514
## ... Similar to previous best
## Run 13 stress 0.08652562
## ... New best solution
## ... Procrustes: rmse 5.390821e-05 max resid 0.0001331352
## ... Similar to previous best
## Run 14 stress 0.08652563
## ... Procrustes: rmse 3.631957e-05 max resid 8.49155e-05
## ... Similar to previous best
## Run 15 stress 0.1550643
## Run 16 stress 0.08652562
## ... New best solution
## ... Procrustes: rmse 3.079686e-05 max resid 7.624284e-05
## ... Similar to previous best
## Run 17 stress 0.1247858
## Run 18 stress 0.08652562
## ... Procrustes: rmse 7.707074e-06 max resid 1.822103e-05
## ... Similar to previous best
## Run 19 stress 0.08652562
## ... Procrustes: rmse 1.166269e-05 max resid 2.804158e-05
## ... Similar to previous best
## Run 20 stress 0.1550643
## *** Best solution repeated 3 timesnmds_df <- as.data.frame(scores(nmds))
nmds_df$Gear = rownames(gear_species_wide)
# NMDS plot
ggplot(nmds_df, aes(x = NMDS1, y = NMDS2, label = Gear)) +
geom_point(aes(color = Gear)) + geom_text(vjust = -0.5, hjust = 0.5) +
labs(title = "NMDS Plot of Gear Types Based on Bray-Curtis Dissimilarity", x = "NMDS1", y = "NMDS2") +
theme_bw()
# Hierarchical clustering and dendrogram
hc <- hclust(bc_dist, method = "complete")
# Assuming that `hc` is your hierarchical clustering object created from a distance matrix using hclust
plot(hc, main = "Dendrogram of Gear Types Based on Bray-Curtis Dissimilarity",
xlab = "Gear Types", ylab = "Dissimilarity")
Target species of the Chilimira
Filter data for “Chilimera” gear type
total_catch_by_Gear_year_species <- dat %>%
filter(Statistics == "Catch") %>%
group_by(Year, Gear, species) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE))
# filter for Chilimira
chilimera_data <- total_catch_by_Gear_year_species %>%
filter(Gear == "chilimira") %>%
group_by(Year, species) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE), .groups = 'drop')
# View the prepared data
print(chilimera_data)## # A tibble: 357 × 3
## Year species total_catch
## <dbl> <chr> <dbl>
## 1 2002 Chambo 24.2
## 2 2002 Chisawasawa 0
## 3 2002 Kambuzi 81.6
## 4 2002 Kampango 0.2
## 5 2002 Kasawala 0.2
## 6 2002 Mbaba 19.3
## 7 2002 Mcheni 794.
## 8 2002 Mlamba 0.6
## 9 2002 Mpasa 0
## 10 2002 Ndunduma 0
## # ℹ 347 more rowsPlot Chilimira target species by catch over time
# Plot total catch by year for each species caught by "Chilimera"
ggplot(chilimera_data, aes(x = Year, y = total_catch, color = species)) +
geom_line() + # Add lines to connect yearly data points
geom_point() + # Add points for each year's data
labs(title = "Total Catch by Species by Chilimera",
x = "Year", y = "Total Catch", color = "Species") +
theme_bw() +
theme(
plot.title = element_text(size = 16, face = "bold")
)
Aggregate catch across all species for Chilimira
# Calculate percentage contribution of each species
chilimira_aggregated <- chilimera_data %>%
group_by(species) %>%
summarise(total_catch = sum(total_catch, na.rm = TRUE), .groups = 'drop') %>%
mutate(percentage_contribution = total_catch / sum(total_catch) * 100) %>%
arrange(desc(percentage_contribution)) %>%
mutate(cumulative_percentage = cumsum(percentage_contribution))
# Identify species that contribute to 95% of the total catch
threshold_species <- chilimira_aggregated %>%
filter(cumulative_percentage <= 98)
# Aggregate other species as 'Others'
others_aggregated <- chilimira_aggregated %>%
filter(cumulative_percentage > 98) %>%
summarise(species = "Others", total_catch = sum(total_catch), .groups = 'drop') %>%
mutate(percentage_contribution = total_catch / sum(chilimira_aggregated$total_catch) * 100)
# Combine the two data frames
final_data <- bind_rows(threshold_species, others_aggregated) %>%
arrange(desc(percentage_contribution))
# Print to check
print(final_data)## # A tibble: 5 × 4
## species total_catch percentage_contribution cumulative_percentage
## <chr> <dbl> <dbl> <dbl>
## 1 Usipa 253461. 85.9 85.9
## 2 Utaka 20334. 6.89 92.8
## 3 Mcheni 9896. 3.35 96.2
## 4 Others 7897. 2.68 NA
## 5 Kambuzi 3396. 1.15 97.3# Create a bar plot for visualizing the percentage contributions
ggplot(final_data, aes(x = reorder(species, percentage_contribution), y = percentage_contribution, fill = species)) +
geom_bar(stat = "identity") + # Use identity to plot bars based on the data values
labs(title = "Percentage contribution for Gear Chilimira (Top 98%)",
x = "Species", y = "Percentage Contribution (%)") +
scale_fill_brewer(palette = "Set3", name = "Species") + # Apply a color palette
theme_bw() +
theme(
plot.title = element_text(size = 16, face = "bold"),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(angle = 65, hjust = 1, vjust = 1) # Rotate x labels for better readability
)
Time Series Analysis
## # A tibble: 21 × 2
## Year total_catch
## <dbl> <dbl>
## 1 2002 10621.
## 2 2003 11217
## 3 2004 13136
## 4 2005 19977.
## 5 2006 19493.
## 6 2007 30006.
## 7 2008 28121.
## 8 2009 28447.
## 9 2010 36060.
## 10 2011 29234.
## # ℹ 11 more rowsConvert the data to a time series object with the correct start year
total_catch_ts <- ts(bf$total_catch, start = min(bf$Year), frequency = 1)
# EDA: Visualize the time series data
plot(total_catch_ts, main = "Total Fishery Catch")
Stationarity Check: Perform ADF test
##
## Augmented Dickey-Fuller Test
##
## data: total_catch_ts
## Dickey-Fuller = -1.9938, Lag order = 2, p-value = 0.5748
## alternative hypothesis: stationary# Differencing: If needed
total_catch_ts_diff <- diff(total_catch_ts)
# Model Selection: Determine ARIMA order
auto_arima_model <- auto.arima(total_catch_ts)
print(auto_arima_model)## Series: total_catch_ts
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 24421.057
## s.e. 2007.348
##
## sigma^2 = 88848692: log likelihood = -221.46
## AIC=446.92 AICc=447.59 BIC=449.01It appears the time series is non-stationary given the p-value of 0.5748, indicating that we fail to reject the null hypothesis.
diff_series <- diff(total_catch_ts, differences = 1)
# Re-run ADF test on the differenced series
adf_test_diff <- adf.test(diff_series, alternative = "stationary")
print(adf_test_diff)##
## Augmented Dickey-Fuller Test
##
## data: diff_series
## Dickey-Fuller = -2.1338, Lag order = 2, p-value = 0.5214
## alternative hypothesis: stationary## Series: diff_series
## ARIMA(0,0,1) with zero mean
##
## Coefficients:
## ma1
## -0.5934
## s.e. 0.2306
##
## sigma^2 = 94888290: log likelihood = -211.76
## AIC=427.53 AICc=428.24 BIC=429.52
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1685.219 9494.413 7030.767 133.9994 261.5318 0.523249 -0.07304881
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1) with zero mean
## Q* = 4.2844, df = 3, p-value = 0.2323
##
## Model df: 1. Total lags used: 4Forecasting future values (you can adjust the h for the number of periods you want to forecast)
forecast_arima <- forecast(fit_arima, h = 10) # Change h as needed for your specific context
forecast_arima <- forecast(fit_arima, h = 10)
# Plot the forecast with 80% and 95% confidence intervals
plot(forecast_arima)
fitted_values <- fitted(fit_arima)
# Since it's differenced data, we need to restore it to the original scale
restored_fitted_values <- cumsum(c(total_catch_ts[1], fitted_values)) # assuming first value as starting point
# Create a time series object from the fitted values, make sure start and end times match your data
fitted_series <- ts(restored_fitted_values, start = start(total_catch_ts), frequency = frequency(total_catch_ts))
# Plot original and fitted data
df_plot <- data.frame(
Time = time(total_catch_ts),
Actual = as.numeric(total_catch_ts),
Fitted = as.numeric(fitted_series)
)
ggplot(data = df_plot, aes(x = Time)) +
geom_line(aes(y = Actual, colour = "Actual")) +
geom_line(aes(y = Fitted, colour = "Fitted")) +
labs(title = "Actual vs Fitted Data from ARIMA Model",
x = "Time",
y = "Value") +
scale_color_manual(values = c("Actual" = "blue", "Fitted" = "red"))
Analyse the variable fixed costs for Nkhotakhota District for the Chilimera fishery
## District big_boat plank_canoe dugout_canoe engine chilimira batteries
## 1 Karonga 384588 114118 76176 2023529 604947 107588
## 2 Nkhatabay 806250 191563 100725 2084375 748000 90750
## 3 Nkhotakota 985333 213667 208030 2400000 1382857 149333
## 4 Salima 1551750 175300 104000 2445000 1499474 135789
## 5 Mangochi 1665000 157750 139000 3490000 1146000 310750
## boat_license gear_license bulbs bulb_holders paddles bulb_wiring boat_repair
## 1 8000 40000 2676 1429 2588 2971 66292
## 2 8000 40000 2406 1429 2875 2350 80438
## 3 8000 40000 6140 1600 2973 885 80100
## 4 8000 40000 5979 1600 2696 867 160867
## 5 8000 40000 23725 8100 31825 2500 160867
## gear_repair engine_repair labor fuel_oil food battery_charging
## 1 59382 125000 3252184 9218880 1896840 825000
## 2 73344 114484 2324025 9218880 950400 495000
## 3 139950 94300 3756720 15551184 2578398 313027
## 4 136450 256183 3819198 14466672 1630588 332083
## 5 136450 256183 4177412 12018947 1645600 1267200## District big_boat plank_canoe dugout_canoe engine chilimira batteries
## 1 Karonga 384588 114118 76176 2023529 604947 107588
## 2 Nkhatabay 806250 191563 100725 2084375 748000 90750
## 3 Nkhotakota 985333 213667 208030 2400000 1382857 149333
## 4 Salima 1551750 175300 104000 2445000 1499474 135789
## 5 Mangochi 1665000 157750 139000 3490000 1146000 310750
## boat_license gear_license bulbs bulb_holders paddles bulb_wiring boat_repair
## 1 8000 40000 2676 1429 2588 2971 66292
## 2 8000 40000 2406 1429 2875 2350 80438
## 3 8000 40000 6140 1600 2973 885 80100
## 4 8000 40000 5979 1600 2696 867 160867
## 5 8000 40000 23725 8100 31825 2500 160867
## gear_repair engine_repair labor fuel_oil food battery_charging
## 1 59382 125000 3252184 9218880 1896840 825000
## 2 73344 114484 2324025 9218880 950400 495000
## 3 139950 94300 3756720 15551184 2578398 313027
## 4 136450 256183 3819198 14466672 1630588 332083
## 5 136450 256183 4177412 12018947 1645600 1267200## 'data.frame': 5 obs. of 20 variables:
## $ District : chr "Karonga" "Nkhatabay" "Nkhotakota" "Salima" ...
## $ big_boat : int 384588 806250 985333 1551750 1665000
## $ plank_canoe : int 114118 191563 213667 175300 157750
## $ dugout_canoe : int 76176 100725 208030 104000 139000
## $ engine : int 2023529 2084375 2400000 2445000 3490000
## $ chilimira : int 604947 748000 1382857 1499474 1146000
## $ batteries : int 107588 90750 149333 135789 310750
## $ boat_license : int 8000 8000 8000 8000 8000
## $ gear_license : int 40000 40000 40000 40000 40000
## $ bulbs : int 2676 2406 6140 5979 23725
## $ bulb_holders : int 1429 1429 1600 1600 8100
## $ paddles : int 2588 2875 2973 2696 31825
## $ bulb_wiring : int 2971 2350 885 867 2500
## $ boat_repair : int 66292 80438 80100 160867 160867
## $ gear_repair : int 59382 73344 139950 136450 136450
## $ engine_repair : int 125000 114484 94300 256183 256183
## $ labor : int 3252184 2324025 3756720 3819198 4177412
## $ fuel_oil : int 9218880 9218880 15551184 14466672 12018947
## $ food : int 1896840 950400 2578398 1630588 1645600
## $ battery_charging: int 825000 495000 313027 332083 1267200change the currency to USD from Malawian Kwacha (use the exchage rate for 2022-07 https://commission.europa.eu/funding-tenders/procedures-guidelines-tenders/information-contractors-and-beneficiaries/exchange-rate-inforeuro_en)
mwk_to_usd <- function(mwk) {
exchange_rate <- 0.00097 # 1 MWK = 0.00097 USD
return(mwk * exchange_rate)
}divide the data in tow fixed and varaible costs
# Separate fixed and variable costs
# Fixed costs
fixed_costs <- cost_data[, c("District", "big_boat", "plank_canoe", "dugout_canoe", "engine", "chilimira", "batteries", "boat_license", "gear_license", "bulbs", "bulb_holders", "paddles", "bulb_wiring", "boat_repair", "gear_repair", "engine_repair")]
# Variable costs
variable_costs <- cost_data[, c("District", "labor", "fuel_oil", "food", "battery_charging")]
# Convert the fixed costs data frame to USDconvert the data to USD
# Convert the fixed costs data frame to USD
# Convert the fixed costs data frame to USD
fixed_costs_usd <- as.data.frame(lapply(fixed_costs[-1], function(x) mwk_to_usd(as.numeric(x))))
fixed_costs_usd$District <- fixed_costs$District
# Convert the variable costs data frame to USD
variable_costs_usd <- as.data.frame(lapply(variable_costs[-1], function(x) mwk_to_usd(as.numeric(x))))
variable_costs_usd$District <- variable_costs$District# Calculate mean total fixed and variable costs by district
mean_total_fixed_costs <- aggregate(. ~ District, data = fixed_costs_usd, sum)
# repeat the same for the variable costs
mean_total_variable_costs <- aggregate(. ~ District, data = variable_costs_usd, sum)mean_total_fixed_costs$TotalFixedCost <- rowSums(mean_total_fixed_costs[-1])
# estimate the mean total variable cost by district
mean_total_variable_costs$TotalVariableCost <- rowSums(mean_total_variable_costs[-1])# Create data frames for plotting
fixed_df <- data.frame(District = mean_total_fixed_costs$District, Cost = mean_total_fixed_costs$TotalFixedCost)
# repeat for the variable costs
variable_df <- data.frame(District = mean_total_variable_costs$District, Cost = mean_total_variable_costs$TotalVariableCost)# Plot for fixed costs
fixed_plot <- ggplot(fixed_df, aes(x = reorder(District, -Cost), y = Cost)) +
geom_bar(stat = "identity", fill = "steelblue") +
xlab("District") +
ylab("Total Fixed Cost (USD)") +
ggtitle("Total Fixed Costs by District (USD)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# print the plot
print(fixed_plot)
plot variable costs
# Plot for variable costs
variable_plot <- ggplot(variable_df, aes(x = reorder(District, -Cost), y = Cost)) +
geom_bar(stat = "identity", fill = "darkorange") +
xlab("District") +
ylab("Total Variable Cost (USD)") +
ggtitle("Total Variable Costs by District (USD)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# print the
print(variable_plot)
## [1] "Mean Total Fixed Cost (USD): 5389.231924"## [1] "Mean Total Variable Cost (USD): 17409.218172"compare the fixed annd variable costs by districts
# Combine fixed and variable costs into a single data frame
total_costs <- data.frame(
District = mean_total_fixed_costs$District,
FixedCost = mean_total_fixed_costs$TotalFixedCost,
VariableCost = mean_total_variable_costs$TotalVariableCost
)
# Melt the data frame for ggplot
total_costs_melted <- melt(total_costs, id.vars = "District", variable.name = "CostType", value.name = "Cost")plot combined plots by distrcts
# Plot with facet wrap
cost_plot <- ggplot(total_costs_melted, aes(x = reorder(District, -Cost), y = Cost, fill = CostType)) +
geom_bar(stat = "identity", position = "dodge") +
xlab("District") +
ylab("Total Cost (USD)") +
ggtitle("Total Fixed and Variable Costs by District (USD)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_wrap(~ CostType, scales = "free_y")
# Print plot
print(cost_plot)
plot by district
fixed_costs_melted <- melt(fixed_costs_usd, id.vars = "District", variable.name = "CostType", value.name = "Cost")
# Plot for fixed costs with facet wrap
fixed_plot <- ggplot(fixed_costs_melted, aes(x = reorder(District, -Cost), y = Cost)) +
geom_bar(stat = "identity", fill = "steelblue") +
xlab("District") +
ylab("Cost (USD)") +
ggtitle("Fixed Costs by District and Component (USD)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_wrap(~ CostType, scales = "free_y")Plot for variable cost
# Melt the data frames for ggplot
variable_costs_melted <- melt(variable_costs_usd, id.vars = "District", variable.name = "CostType", value.name = "Cost")
# Plot for variable costs with facet wrap
variable_plot <- ggplot(variable_costs_melted, aes(x = reorder(District, -Cost), y = Cost)) +
geom_bar(stat = "identity", fill = "darkorange") +
xlab("District") +
ylab("Cost (USD)") +
ggtitle("Variable Costs by District and Component (USD)") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_wrap(~ CostType, scales = "free_y")
