A Comprehensive Approach through Robust Regression and Gaussian/Mixed-Markov Graphical Models on the Example of Maritime Transportation Accidents: Evidence from a Listed-in-NYSE Shipping Company

Abstract

The main objective of this article is to determine the internal factors of maritime transportation accidents using a comprehensive approach through robust regression and Gaussian/mixed-Markov graphical models. Globally, this could be a strong incentive for the employees to negotiate higher compensation and for the insurance companies to impose higher premiums to cover the risk for these kinds of accidents. The article uses a dataset consisting of 166 real cases (human injuries) in the period 2014–2022 in different ships owned by a shipping company indexed in the New York Stock Exchange. The results of the study support the hypotheses as have been set in the article, connecting the internal factors with the injuries of any type. The practical implementation of the study is its ability to be used by policy makers in shipping to compensate employees depending on the risk of their work on board and at the same time to calculate the insurance premiums in a more accurate way. The originality of the research lies in the fact that this is a unique study in maritime transportation related to human accidents and not on ship or cargo casualties. The idea came from the results of another study conducted on a bibliometric analysis of the factors related to maritime transportation accidents. The findings of the current study can provide valuable insights to stakeholders and shipping planners in formulating effective policies for better wage packages and insurance premiums.

1. Introduction

Maritime transportation accidents (MTA) are distinguished in several types depending on the factors associated with them. The internal factors are related with causes on board, either “on-duty” or “off duty”, such as the nationality of the employee, the work location, the rank of the injured person, the working period in the same position, etc., and the external factors which are related to the sea trip such as maritime disasters, weather conditions, mechanical failures, collisions, etc. (Roberts et al. 2014).

Mental health and sickness of any kind cannot be considered as MTA because they are not related to an accident. COVID-19 disease cases on board must be separated from MTA because they are related to a sudden phenomenon and not to working conditions on board (Shan 2021), although as Lefkowitz and Slade (2019) stated, the COVID-19 pandemic embedded depression, anxiety or even suicide among maritime workers, and was associated with a high rate of insecurity, increasing the possibilities of being injured on board.

Based on a bibliometric analysis conducted by the same authors, there is limited literature highlighting the issue of MTA worldwide. The research is spread among few countries, with limited interconnections among researchers and few citations (Zampeta and Chondrokoukis, forthcoming).

The present research aims to distinguish the maritime transportation accidents with human casualties based on a dataset of 166 cases in the period 2014–2022 in different ships owned by a shipping company indexed in the NYSE. It describes the accidents taking into consideration the internal factors as described above, while the external factors are not considered in this article since they are the subject matter of a future research by the same authors.

The results could be used by the employees working on board for better wages and by the insurance companies for higher-risk premiums. The study has been conducted using advance econometric methodologies such as robust regression modelling and the Gaussian/mixed-Markov graphical models.

The results, as they are discussed in more detail in Section 4 and Section 5, are in accordance with a previous study which shows that the main internal factors behind an MTA are the following: the nationality of the employee, the rank of the injured person, the working deck, and the years of employment (Zampeta and Chondrokoukis 2022).

2. Literature Review

World trade is dependent on shipping (Castells 2011). Shipping is one of the most globalised industries in the world economy and is the leading means of transport (Ljung 2010; Tang and Gekara 2018). In total, about 80 percent of foreign trade is made by marine transport (European Union 2009). Shipping is a highly international, multicultural, and technological industry, and it faces strong demands on economic efficiency and profitability (Hanzu-Pazara et al. 2010; Ljung 2010). This has led to a globalized labour market of seafarers and to ship crews that are more and more multinational.

Multiculturalism is a general feature of crews today, and in this, language play a crucial role (Silos et al. 2012). About 70–80% of the world’s merchant fleet has multicultural crews (Magramo and Cellada 2009; Pyne and Koester 2005). Multicultural crews and a possible lack of a common language have produced a rising worry of the competence of ship crews. The worry of maritime safety has caused a growing demand for research in what kind of competences the crews operating the seas have. The question is inevitable, especially when it concerns areas with a high risk of accidents.

Maritime transportation accidents (MTA) and maritime emissions are gaining importance in recent years because of the upward trend on these types of incidents (Hussain et al. 2022).

Harrald et al. (1998) stated that human error is cited as the predominant cause of maritime transportation accidents. They described the modelling of human-error-related accident event sequences in a risk assessment of maritime oil transportation in Prince William Sound, Alaska. The risk analysts were confronted with incomplete and misleading data that made it difficult to use theoretical frameworks.

Maritime safety has been a core subject in maritime studies because it is coupled with transport safety, shipping efficiency, distribution reliability and loss prevention. Maritime accidents have often been attributed to human error, and discussion of human error and maritime accidents can be found in Millar (1980).

Havold (2010) demonstrated the use of safety culture to improve safety in the maritime environment. Chin-Shan and Chaur-Luh (2008) analysed the safety climate in the container shipping. Analyses of maritime accidents often ignore the link to human error.

Celik and Cebi (2009) proposed an analytical framework for identifying human errors in shipping accidents. The benefit of having an analytical framework is to provide a consistent manipulation of data and information of shipping accidents. They found that the primary root causes of shipping accidents are skill-based human errors and the shortfalls of execution of organizational processes.

Yip et al. (2015) provided empirical evidence that the maritime safety can be improved via the training of vessel crew members. Unver and Kocatepe (2019) and Akpinar and Sahin (2019) analysed the failures in maritime sector by providing detailed root causes.

