ML For Hatchery

Author

Nicole Steffer 

library(keras)
library(tensorflow)


library(dplyr)
library(caret)

tensorflow::set_random_seed(123)
set.seed(123)
rm(list = ls())
D0 <- readRDS("datacompile_2024-02-14-N.RDS")
predictors <- readRDS("predictors.rds")


predictors_numeric = sapply(D0[, predictors], is.numeric)
predictors_numeric = predictors[predictors_numeric]


D <- D0 %>%
    dplyr::select(all_of(c("Yield", predictors_numeric))) %>%
    na.omit() %>%
    scale() %>%
    as.matrix()

Data

We used oyster hatchery data gathered from The Horn Point Laboratory Oyster Hatchery (HPLOH) to examine how various observations impact oyster larvae. The observations included where the broodstock collection was collected from, conditioning of the water, its pH, salinity and average temperature. Spawning, which looked at average shell height for male and female, how many millions of eggs each female had (fecundity), the gonadal index and if there was stimulation used as well as larval culture when production was started and if the carbonate tank was buffered.
The data also includes environmental data that was collected from Goose’s Reef buoy which observed chlorophyll concentration, salinity, turbidity, water temperature, and dissolved oxygen concentration on the day the tanks were filled for the hatchery. This data also included the atmospheric conditions, air temperature, wind direction and wind speed. We considered a total of 144 observations.

See Figure 1 for the histogram of the response variable

hist(D0$Yield/100, xlab = "yield", ylab = "frequency", main = "")
Figure 1: Histogram

Model

We used a neural network model using Keras and Tensorflow in Rstudio to build our model. The hyperparameters we started with were using 60% of the data as the training set, we used the activation function “relu” with two dense layers, one at 40 the other 60 and an another activation layer of 1 with the activation function “linear” for our final output. We used a batch size of 6, epochs 40 and validation split of 0.2. Using these hyperparameters we did not get the results we desired so we changed the dense layers, batch size, epochs, validation splits and percent of data in the training set but continued to get an \(R^2\) around .08-.1. We changed the activation functions for the first two dense layers to “ELU” , increased the training set percentage and with trial and error we were able to settle on using 80% of the data for training, two dense layers at 40 and 60 with activation function “ELU” and one dense layer of 1 with an activation function of “linear”. Batch size 10, epochs 50 and validation split of .1 to obtain an \(R^2\) of .25 which was the best output we got.

n <- nrow(D) 
n_train <- round(0.8 * n) 
i_train <- sort(sample(n, n_train))

The sample size is 437. The size of the training set is 350.

Y_train = D[i_train, "Yield"]
Y_test = D[-i_train, "Yield"]
X_train = D[i_train, predictors_numeric]
X_test = D[-i_train, predictors_numeric]


model <- keras_model_sequential(input_shape = length(predictors_numeric)) %>%
    layer_dense(40, activation = "ELU") %>%
    layer_dense(60, activation = "ELU") %>%
    layer_dense(1, activation = "linear")


summary(model)
#> Model: "sequential"
#> ________________________________________________________________________________
#>  Layer (type)                       Output Shape                    Param #     
#> ================================================================================
#>  dense_2 (Dense)                    (None, 40)                      5680        
#>  dense_1 (Dense)                    (None, 60)                      2460        
#>  dense (Dense)                      (None, 1)                       61          
#> ================================================================================
#> Total params: 8201 (32.04 KB)
#> Trainable params: 8201 (32.04 KB)
#> Non-trainable params: 0 (0.00 Byte)
#> ________________________________________________________________________________

Train the model

model %>% compile(
    loss = "mean_squared_error", 
    optimizer = optimizer_adam(learning_rate = 0.001)
)
BS = 10
EP = 50
history <- model %>%
    fit(x = X_train,
        y = Y_train,
        batch_size = BS,
        epochs = EP,
        validation_split = 0.1,
        verbose = 0)


predictions <- predict(model, X_test)
#> 3/3 - 0s - 82ms/epoch - 27ms/step

Compute performance metrics on the testing set

Perf=caret::postResample(obs = Y_test,
                    pred = predictions)
Perf
#>      RMSE  Rsquared       MAE 
#> 1.1106142 0.2483886 0.8312196

See Figure 2 for the results

plot(x = Y_test, y = predictions)
Figure 2: Predictions vs Observed

Results

Summary: Using the data set, we developed a neural network to predict hatchery yield. We used 80% of the data in the training set which was a sample size of 350. We used the activation function “ELU” with two dense layers, 40 and 60, one activation function “linear” at 1. Batch size was 10, epochs 50 and validation splits .1 to achieve an \(R^2\) of 0.25.