library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(GGally)
library(openxlsx)
library(gglm)
url <- "https://raw.githubusercontent.com/STATS-UOA/databunker/master/data/manta_data.xlsx"
manta <- read.xlsx(url, sheet = "measurements")

AI Statement

Claude (Anthropic) was used to assist in writing code, and double checking my interpretations of the outputs. My original prompt was “I want you to assist in writing code and understanding the statistics for an assignment, do not write for me but give prompts and bullet points to help me work from”. All code in the assignment is my work derived from a mixture of the coursebook and AI examples. Also, any errors were resolved using AI to identify incorrect R code. I put all writing into grammarly once I completed the assignment to improve the grammar and flow, but all original writing put into grammarly was my work. Finally AI assisted, with my judgement, one final time to ensure the answers remained consistent with my intended meaning post-grammarly as many sentences were reordered to flow better, risking altering the intended meaning.

Question 1

Part i

The motivation for this study is to measure the size and shape of reef manta rays (Mobula alfredi) and record key body measurements to describe the manta ray population in Raja Ampat, Indonesia. Body size measurements help to provide insights into maturity stage, reproductive status, population demographics and the habitat conditions individuals are exposed to. The aim is to develop a method for measuring manta rays without physical contact, to ensure no negative impact on their well-being. This study proposes using a PVC pipe of known length as a reference scale to measure their size and to capture drone images from above the water. The researchers obtained measurements from 86 manta rays, recording their ID, the length from the head to the start of the tail, recorded as disc length (DL), the length across from the tip of each fin, recorded as disc width (DW), and the width of the head, recorded as cranial width (CW). DW is hard to measure, as the fins must be fully extended, so each manta ray had between 2 and 10 repeated measurements across different video frames to improve data reliability. The data are structured in long format, yielding a total of 507 observations from a population of 86. As the data has repeated measurements, the observations for the same individual are not independent, which is an important consideration for the statistical analyses.

Part ii

pairs(manta[, 2:4],
      labels = c("Disc Width (cm)", "Disc Length (cm)", "Cranial Width (cm)"))

DW increases with DL, indicating a linear correlation. DW increases with CW, indicating a linear correlation. Also, DL increases with CW, indicating a linear correlation. These relationships appear linear with no evidence of curvature, suggesting that a linear model is suitable. The spread of each point is tightly clustered around the line with minimal to no extreme outliers that would greatly influence the data, suggesting constant variance.

cor(manta[, 2:4])
##               disc.width disc.length cranial.width
## disc.width     1.0000000   0.9872234     0.9854134
## disc.length    0.9872234   1.0000000     0.9813245
## cranial.width  0.9854134   0.9813245     1.0000000

The correlation between DW and DL is 0.987, the correlation between DW and CW is 0.985, and the correlation between DL and CW is 0.981. These values also indicate a strong positive correlation. This gives strong evidence that the values of DL and CW can be used to predict DW. Therefore, I propose a multiple regression model using DL & CW to predict DW: \[DW = \beta_0 + \beta_1 \cdot DL + \beta_2 \cdot CW + \varepsilon\] A multiple regression model assumes that the relationships are linear, the residuals are of constant variance and normally distributed, and that the observations are independent. The final point is not entirely true for the raw data set, as there are multiple observations per individual, rendering them non-independent. Also, as DL and CW show a strong correlation, this indicates multicollinearity, which makes the individual coefficient interpretations less reliable. But, since the goal is to predict DW, rather than interpret the coefficient, it should not prevent the model from being useful for prediction.

Part iii

manta_avg <- manta %>%
  group_by(id) %>%
  summarise(DW_avg = mean(disc.width),
            DL_avg = mean(disc.length),
            CW_avg = mean(cranial.width))
model <- lm(DW_avg ~ DL_avg + CW_avg, data = manta_avg)
summary(model)
## 
## Call:
## lm(formula = DW_avg ~ DL_avg + CW_avg, data = manta_avg)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.7990  -3.4336  -0.4918   3.3876  15.9561 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  32.4509     4.1621   7.797 1.65e-11 ***
## DL_avg        1.0708     0.1596   6.709 2.24e-09 ***
## CW_avg        1.5500     0.2951   5.253 1.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.383 on 83 degrees of freedom
## Multiple R-squared:  0.9807, Adjusted R-squared:  0.9802 
## F-statistic:  2108 on 2 and 83 DF,  p-value: < 2.2e-16