Ship accidents are caused by various types of failures, e.g., deck officer error (26%), equipment failure (9%), structural failure (9%), crew error (17%), and mechanical failure (5%), among others. The factor that influences the risk level of maritime transportation accidents is defined as risk influence factor (RIF). To determine the risk factors of maritime transport, the latest related literature and maritime accident reports during 2012–2017 have been reviewed (Fan et al. 2020).

Since the United Kingdom Maritime and Coastguard Agency (UK MCA) proposed the formal safety assessment (FSA) framework to the International Maritime Organization, maritime accident risk models have been fast developed because of the goal-setting risk regime.

Based on the literature review above, the research hypotheses of the study is formulated as follows, while their test and discussion are presented in Section 4.4:

H₀₁. 

The mean ranks of work activity are the same across categories of parts of the body injured.

H₀₂. 

The mean ranks of work location are the same across categories of parts of the body injured.

H₀₃. 

The mean ranks of the number of months are the same across categories of parts of the body injured.

H₀₄. 

The mean ranks of the number of ranks are the same across categories of parts of the body injured.

3. Materials and Methods

3.1. Multiple Robust Regression

In least-squares techniques, one of the difficulties is that combinations of values of the explanatory variables can give some observations with far greater influence in the dependent variable than others. In 1970, the research on robust estimation of least squares provided new proposals that aimed at least for the protection against distortion by anomalous data and good efficiency when the data come from the ideal Gaussian model (Li 2006).

One of the progresses of robust regression is accelerating the analysis process by limiting the effect of some types of outliers and calling attention to unusual data. In the regression model, if there are unusual observations, they can sometimes severely distort estimates from regression with OLS (Andersen 2022). Although less common, unusual observations can also cause havoc for generalized linear models. This underscores the importance of detecting and properly handling outliers.

Unusual observations, albeit less frequent, can also ruin generalized linear models. This emphasizes how crucial it is to identify and manage outliers correctly. Additionally, robust estimators offer a powerful way to identify outliers or inconsistent substructures in data collection (Western 1995).

3.1.1. The Structure of Robust Regression Models

After the recognition that parametric models are rarely absolute, precise robust estimation of location has become an important tool (Kafadar 1983). Robust methods are based on the idea of redescending M-estimators. To clarify M-estimators, consider a classical simple linear regression model using the notation of Western (1995), as follows:

$$y=b_0+b_1{x}_1+\epsilon$$

(1)

Here, if we define $r_i=y_i-{\widehat{y}}_i$, the OLS technique aims to minimize the sum of r squared as follows:

$$Minimize{{\displaystyle \sum }}_i{r}_i^2$$

(2)

We named this method the least square residuals, and it is called $L_1$.

The objective function generally determines the shape of a statistic’s influence curve. If the median minimizes the sum of absolute residuals, the objection function is as follows:

$$Minimize{{\displaystyle \sum }}_i\left|r\right|$$

(3)

We named this method least absolute residuals, and it is noted as $L_2$.

One approach to robust regression involves devising estimators whose influence curves resemble the b_i weight. We use ρ for the weight and use a linear function instead of a quadratic one to make it simple; thus, we have the function as follows:

$$Minimize{{\displaystyle \sum }}_i\rho \left(r_i\right)$$

(4)

Robust estimation with the b_i weight uses an iterative weighted least squares algorithm. This algorithm consists of four steps (Western 1995):

➢

Obtain a set of start values from least squares or $L_1$ fit;

➢

If $r_i$ is a residual, $c$ is a tuning constant, and S is the current robust estimate of dispersion; scaled residuals will be created, which are represented by $u_i$ ($u_i=\frac{r_i}{cS}$);

➢

Form a set of weights $w_i$;

$w_i={(1-u_i^2)}^2$, for $u_i^2<1$, =0; otherwise;

➢

Estimate the model again with weighted least squares (WLS) using $w_i$;

➢

Update the weights with the residuals from the WLS fit, repeating until the coefficients show little change.

3.1.2. The Robust Regression Estimators

There are many robust regression estimators, such as L-estimator, M-estimator, MM-estimator, S Estimator, and others. The estimators are constructed on different bases; L-estimator is based on linear combinations of order statistics, M-estimator is extended from M-estimates of location by considering the size of the residuals, S-estimator minimizes a robust M-estimate of the residual scale, and MM-estimators build on both M-estimation and S-estimation to achieve a high breakdown point with high asymptotic efficiency (Andersen 2022).

When the linear regression given in matrix notation by $y=X\theta +\epsilon$ and the residuals are defined as $r_i=y_i-{\widehat{y}}_i$ for i = 1$\le i\le n$, and p is the number of independent variables, the estimators mentioned below are defined in the following paragraphs.

Edgeworth (1887) suggested a technique that involves reducing the sum of the residuals’ absolute values rather than the sum of their squares and the estimator ($L_1$) as follows:

$$\widehat{Q}{L}_1=arg \underset{Q}{\min }{{\displaystyle \sum }}_{i=1}^n\left|r_i\left(Q\right)\right|$$

(5)

Huber (1992) proposed an estimator which minimizes a function ρ of the errors rather than minimizing the sum of squared errors. The estimator keeps robustness with respect to vertical outliers and increases in Gaussian efficiency. The estimator (M) of Huber is as follows:

$${\widehat{Q}}_M=arg \underset{Q}{\min }{{\displaystyle \sum }}_{i=1}^{n}p\left|\frac{r_i\left(Q\right)}{\sigma }\right|$$

(6)

Rousseeuw and Yohai (1984) aim to find the smallest possible dispersion of the residuals and provide S estimator as follows:

$${\widehat{Q}}^S=arg \underset{Q}{\min }{\widehat{\sigma }}^S\left\{r_1\left(\theta \right),\dots ,r_n\left(\theta \right)\right\}$$

(7)

Yohai (1987) proposed MM estimation, which uses iteratively reweighted least squares (IRLS), the estimator is as follows:

$${\widehat{Q}}_{MM}=aarg \underset{Q}{\min }{{\displaystyle \sum }}_{i=1}^{n}p\left|\frac{r_i\left(Q\right)}{{\hat{\sigma }}^S}\right|$$

(8)

Robust estimators are not limited to what we gave above; all of them are proposed based on some pros and cons of the other estimators. In the empirical section, we obtained the results after processing the multiple regression models through robust regression with Huber and b_i-weight iterations. Huber (1992) proposed a class of estimators called M-estimators which satisfy three criteria (Huber 2011):

• reasonably efficient at the assumed model;
• large changes in a small part of the data or small changes in a large part of the data should cause only small changes in the result (resistant);
• gross deviations from the model should not severely decrease its efficiency (robust).

M-estimators are more sensitive to scaling and warn of possible problems in convergence. In the literature, there are debates on the practical usefulness of the b_i-weight and of redescending M-estimates in general (Kafadar 1983).

The method described above has been used in this study as one of the methods fulfilling the research requirements and the limitations regarding the sample dataset. The results are presented in Section 4.4.

3.2. Structural Equation Modelling (SEM)

In academic research, regression-based approaches, which are named first-generation techniques, are used to test hypotheses, such as multiple regression models, discriminant analysis, logistic regression, and ANOVA. These methods have three limitations, which restrict their applicability in some circumstances (Haenlein and Kaplan 2004). These three limitations are:

• In these models, we need one dependent and several independent variables, namely,
• the postulation of a simple model structure.
• In these kinds of models, we have the assumption that all variables can be considered observable.

The conjecture that all variables are measured without error.

Compared to real-life problems, many researchers have stated arguments about these limitations. Jacoby (1978) and Shugan (2002) address the issues of studying the impact of one or two variables and studying defined variables that imply omitting some aspect of reality. The authors would like to remark on the mediating or moderating effects that we do not have in regression-based approaches.

Regarding the second restriction, the theories on unobservable characteristics can only be considered once they have prior stand-alone validation, such as confirmatory factor analysis (Hair et al. 2021). The third limitation is well-known information from econometrics and statistics lectures; each observation has two errors, which are random error and systematic error.

Structural equation modelling (SEM), the second-generation technique, was proposed to overcome these limitations. Whereas regression-based approaches have only one dependent and many independent variables, SEM allows for simultaneous modelling with multiple dependent and independent variables. In SEM, researchers can use unobserved variables, and measurement errors take part in the model.

Structural Equation Modelling Approaches

There are two approaches to construct SEM, which are the covariance-based approach and the variance-based approach. Moreover, Westland (2019) categorized the products of SEM statistical analysis algorithms into three groups: pairwise canonical correlations between pairs, multivariate canonical correlation matrices and systems of regression approaches that fit data to networks of observable variables. Although the author stated that there are three groups, he mentioned the fourth category, developed by new social network analysis, which allows for both visualization and network-specific statistics.

If we present multiple regression analysis as follows:

$$Y_1=X_1+X_2+X_3+\cdots +X_n$$

(9)

Then, the structural equation modelling is presented below (Hair et al. 2021)

$$\begin{gathered} Y_1=X_{11}+X_{12}+X_{13}+\cdots +X_{1n} \\ {Y}_2=X_{21}+X_{22}+X_{23}+\cdots +X_{2n} \\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \\ {Y}_m=X_{31}+X_{32}+X_{33}+\cdots +X_{3n} \end{gathered}$$

(10)

In the literature, mainly two models are used: covariance-based SEM (CB-SEM) and partial least squares SEM (PLS-SEM, also called PLS path modelling) (Westland 2019). Path models can then be represented graphically by a path diagram (also called an arrow scheme). These diagrams present the relationship among the variables visually. A PLS path model consists of two elements, which are the structural model and measurement model. The structural model displays the relationships (paths) between the constructs. The measurement model specifies the relationships among observed variables underlying the latent variables. Using standard notations (Bollen 1989; Stein et al. 2012), we represent general SEM by the following equations.

$$\begin{gathered} x_1={\lambda }_1{\xi }_1+{\delta }_1 x_2={\lambda }_2{\xi }_1+{\delta }_2 x_3={\lambda }_3{\xi }_1+{\delta }_3 \\ {y}_1={\lambda }_3{\eta }_1+{\epsilon }_1 y_2={\lambda }_4{\eta }_1+{\epsilon }_2 y_3={\lambda }_5{\eta }_1+{\epsilon }_3 \\ {\eta }_1={\gamma }_{11}{\xi }_1+\zeta \end{gathered}$$

(11)

where $x_i$ and $y_i$ are observed indicators for latent variables, the ${\xi }_1$ and ${\eta }_1$ are latent variables, the ${\lambda }_i$ are factor loadings, the ${\delta }_i$ and ${\epsilon }_i$ are error terms, and the covariance between error terms is zero.