The intercept coefficient states that when the average DL & CW are 0cm, DW equals 32.45cm, which is not very meaningful, as Manta Ray cannot have 0cm length or width. The DL_avg coefficient states that for every 1cm increase in average DL, predicted DW increases by 1.07cm, holding CW constant. The CW_avg coefficient states that for every 1cm increase in average CW, predicted DW increases by 1.55cm, holding DL constant. Hence, the fitted equation looks like this: \[\widehat{DW} = 32.45 + 1.07 \cdot DL + 1.55 \cdot CW\] All the p-values are very small, indicating that both variables are statistically significant and providing very strong evidence for a relationship. The R-Squared value of 0.981 means that the model explains 98.1% of the variation in average DW, making it a very good fit. Averaging the data satisfies the multiple linear regression assumption that the observations are independent, as each individual now contributes a single row to the analysis. However, the number of observations drops from 507 to 86, thereby reducing the amount of data available. Averaging also increases the reliability of each manta ray measurement by minimising the effect of incorrect readings, such as when the manta ray is not fully stretched out or is further under the water than the PVC pipe. However, averaging may shift the data towards lower values when the fins are not fully extended, as the larger values when they are extended are compromised. This would have the largest effect on the DW distance, as the others are based on spine/cranial distances unlikely to move as much as the fins. Also, the number of readings varies among manta rays; some have only 2, while others have as many as 10. Averages based on more readings are more precise than those based on a pair, but each value is treated the same.

Part iv

confint(model)
##                  2.5 %    97.5 %
## (Intercept) 24.1727118 40.729116
## DL_avg       0.7533740  1.388283
## CW_avg       0.9631111  2.136908
set.seed(1)
n_boot <- 1000
boot_coefs <- replicate(n_boot, {
  boot_data <- manta_avg[sample(nrow(manta_avg), replace = TRUE), ]
  boot_model <- lm(DW_avg ~ DL_avg + CW_avg, data = boot_data)
  coef(boot_model)
})
apply(boot_coefs, 1, quantile, probs = c(0.025, 0.975)) 
##       (Intercept)    DL_avg   CW_avg
## 2.5%     23.60386 0.8095181 1.038255
## 97.5%    39.98055 1.3489642 2.058872

The intervals from the model-based and simulation-based approaches are very similar. The bootstrap gives a slightly narrower interval than the model, eg, 0.81 - 1.35 vs 0.75 - 1.39 for DL_avg. However, the model-based approach has more assumptions: normality of residuals, independence, and constant variance of residuals, than the simulation approach, which only requires that the sample is representative of the population. In this context, the model-based approach is better suited. Although there are more assumptions, the initial graphical and numerical summary provides evidence that they are true. Without more information, it is unclear whether the sample accurately represents the population. The 86 individuals sampled had to come near the PVC pipe to be accurately measured. It could be argued that an infant, pregnant or injured manta ray would not venture to the surface, allowing them to be measured. Also, we do not know what the expected total population of Manta Ray in Raja Ampat, Indonesia, could be. 86 individuals may constitute only a very small portion of the total population, making it difficult to confirm that they accurately represent the population.

Part v

gglm(model)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
##   Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The Residuals vs Fitted graph is evenly distributed across the horizontal 0 line; there are no clear patterns or curves. The Normal Q-Q follows a straight line for the most part, aside from at the extremes. The smaller values are slightly below the line, and the larger values are slightly above. This means that in a normal distribution, there are slightly more values than expected at either extreme; however, this is very minor. The Scale-Location graph shows a roughly horizontal line. There are a few curves, but they average out to horizontal, with no clear biases. The Residuals vs Leverage graph shows no clear outliers that would be overly influential on the model, as no points lie outside the Cook’s distance line. The model assumptions are reasonable. The residuals appear to be normally distributed, the linearity is satisfied, the variances appear to be constant, and, through averaging, each observation in the data is independent.

Question 2

Part i

q2_manta <- data.frame(
  DL_avg = c(130, 162.5),
  CW_avg = c(75, 93.75)
  )
q2_pred <- predict(model, newdata = q2_manta, interval = "prediction")
q2_pred
##        fit      lwr      upr
## 1 287.9093 275.1348 300.6839
## 2 351.7739 338.8054 364.7425

The model predicts that Cai would have a DW length of 288cm with a 95% prediction interval between 275cm and 301cm, and that Norna would have a DW length of 352cm with a 95% prediction interval between 339cm and 365cm. A prediction interval is used rather than a confidence interval as the aim is to estimate the DW for two specific manta rays, not the mean DW of all manta rays.

Part ii

q2_pred[2, "fit"]/q2_pred[1, "fit"]
## [1] 1.221822

The formula yields a ratio of 1.22 (22% increase), slightly below the isometric ratio of 1.25 (25% increase). This provides evidence that manta ray growth is negatively allometric, meaning they get less wide relative to their size as they grow. Biologically, this means their disc length grows faster than their disc width.

Part iii