Traditional regression procedures are robust when it comes to measurement errors in the outcome but not in the predictors. Additionally, the relationship between error terms for two independent outcomes cannot be modelled using univariate regression techniques. SEM gives us the ability to model measurement error for both the predictor and the outcome. The method described above has been used in this study as one of the methods fulfilling the research requirements and the limitations regarding the sample dataset. The results are presented in Section 4.4.

3.3. Gaussian/Mixed-Markov Graphical Models (GGMs, MGMs)

Graphical models are used to designate relationships among a set of variables (Wermuth and Cox 2015). Graphical models bring together graph theory and probability theory for multivariate statistical modelling in a potent formalism (Wainwright and Jordan 2008). In these graphs, each variable is represented by a node, and any pair of nodes may become coupled, such as an edge. Edges represent corresponding conditional dependence; if the edges are missing, it means some form of conditional independence between the pair of variables. Edges can be drawn, directed or undirected, which show the direction of dependence of response on an explanatory variable and an equal standing, thus the edge between two variables.

One of the types of widely used graphical models are the graphical Markov models. Although graphical Markov models started to be developed after 1970 (Wermuth and Cox 2015), the history of the model started with research in genetics (Wright 1921), in physics (Gibbs 2010) and in probability theory (Markov 1912). Graphical Markov models are special subclasses of log-linear models for contingency tables and joint Gaussian distributions.

Altenbuchinger et al. (2020) defines Gaussian graphical models (GGMs) as tools to infer dependencies between biological variables with the assumptions of multivariate normal distributed data. Mixed graphical models (MGMs) can be a better choice if the data are not normally distributed. MGMs combine characteristics of Gaussian graphical and the Ising model. The Ising model uses discrete data.

MGMs are probabilistic graphical models which reflect the joint probability density function of a set of variables following two or more different data distributions. Here, if we give an example for a set of variables following two or more different data distributions, one set of variables distributed as a Gaussian, another set of variables distributed as a multinomial, we may have three different distributions, while one set may be a Poisson distribution.

A Gaussian graphical model (GGM) is a probability distribution. The distribution in n dimensions with p density is as follows (Kelner et al. 2020):

$$p_x\left(x\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^{n}det}}\exp (-{(x-\mu )}^T{{\displaystyle \sum }}_{}^{-1}\left(x-\mu \right)/2)$$

(12)

where μ is the mean and Ʃ is the covariance matrix. GGMs are one of the most widely used methods to model statistical relationships between observable variables in the natural and social sciences, machine learning, and other fields. In most of the settings in which GGMs are applied, the dimension is greater than the sample size (Kelner et al. 2020; Liu 2013).

The typical way of GGM estimation depends on regularized optimizations which depend on tuning parameters. If tuning parameters are large, they are powerless to find the edges with small weight; if the tuning parameters are small, they will generate many false edges resulting in high false discovery rates.

For the empirical analysis of the paper, GGMs and MGMs will be used because they have the benefit of producing reliable results regardless of the indicators’ measurement units or the kinds of variables employed. Full-order partial correlations are correlations between two variables corrected for all other variables under investigation. They make it possible to distinguish between direct and indirect effects. The foundation for estimating them is provided by Gaussian graphical models (GGMs) (Altenbuchinger et al. 2020; Bishop 2006).

The method described above has been used in this study as one of the methods fulfilling the research requirements and the limitations regarding the sample dataset. The results are presented in Section 5.

4. Research Results with Robust Regression Models

4.1. Descriptive Statistics

This study uses data based on research in the “2014–2022_Personal Accidents Metrics.xls” as presented in Appendix A. The aim of the analysis is to enhance and assess the main credentials influencing maritime workers’ health and the injuries they suffer in case of maritime accidents. The variables considered in this research are the following.

Nationality—there are some main nationalities but also some rare nationalities; for example, there are many “Filipino” but only one “Latvian”. Therefore, for Latvian, Romanian, etc., we generated one category: Central Eastern European (CEE).

Work activity level—we have both “Maintenance on Deck” and “Deck Maintenance”, and we merged them to create one category. This is why the study has new variables in the following analysis which are not in the research file. This is also one of the reasons why in this section, we define the variables; some came from the research file, and some are generated.

In the following section, we present tables and figures generated from the research file using the SPSS statistical package. Figures and tables generated from SPSS are more reliable than using the Excel routines.

The qualitative (categorical) variables and their categories are presented below:

Category:
1. FAC—First Aid Case; 2. LWC—Lost Workday Case; 3. MTC—Medical Treatment Cases; 4. Other; 5. RWC—Restricted Work Case.
Rank:
1. AB; 2. Bosun; 3. Cadet; 4. Cook; 5. Electrician; 6. Engineer; 7. Fitter; 8. Officer; 9. Oiler; 10. Ordinary Seaman (OS); 11. Pumpman; 12. Steward; 13. Wiper.
Nationality:	Nationality M (Nationality Merged):
1. Brazilian	1. Other
2. Filipino	2. Filipino
3. Georgian	3. Other
4. Greek	4. Greek
5. Hellenic	4. Greek
6. Latvian	5. CEE–Central Eastern European
7. Romanian	5. CEE–Central Eastern European
8. Russian	3. Other
9. Ukrainian	5. CEE–Central Eastern European
Work Location:
1. Accommodation; 2. Deck; 3. Engine; 4. Galley; 5. Other (Bridge, Cabin, Cargo Control Room, L/B DECK, Manifold, Pump Room, S/G Room, Workshop).
Work activity:
1. Deck Maintenance; 2. Deck Operation; 3. During Work; 4. Engine;
5. Engine Operation; 6. Mooring Operation; 7. Other.
Body part:
1. Arms; 2. Back; 3. Burns; 4. Chest; 5. Eye; 6. Feet; 7. Finger; 8. Hand; 9. Head; 10. Leg; 11. Others.
Body part Merged (Based on Panel Data Analysis):
1. Hand injuries, fingers, hand, wrist (HIFHW);
2. Foot injuries, ankle, knees, and legs (FAKL);
3. Body injury, back, chest, shoulder, ribs (BIBCSR) + burns + eye + head + other.

To determine the relationship between the body part of the injury and specific coordinates related to the maritime activity, such as rank, nationality of workers, work location, type of work activity and period on board, we calculate the correlation between variables and Chi-Square test. Because variables are nominal, we used a Cramer V value to see the correlation among them (Prematunga 2012).

The result shows that there is no correlation between the variables. For the second analysis, we use the Chi-Square independence test. Nevertheless, before the test, we must merge categories of rank to provide expected frequencies, which should be at least five for the majority (80%) of the cell’s requirement. The new rank of variables is defined below in

Table 1. Definition of new variables.

Old Value: Rank	New Value: New Rank
1. AB	1 Deck Dept
2. Bosun	1 Deck Dept
3. Cadet	1 Deck Dept
4. Cook	3 Catering/Steward Dept
5. Electrician	2 Engine Dept
6. Engineer	2 Engine Dept
7. Fitter	2 Engine Dept
8. Officer	1 Deck Dept
9. Oiler	2 Engine Dept
10. OS	1 Deck Dept
11. Pumpman	2 Engine Dept
12. Steward	3 Catering/Steward Dept
13. Wiper	2 Engine Dept

Source: own study.

. The work location of the positions presented above is shown in Figure 1.

4.2. Cross Tabulation

Because we have merged the variables to make them homogeneous and to avoid problems related to the sample size, it is required to consider the cross tabulation technique. Cross tabulation will give a more detailed view of the data. For example, the common nationality is Filipino, which is why the most common injured nationality is Filipino, but maybe the number of injuries within the Filipino population is less than in other nationalities. We can use crosstabs to investigate this kind of situation.

We have different ranks; if there is any pattern in rank to have the same injury, we can investigate this.

Table 2. Crosstab RankM/BodyPart.

		BodyPart											Total
		Arms	Back	Burns	Chest	Eye	Feet	Fingers	Hand	Head	Leg	Other
AB	Count	1	3	0	1	6	3	5	2	3	0	3	27
	%	3.7%	11.1%	NA	3.7%	22.2%	11.1%	18.5%	7.4%	11.1%	NA	11.1%	100.0%
Bosun	Count	0	1	0	0	0	0	3	0	4	3	1	12
	%	NA	8.3%	NA	NA	NA	NA	25.0%	NA	33.3%	25.0%	8.3%	100.0%
Cadet	Count	0	2	0	0	1	1	3	1	1	0	1	10
	%	NA	20.0%	NA	NA	10.0%	10.0%	30.0%	10.0%	10.0%	NA	10.0%	100.0%
Cook	Count	0	0	0	0	0	0	3	0	0	0	0	3
	%	NA	NA	NA	NA	NA	NA	100.0%	NA	NA	NA	NA	100.0%
Electrician	Count	0	1	1	0	0	0	3	1	2	0	1	9
	%	NA	11.1%	11.1%	NA	NA	NA	33.3%	11.1%	22.2%	NA	11.1%	100.0%
Engineer	Count	0	2	5	0	4	2	15	2	7	3	2	42
	%	NA	4.8%	11.9%	NA	9.5%	4.8%	35.7%	4.8%	16.7%	7.1%	4.8%	100.0%
Fitter	Count	0	0	0	0	1	0	0	0	1	1	1	4
	%	NA	NA	NA	NA	25.0%	NA	NA	NA	25.0%	25.0%	25.0%	100.0%
Officer	Count	0	0	1	0	1	2	1	1	0	1	2	9
	%	NA	NA	11.1%	NA	11.1%	22.2%	11.1%	11.1%	NA	11.1%	22.2%	100.0%
Oiler	Count	0	0	0	1	1	0	3	1	1	3	0	10
	%	NA	NA	NA	10.0%	10.0%	NA	30.0%	10.0%	10.0%	30.0%	NA	100.0%
OS	Count	1	2	0	0	5	2	6	1	1	2	1	21
	%	4.8%	9.5%	NA	NA	23.8%	9.5%	28.6%	4.8%	4.8%	9.5%	4.8%	100.0%
Pumpman	Count	0	2	0	0	1	0	0	0	0	2	0	5
	%	NA	40.0%	NA	NA	20.0%	NA	NA	NA	NA	40.0%	NA	100.0%
Steward	Count	0	1	0	0	0	0	1	1	1	0	0	4
	%	NA	25.0%	NA	NA	NA	NA	25.0%	25.0%	25.0%	NA	NA	100.0%