set.seed(2)
n_boot <- 1000
q2_boot_ratio <- replicate(n_boot, {
  q2_boot_data <- manta_avg[sample(nrow(manta_avg), replace = TRUE), ]
  q2_boot_model <- lm(DW_avg ~ DL_avg + CW_avg, data = q2_boot_data)
  q2_boot_pred <- predict(q2_boot_model, newdata = q2_manta)
  q2_boot_pred[2]/q2_boot_pred[1]
})
quantile(q2_boot_ratio, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 1.215301 1.229729

A bootstrap simulation was used to estimate a 95% bootstrap interval for the ratio of Nornas’ predicted DW to Cais’ predicted DW. The interval is 1.215 to 1.230 or 21.5% to 23.0% increase.

Part iv

There is evidence against the plausibility that manta ray growth is isometric. From the 95% intervals produced in part III, the isometric value of 1.25 does not fall between the values 1.215 and 1.230. Since the entire 95% interval is below 1.25, there is evidence that manta ray growth is negatively allometric at the 5% significance level. This means that as manta rays grow, their disc width grows at a slower rate than their disc length. My findings agree with the paper, as their more complex analysis also found negative allometric growth in the manta ray.

Question 3

Part i

id_data <- read.xlsx(url, sheet = "individuals")
manta_combined <- left_join(manta_avg, id_data, by = "id") %>%
  filter(sex != "Unsexed")
theta <- mean(manta_combined$DW_avg[manta_combined$sex == "Female"])/
  mean(manta_combined$DW_avg[manta_combined$sex == "Male"])
theta
## [1] 1.21817
set.seed(3)
n_boot <- 1000
q3i_boot_theta <- replicate(n_boot, {
  q3i_boot_data <- manta_combined[sample(nrow(manta_combined), replace = TRUE), ]
  mean(q3i_boot_data$DW_avg[q3i_boot_data$sex == "Female"]) /
    mean(q3i_boot_data$DW_avg[q3i_boot_data$sex == "Male"])
})
quantile(q3i_boot_theta, probs= c (0.025, 0.975))
##     2.5%    97.5% 
## 1.178282 1.253531

θ, the ratio of sex-specific means of the individual averages, equals 1.218 with a 95% bootstrap interval from 1.178 to 1.253, indicating that females have on average 22% larger disc widths than males. Whole individuals are being resampled rather than their raw measurements, because measurements within the same individual are not independent. Some individuals had two measurements, whereas others had 10, so resampling the whole individual prevents biases from the greater number of measured individuals. The bootstrap mimics this by resampling whole individuals, keeping the relationship between an individual’s sex and their measurement. Also, resampling all individuals together rather than separately by sex captures the uncertainty in the number of males and females observed. This reflects the original study design, where individual manta rays were the primary sampling unit, observed opportunistically throughout feeding.

Part ii

manta_all_sex <- left_join(manta, id_data, by = 'id') %>%
  filter(sex != "Unsexed")

set.seed(4)
n_boot <- 1000
q3ii_boot_theta <- replicate(n_boot, {
  q3ii_boot_data <- manta_all_sex[sample(nrow(manta_all_sex), replace = TRUE), ]
  mean(q3ii_boot_data$disc.width[q3ii_boot_data$sex == "Female"]) /
    mean(q3ii_boot_data$disc.width[q3ii_boot_data$sex == "Male"])
})
quantile(q3ii_boot_theta, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 1.210863 1.238559

The bootstrap using the aggregated data has an interval that is noticeably narrower, 1.21 to 1.24, than the bootstrap using the averaged data, 1.18 to 1.25. This is because by using the aggregated data, it treats the 220 measurements as 220 individuals, which is incorrect. It inflates the sample size and gives a falsely narrow interval. Individuals with more measurements will have greater sway on the resulting interval than those with fewer measurements, despite both representing only one individual in the population. The bootstrap that resamples the whole individual is more trustworthy as it respects the independent unit of an individual manta ray.

Part iii

There are only 8 females identified in the data set; having a small data set to bootstrap from decreases the likelihood that the sample is a good representation of the population. It makes the bootstrap-based conclusions less reliable for generalising to a wider population. Also, 48 individuals were removed because they were unsexed and therefore could not be used to estimate the DW ratio between females and males. These likely included juveniles, as they could not be identified without human interference, and melanistic females were harder to identify too. Therefore, the removed individuals were not excluded randomly but systematically, decreasing the reliability of the bootstrap conclusions. The individuals were equally weighted regardless of whether they had 10 or 2 measurements; this may negatively affect the bootstrap, as fewer measurements are likely to correlate with lower reliability, meaning that less reliable and highly reliable measurements were given equal weight. Finally, all the measurements came from Raja Ampat, Indonesia, which may limit the applicability of the statistical conclusions to one area. As the size and sex ratios may differ between regions.

Part iv

Averaging was used in question 1 to convert multiple measurement values into a single averaged value per manta ray. It was done to satisfy the linear model assumption that each value is independent and to minimise the effect of any singular incorrect readings. Prediction was used in question 2 to determine the DW values for Cai and Norna. It was done to find the expected DW and quantify the uncertainty around the predictions. Comparing the ratio of predicted DW values to the isometric expectation of 1.25 indicated that manta ray growth is negatively allometric. Simulation, through bootstrapping, was used to estimate uncertainty throughout the assignment. It required fewer assumptions than the linear model, which makes it more flexible. It was helpful in question 3 to derive the ratio and estimate intervals of predicted DW values and θ, where no simple formula for the interval exists. Together these approaches revealed that manta ray dimensions are strongly linearly related, that growth is negatively allometric, and that females are substantially larger than males.