Wiper	Count	2	2	0	0	2	1	2	1	0	0	0	10
	%	20.0%	20.0%	NA	NA	20.0%	10.0%	20.0%	10.0%	NA	NA	NA	100.0%
	Count	4	16	7	2	22	11	45	11	21	15	12	166
	%	2.4%	9.6%	4.2%	1.2%	13.3%	6.6%	27.1%	6.6%	12.7%	9.0%	7.2%	100.0%

Source: own elaboration.

below shows that the most common injury is finger injury at 27%, and the most common crew which had an injury is engineer. Most common injuries within rank are:

➢

Eye (18%) for AB;

➢

Head (33%) for Bosun;

➢

Finger (30%) for Cadet;

➢

Finger (100%) for Cook;

➢

Finger (30%) for Electrician;

➢

Finger (15%) for Engineer;

➢

Finger (28%) for OS.

The rest of the ranks have one or two injuries with the same ratios within their injury rate in many cases, which is why we did not present the ratios of other ranks. In summary, although the most common injury is a finger injury, for AB, the common injury is the eye, and for Bosun, the common injury is the head.

Another question is if any rank has more accidents in any specific location.

Table 3. Cross tab RankM/WorkL.

			WorkL					Total
			Accommodation	Deck	Engine	Galley	Other
RankM	AB	Count	5	20	1	1	0	27
		%	18.5%	74.1%	3.7%	3.7%	NA	100.0%
	Bosun	Count	0	12	0	0	0	12
		%	NA	100.0%	NA	NA	NA	100.0%
	Cadet	Count	2	3	5	0	0	10
		%	20.0%	30.0%	50.0%	NA	NA	100.0%
	Cook	Count	0	0	0	2	1	3
		%	NA	NA	NA	66.7%	33.3%	100.0%
	Electrician	Count	0	2	7	0	0	9
		%	NA	22.2%	77.8%	NA	NA	100.0%
	Engineer	Count	3	7	26	3	3	42
		%	7.1%	16.7%	61.9%	7.1%	7.1%	100.0%
	Fitter	Count	0	1	3	0	0	4
		%	NA	25.0%	75.0%	NA	NA	100.0%
	Officer	Count	0	5	1	1	2	9
		%	NA	55.6%	11.1%	11.1%	22.2%	100.0%
	Oiler	Count	0	1	6	0	3	10
		%	NA	10.0%	60.0%	NA	30.0%	100.0%
	OS	Count	3	17	1	0	0	21
		%	14.3%	81.0%	4.8%	NA	NA	100.0%
	Pumpman	Count	0	2	0	0	3	5
		%	NA	40.0%	NA	NA	60.0%	100.0%
	Steward	Count	1	0	0	3	0	4
		%	25.0%	NA	NA	75.0%	NA	100.0%
	Wiper	Count	1	2	7	0	0	10
		%	10.0%	20.0%	70.0%	NA	NA	100.0%
Total		Count	15	72	57	10	12	166
		%	9.0%	43.4%	34.3%	6.0%	7.2%	100.0%

Source: own elaboration.

shows that the most common location of injuries is the deck. Within rank variables, the most common locations are:

➢

Deck (74%) for AB;

➢

Deck (100%) for Bosun;

➢

Galley (67%) for Cook;

➢

Engine (78%) for Electrician;

➢

Engine (62%) for Engineer;

➢

Engine (75%) for Fitter;

➢

Deck (56%) for Officer;

➢

Engine (60%) for Oiler;

➢

Deck (81%) for OS;

➢

Other (60%) for Pumpman;

➢

Galley (75%) for Steward;

➢

Engine (70%) for Wiper.

The results show that for different work expertise, there are different locations regarding most common injury. Generally, technical workers have accidents in the engine location and the deck for servants.

Based on the within percentages, we saw that some experts have more injuries in their respective places.

Table 4. Cross tab RankM/WorkL.

		WorkA							Total
		Deck Maintenance	Deck Operation	During Work	Engine Maintenance	Engine Operation	Mooring Operation	Other
AB	Count	3	7	7	1	0	3	6	27
	%	11.1%	25.9%	25.9%	3.7%	NA	11.1%	22.2%	100.0%
Bosun	Count	4	4	2	0	0	1	1	12
	%	33.3%	33.3%	16.7%	NA	NA	8.3%	8.3%	100.0%
Cadet	Count	1	1	3	2	1	0	2	10
	%	10.0%	10.0%	30.0%	20.0%	10.0%	NA	20.0%	100.0%
Cook	Count	0	0	0	0	0	0	3	3
	%	NA	NA	NA	NA	NA	NA	100.0%	100.0%
Electrician	Count	0	1	0	5	1	0	2	9
	%	NA	11.1%	NA	55.6%	11.1%	NA	22.2%	100.0%
Engineer	Count	0	2	8	14	6	2	10	42
	%	NA	4.8%	19.0%	33.3%	14.3%	4.8%	23.8%	100.0%
Fitter	Count	0	0	0	3	0	0	1	4
	%	NA	NA	NA	75.0%	NA	NA	25.0%	100.0%
Officer	Count	1	3	2	1	0	0	2	9
	%	11.1%	33.3%	22.2%	11.1%	NA	NA	22.2%	100.0%
Oiler	Count	1	0	2	2	4	0	1	10
	%	10.0%	NA	20.0%	20.0%	40.0%	NA	10.0%	100.0%
OS	Count	3	5	3	1	0	4	5	21
	%	14.3%	23.8%	14.3%	4.8%	NA	19.0%	23.8%	100.0%
Pumpman	Count	0	0	2	0	0	0	3	5
	%	NA	NA	40.0%	NA	NA	NA	60.0%	100.0%
Steward	Count	0	0	0	0	0	0	4	4
	%	NA	NA	NA	NA	NA	NA	100.0%	100.0%
Wiper	Count	0	0	1	5	2	0	2	10
	%	NA	NA	10.0%	50.0%	20.0%	NA	20.0%	100.0%
	Count	13	23	30	34	14	10	42	166
	%	7.8%	13.9%	18.1%	20.5%	8.4%	6.0%	25.3%	100.0%

Source: own elaboration.

shows the main activity of the ranks when they have an injury. Like work location cross-tabulation analysis, technical workers’ most common injury activities are the same, which is engine maintenance.

4.3. Research Questions and Variables Selection

In this section, we select the variables based on the previous analysis. In the previous analysis, body part injured is the dependent variable, and rank, nationality, work location, work activity and period on board in months (POBM) are the independent variables. We use all these independent variables except for Nationality because there is a high weight on Filipino nationality, which will affect the results. We aim to analyse if these independent variables significantly differentiate injured body parts. The research questions which will be analysed below are the following:

• work activity is the main coordinate that determines the injuries of maritime workers;
• work location determines an increase in body injuries of maritime workers;
• number of months spent on the ship (period on board) leads to a decrease in the number of injuries.

Moreover, in panel data models and SEM models, the dependent variable “parts of body injured” is used in different categories and with different estimation methods. In this section, we similarly use two different “parts of body injured” variables, one has 3 categories, and one has 23 categories.

To investigate whether there is a statistically significant difference between parts of the body injured, we select the non-parametric method to check the normality of the variables. The test’s null hypothesis is that the mean ranks of the groups are the same. Under this test, we can check only the existence of differences among the groups, but we cannot say which group is more important than the other.

The subsection results below show that the variable “body injuries” in three categories has a significant relationship between work activity and workplace, but not between rank and POBM. In 23 categories, the variable “injured parts of the body” has a relationship between workplace and sequence, but there is no significant relationship between work activity and POBM.

4.4. Research Hypotheses Testing and Discussion

Following the analysis as described above in the section on cross-tabulation (Section 4.2), we will proceed with the merged data on body part injuries and work activity. The dataset consists of 166 observations (

Table 5. Descriptive statistics: minimum and maximum.

	N	Minimum	Maximum
BodyPartM	166	1	3
Work Activity	166	1	73
Valid N (listwise)	166

Source: own elaboration.

), and the Kruskal–Wallis test for work activity is presented in

Table 6. Kruskal–Wallis test results for work activity.

Test Statistics a,b
	Work Activity
Chi-Square	5.174
df	2
Asymp. Sig.	0.075

^a Kruskal–Wallis Test; ^b Grouping Variable: BodyPartM. Source: own elaboration.

. The Kruskal–Wallis test, proposed by Kruskal and Wallis in 1952, is a nonparametric method for testing whether samples originate from the same distribution. The first research hypothesis is as follows:

H₀₁. 

The mean ranks of work activity are the same across all categories of parts of the body injured.

Table 6 shows that the null hypothesis is not rejected at the 10% significance level but is rejected at the 5% level (Asymp. Sig. 0.075). We conclude that work activity injuries are different across body part groups. Different work activities injure different body parts.

Table 7. Descriptive statistics: means.

Ranks
	BodyPartM	N	Mean Rank
Work Activity	BIBCSR	80	92.27
	FAKL	28	75.38
	HIFHW	58	75.33
	Total	166

Source: own elaboration.

presents the work activity and the ranks of the merged body parts with their means, while

Table 8. Hypothesis test summary.

Null Hypothesis	Test	Significance	Decision
The distribution of work activity is the same across categories of BodyPartM	Independent samples Kruskal–Wallis Test	0.075	Reject the null hypothesis

Note: Asymptotic significances are displayed. The significance level is 0.10. Source: own elaboration.

presents the hypothesis test summary as defined above in H₀₁.

To verify the results, we used the independent samples Kruskal–Wallis test, a nonparametric method for testing whether samples originate from the same distribution, as shown in Figure 2.

The second hypothesis is stated among work location and parts of the body injured as follows:

H₀₂. 

The mean ranks of work location are the same across categories of parts of the body injured.

Table 9. Descriptive statistics: minimum and maximum.

	N	Minimum	Maximum
BodyPartM	166	1	3
Work Location	166	1	13
Valid N (listwise)	166

Source: own elaboration.

presents the descriptive statistics of merged body parts and work location,

Table 10. Kruskal–Wallis test results for work location.

	Work Location
Chi-Square	7.417
df	2
Asymp. Sig.	0.025

Source: own elaboration.

shows that the null hypothesis is rejected at the 5% significance level. We concluded that work location injuries are different across body part groups.

Table 11. Descriptive statistics: means.

Ranks
	BodyPartM	N	Mean Rank
Work Location	BIBCSR	80	74.71
	FAKL	28	82.71
	HIFHW	58	96.01
	Total	166

Source: own elaboration.

presents the work location and the ranks of the merged body parts with their means, while

Table 12. Hypothesis test summary.

Hypothesis Test Summary
Null Hypothesis	Test	Significance	Decision
The distribution of work location is the same across categories of BodyPartM	Independent samples Kruskal–Wallis test	0.025	Reject the null hypothesis

Note: Asymptotic significances are displayed. The significance level is 0.10. Source: own elaboration.

presents the hypothesis test summary as defined above in H₀₂:

To verify the results, we used the independent samples Kruskal–Wallis test, a nonparametric method for testing whether samples originate from the same distribution as shown in Figure 3.

The third research hypothesis is stated among work period (number of months in the same position) and parts of the body injured as follows:

H₀₃. 

The mean ranks of the number of months are the same across categories of parts of the body injured.

Table 13. Kruskal–Wallis test results for period on board.

	Period on Board (Months)
Chi-Square	0.415
df	2
Asymp. Sig.	0.812

Source: own elaboration.

shows that the null hypothesis cannot be rejected. We conclude that the mean ranks of the number of working months in the same position are the same across categories of parts of the body injured (Asymp. Sig. statistic 0.812).

The fourth research hypothesis is stated among the number of ranks and parts of the body injured as follows:

H₀₄. 

The mean ranks of the number of rank are the same across categories of parts of the body injured.

Table 14. Descriptive statistics: minimum and maximum.

Descriptive Statistics
	N	Minimum	Maximum
BodyPartM	166	1	3
Rank	166	1	23
Valid N (listwise)	166

Source: own elaboration.

presents the descriptive statistics of the variables, while

Table 15. Kruskal–Wallis test results for rank.

	Period on Board (Months)
Chi-Square	0.41
df	2
Asymp. Sig.	0.810

Source: own elaboration.

shows that the null hypothesis cannot be rejected. We concluded that the rank is the same across body part groups.

5. Results of Gaussian and Mixed Markov Graphical Models (GGMs, MGMs)

The research endeavour is complemented by a network analysis performed through Gaussian and mixed-Markov graphical models (GGMs, MGMs) processed through partial correlations. GGMs and MGMs substantiate previous results and strengthen the robustness of the empirical analysis. The main purpose of these advanced modern econometric models is to evidence the existence and intensity of the connections between all variables in a comprehensive approach and to enhance the linkages between specific coordinates related to the maritime activity and workers injuries. The Gaussian graphical models (GGMs) are entailed in Figure 4, and the mixed-Markov models (MGMs) are presented in Figure 5.

The configuration of a GGM under the format of a network, which is usually titled a partial correlation network, presents positive partial correlations with blue edges, while negative partial correlations are entailed by red edges (Figure 4 and Figure 5). In addition, the absolute strength of a partial correlation is highlighted by the width and saturation of an edge (Epskamp et al. 2018). If there is no edge between two nodes/variables, it means that the partial correlation is zero and that those two variables are independent after conditioning on all other variables in the dataset. In this case, the GGM can be seen as a network model of conditional associations.

Both GGMs presented in Figure 4 and MGMs entailed in Figure 5 show that there is a very strong positive relation between work activity and parts of body injured (configured here as an encoded combined variable from hand, foot, and body injuries). At the same time, a positive correlation is revealed between work location and parts of body injured, work location also being strongly and positively linked with rank and nationality of maritime workers. Hand injuries, foot injuries and body injuries are positively correlated with work location and are inversely correlated with work activity, yet in the case of hand injuries, these linkages are more pronounced.

6. Conclusions

To comply with the requirements regarding sample size, number of variables and type of research, we have merged some of the variables as shown in Section 4.1. We have presented the cross-tabulation (Section 4.2) to determine the merged variables used in the analysis.

Robust regression presents the advantage of undertaking robust estimates, thus avoiding spurious regression and coping with possible outliers within the sample. The analysis above has shown the following.

The first hypothesis “The mean ranks of work activity are the same across categories of parts of the body injured” is rejected at the 5% significance level. Different work activities injure different body parts.

The second hypothesis “The mean ranks of work location are the same across categories of parts of the body injured” is also rejected at the 5% significance level. Work location injuries are different across body part groups.

The third hypothesis “The mean ranks of the number of months are the same across categories of parts of the body part injured” is not rejected. Number of months is an irrelevant variable to the dependent variable (body injuries).

The fourth hypothesis “The mean ranks of the number of ranks are the same across categories of parts of the body injured” is also not rejected. Rank is an irrelevant variable to the dependent variable (body injuries).

Therefore, this research has shown that work activities on board and the work location are important factors for body injuries, while the time of employment and the rank of the worker are irrelevant variables on body injuries. Compensation packages and insurance premiums must be different for these positions (work activities and work location on the ship).

Spurious correlation is also avoided with the use of Gaussian and mixed-Markov graphical models, which provide a comprehensive view of the interlinkages between all considered variables. These graphical models employed in the current research have shown strong connections between work activity, work location, rank, nationality, and the injuries suffered by maritime workers.

GGMs and MGMs also present the advantage of providing robust results regardless of the measurement units of indicators/type of variables used in the empirical analysis. Moreover, structural equation modelling strengthened the research endeavour and showed consistency of all results through robust estimates captured by the maximum likelihood procedure (MLE).

Overall, the findings of the current study can provide valuable insights to stakeholders, shipping companies, insurance companies and policy planners in formulating effective policies for insurance premiums, compensation packages and human resource evaluations.

Future studies could focus on analysing MTA in terms of external factors, as the present study referred only to internal factors, (one limitation of this study), as described in the introduction or/and in the mix of two, internal and external factors. It could also be interesting to analyse a bigger sample with ships from different shipping companies (another limitation of this study).