A multifactorial study on duration of temporary disabilities in Spain

ABSTRACT

The extreme variability of temporary disability duration has a deep effect in public health. We tried to understand what factors duration of disability depends on. Through cohort study with data of temporary disabilities collected by Ibermutuamur from 2008 to 2012, we used statistical multivariate methods. The most reliable and convenient algorithm to predict duration was a categorical classification tree that distinguished between brief and long disabilities, taking into account both medical-biological and socioeconomic factors. The influence of socioeconomic factors in the disability process made numeric predictive models not accurate enough. Some of these socioeconomic factors were isolated and their influences were quantified. In particular, the one we named factor unemployment could explain a huge increase in duration for certain common diagnoses such as anxiety, low back pain, headache, and depression.

KEYWORDS:

• Classification tree
• multifactorial ANOVA
• prediction model
• temporary disability
• unemployment

Many efforts have been made in different countries to understand what factors explain duration of disability.^1–3 In this work, we analyze data on temporary disability (TD) due to common disease. These data were collected in Spain from 2008 to 2012 by Ibermutuamur, one of the main widespread companies (Mutuas) that cooperate with National Institute of Social Security in managing worker's compensation insurance and medical attention in case of temporary disability, both strictly work related (due to a work accident or a typical occupational disease) and not work related (due to a common disease or an accident not related to work). In fact, about 78% of temporary disabilities in Spain are managed by Mutuas, so that the data we studied could be considered to be representative in order to understand what factors affect duration of temporary disabilities in Spain. Besides biological and medical factors, we analyzed other factors included in the data (according to Ibermutuamur's data-gathering protocol). Some of them, such as socioeconomic and occupational factors, were considered to be relevant in previous papers.^4–5 In Spain, recent studies have been carried out from small samples.^6–9 In this study, we tried to take advantage of the huge size of data and powerful data mining in order to isolate the relevant factor as clearly as possible.

The study of this issue depends obviously on the country we focus on since we must take into account its specific rules. In this case, we point out that, in Spain,⁷ employees receive 60% of salary from the fourth day to the 20th day of TD, and 75% during the rest of TD. Moreover, money is paid by the company during the first 15 days, and it is paid (directly or indirectly) by Mutua from then on. Therefore, it is important to discriminate whether TD lasts more than 15 days. Self-employed workers have the same benefit but it is paid directly by Mutua from the fourth day until the end of TD.

Methods

Original data were composed of N = 1,085,824 entries, each one consisting of 40 values (most of them categorical) corresponding to a specific person who received medical release after a temporary medical leave that occurred between January 1, 2008 and December 31, 2012. The duration of TD (total number of days of medical leave) was the response variable. In other words, every entry consisted of one response variable and 39 predictors. Duration was also transformed into another binary categorical variable to distinguish whether duration was more than 15 days. We discriminated between two main types of predictors: medical-biological and socioeconomic. Prior to searching a suitable statistical model, data were simplified by considering just the most common diagnoses in workers according to Pareto rule (see Figure 1), thus reducing the total to N = 370,076. Moreover, some original predictors were removed after applying ANOVA model (as explained in the following), so that we eventually considered just 10 final predictors: gender, age, occupation, geographical zone (we distinguished between 4 different zones as heterogeneously as possible), main diagnosis (we considered just the 29 most common diagnoses in Figure 1), existence of codiagnosis, type of employment regime (employee or self-employed), payment (indirect if Ibermutuamur pays through company or direct if Ibermutuamur pays directly to the worker), civil status, and children (having or not). We can see in Figure 1 and Tables 1 and 2 further details about these variables. Other variables not considered in the study were, for example, number of different episodes during TD, smoking status, drinking status, medical background, town, type of contract, beginning date, ending date, and so on.

Figure 1. Distribution of the main diagnostics.

Table 1. Duration of temporary disability by diagnosis.

Diagnosis	N	Mean	SD	Median	Over 15 days (%)
Gastroenteritis	40,780	3.8	8.0	3	1.9
Noninfective gastroenteritis	16,342	4.5	11.5	3	2.4
Tonsillitis strep throat	10,374	4.8	5.7	4	2.2
Acute pharyngitis	29,449	4.8	8.3	3	2.5
Tonsillitis not strep throat	16,373	5.4	6.5	4	2.7
Cold	26,171	5.4	9.0	4	3.5
Flu	24,169	6.3	7.5	5	3.9
Diarrhea	3,909	6.4	19.4	3	4.4
Others flu types	7,137	6.6	7.7	5	4.2
Bronchitis not acute or chronic	3,461	10.8	20.6	7	13.4
Migraine	3,202	10.8	30.1	2	11.7
Acute or chronic bronchitis	11,482	11.7	20.6	8	15.9
Headache	3,622	17.8	40.2	4	20.6
Renal colic	7,850	20.0	36.0	8	29.0
Ankle strain	6,571	24.1	30.7	15	48.7
Giddiness	9,379	24.7	46.8	8	30.2
Dorsalgia	4,144	32.1	49.8	11	42.5
Low back pain	48,933	35.2	49.9	15	49.7
Cervical pain	17,886	50.5	58.0	31	64.5
Inguinal hernia	6,044	52.8	41.5	43	95.0
Other maternal disorder related to pregnancy	3,822	58.2	46.6	46	88.3
Sciatica	21,801	59.9	71.9	32	65.8
Cervical strain	4,156	61.1	50.2	50	84.1
Anxiety	19,857	61.8	73.9	32	66.7
Threatened abortion	6,986	69.0	72.3	38	74.7
Carpal tunnel syndrome	3,859	78.0	67.2	57	93.4
Depression	6,437	82.4	87.1	49	74.3
Medical meniscus injury	2,801	90.2	75.6	64	94.7
Adjustment disorders	3,079	107.0	86.6	82	92.3
Total	370,076	7.0	48.9	7	30.7

Note. SD = standard deviation.

Table 2. Duration of temporary disability by the remaining relevant variables.

Predictor	Category	N	Mean	SD	Median	Over 15 days (%)
Payment	Indirect	335,593	19.4	38.9	5	24.6
	Direct	34,483	91.7	78.9	67	89.9
Regime	Employed	326,127	22.6	45.0	6	26.6
	Self-employed	23,814	66.1	69.7	43	76.0
Civil status	Married	73,227	42.6	62.3	14	48.7
	Divorced	7,283	49.1	67.6	18	52.4
	Couple not married	3,062	44.4	64.4	14	48.5
	Single	37,620	34.9	56.2	9	40.9
	Widow	1,156	56.9	74.6	26	59.0
	Not available	247,728	18.9	39.6	5	22.8
Children	Yes	55,647	48.6	66.6	19	53.9
	No	161,703	26.3	48.8	7	31.2
	Not available	152,726	17.8	37.56	5	21.6
Codiagnosis	Yes	4,647	99.8	77.2	75	98.1
	No	365,429	25.2	47.7	6	29.8
Occupation	Armed forces occupations	122	55.6	69.1	32	63.9
	Manager	2,498	62.9	71.4	39	72.5
	Professional	9,378	32.6	50.8	11	43.4
	Technician & associated professional	17,741	34.8	54.9	11	43.4
	Administrative secretary & related associate professional	9,774	28.4	49.2	8	35.4
	Service and sales worker	21,501	36.0	55.6	11	43.4
	Skilled agricultural, forestry & fishery	1,255	36.4	58.0	9	41.8
	Craft and related trades worker	16,282	31.9	55.0	8	37.7
	Plant & machine operator, assembler	12,055	29.0	52.7	8	34.8
	Elementary occupation	17,453	33.0	56.6	8	37.6
	Not available	262,017	23.1	45.8	6	26.6
Age	<32	87,922	17.7	36.9	5	21.1
	32–37	92,793	24.3	44.2	6	30.4
	38–42	62,905	25.7	47.8	7	31.0
	43–51	68,894	28.8	52.6	8	32.9
	>51	56,550	39.8	64.0	11	43.1
Zone	Zone 1	32,770	38.9	58.6	13	46.2
	Zone 2	131,094	33.4	55.5	9	38.1
	Zone 3	73,513	21.3	42.9	5	25.5
	Zone 4	132,699	18.5	40.0	5	22.4
Gender	Male	183,683	22.7	46.3	5	26.1
	Female	185,103	29.5	51.1	8	35.2

Note. SD = standard deviation.

Statistical analysis

The statistical software SPSS 19.0 (Chicago, IL, USA) and R 3.2.2 (R Foundation for Statistical Software Platform) were used to analyze the data. Due to the large size of the sample, the main decisions throughout the statistical process depended on the performance of the different models according to a random subsample joined by 30% of data. In other words, 258,711 entries were used as a training set to fit the different models, and 111,365 entries were used as a validation set.

For quantitative prediction, a WLS multi-factor ANOVA¹⁰ and a gamma generalized linear model¹¹ (with a logarithmic linkage function) were applied to fit the variable log-duration and duration, respectively, considering the categorical predictors and the training sample. The logarithmic transformation was considered in both models due to the strong skewness of the variable duration. Skewness also made us consider median instead of average in some parts of study.

Sources of variability in the ANOVA model (we are considering here variables as well as interactions between them) were carefully chosen according to corrected R² so that we could explain the maximum percentage of variance with the minimum degrees of freedom. In addition, the influence of each source of variability was measured by eta-squared partial coefficient (η²), so that those with η² under 0.001 were not considered in any case (see Figure 2). In fact, the sources of variability removed from ANOVA model were no longer considered in the rest of the study. The remainder of the variables are described in Table 2, where they appear sorted by their respective η² (when the interaction of a variable with diagnosis was big enough, η² corresponding to interaction was added up to the one corresponding to the variable itself in Figure 2).

Figure 2. Importance of main variables in categorical and numerical models.

For categorical prediction, we carried out a logistic regression¹² and a classification tree,¹³ in which we considered CRT method, the Gini index as measure of error, a prune coefficient of 0.5, and a maximum depth of 10 levels in the growth of the tree. The influence of each variable on categorical prediction was measured by its capability to reduce the Gini index. We also tried to improve the classification tree by carrying out AdaBoost,¹⁴ the most popular boosting algorithm,¹³ implemented in R-package Adabag.¹⁵ We considered that the other possible improvement of classification tree, Random Forest,¹⁶ was not appropriate in this case because of the particular conditions of our predictors.¹³ Finally, a neural network algorithm¹⁷ was applied.

To illustrate the joint correlation between all the categorical variables involved in categorical prediction, we applied a multiple correspondence analysis¹⁸ (MCA), which provided the Figure 3. There, each category is identified by a symbol specific to the involved variable (circles for both types of categorical duration, stars for diagnostics, pentagon for occupations, etc.). Proximity between categories may be roughly understood as direct statistical association if they belong to different variables and as a similar statistical behavior if they belong to the same variable.

Results

Quantitative prediction

Duration of TD had a strong skewness and median 7. In other words, temporary disabilities were usually short although many of them were long or very long. Our challenge was to explain the behavior of data by a suitable statistical model. ANOVA was applied to fit the variable log-duration from the categorical predictors considering just the training sample. A value of 0.635 was reached for corrected R² considering the 10 variables we explained in the preceding and some interactions between them (112 degrees of freedom in all). The behavior of the model residuals turned out to be appropriate. Taking into account the validation sample, ANOVA predicted 10.8 days less than the real duration in mean with a standard deviation (SD) of 43.8. However, its error was close to 0 in terms of the median. The main factor was clearly diagnosis, followed by payment. The interaction between diagnosis and some other variables (especially payment) turned out to be also relevant. In Figure 2, we can value the influence of each variable according to its respective η² coefficient (numerical model). A Gamma GLM was also applied to fit the variable duration considering the same sources of variability in the model above, but it performed quite similar to ANOVA.

Categorical prediction

To predict whether TD would last more than fifteen days, we tried first the logistic regression model, considering as predictors the sources of variability included in the ANOVA model and, as response, the categorized duration of TD. We obtained a Nagelkerke R² coefficient 0.660 with the training sample, and 17.7% of misclassification error rate with the validation sample (12.0% for brief TD and 19.2% for long TD). Second, we applied the method that is, from our point of view, the most natural and convenient to deal with this problem: the decision tree applied to the variables selected in ANOVA as predictors. It led to 14.2% of misclassification error rate with the validation sample (9.5% for brief TD and 24.9% for long TD). The importance of the variables according to their ability to reduce the Gini index is also shown in Figure 2 (categorical model) together with η² (Gini index expresses the importance of a variable at predicting whether a TD will last more than 15 days, while η² expresses the capability of a variable at predicting an exact duration). Adaboost algorithm did not improve the fitting of the classification tree since the misclassification error rate was almost the same as that in the simple tree, and the same thing happened with the neural network algorithm.

The main variables

There exists a strong consistency between the different models considered in the preceding since all of them led to similar conclusions. Due to its better interpretability, we would use the simple classification tree for categorical predictions and would take into account the ANOVA model for numerical predictions.

Diagnosis was the most important variable in both cases, as described in Table I. The other relevant variables were, sorted according their influence in the classification tree, payment, regime, civil status, children, codiagnosis, occupation, age, geographical zone, and gender. They are described in Table 2. The purpose of Figure 4 is to illustrate the data from a numeric point of view as simply as possible. We can check there the medians of the duration depending on the two main factors in the numeric model: diagnosis and payment. We should clear up that regime showed also an important influence on duration (see Figure 2) although it is strongly correlated with payment (namely, we checked that direct payment corresponded to 82.1% of the self-employed and just 3.8% of employees). To summarize data from a categorical point of view, all the relevant categories are represented by MCA in Figure 3.

Figure 3. Categorical relationships between the relevant variables.

Figure 4. Median of duration by diagnosis and payment (direct or indirect).

Comment

While categorical prediction (classification tree) achieved an acceptable fitting, we think that quantitative prediction turned out to be quite ambitious in this case, since both ANOVA and gamma models are not accurate enough to make reliable predictions. We also think that, in general and for very big samples like ours, the lack of accuracy cannot be explained by the statistical model chosen but by the lack of relevant information (variables) in the data. Let us see the variables we have (although we prefer to talk about factors).

Diagnosis

Diagnosis was by far the most important factor, as we can see in Figure 2 and Table 1. Figure 4 shows which diagnoses are in general associated with longer TD, namely, those at the right part of the straight line. In fact, the first step in the classification tree roughly lies in knowing in which side of Figure 3 the diagnoses are located. Despite the large variability interdiagnosis, we also found a huge variability intradiagnosis that could be partially explained by the other factors.

Payment regime

Payment is the primary other factor. It is quite surprising how duration increased considerably with direct payment, especially for some specific diagnoses (see Figure 4) such as depression (one month for indirect and four months for direct in median), anxiety, threatened abortion, low back pain, headache, and migraine. Nevertheless, we found a smaller difference in other diagnoses such as inguinal hernia or carpal tunnel syndrome. We should study payment together with regime since they were bound from a statistical point of view (see Figure 3). In addition, the huge standard deviation of duration for direct payment suggests a heterogeneous behavior in this group. This fact has been pointed out in previous studies from smaller samples.^6,9 We will try to explain this subject better in the section devoted to unemployment and other irregularities.

Age and familial responsibilities

Civil status showed an important influence in the model, since medians went from 6 days for singles to 26 for widows/widowers. We can also see in Table 2 that duration of TD was greater for people with children. Moreover, it was clear that duration of TD was greater as age increased. We can see in Figure 3 that these three variables were correlated in the sense that youth was associated with single status and lack of children, and these circumstances led to brief TD. Therefore, this fact could be explained as the influence on duration of two correlated factors: a biological factor, namely age, and a sociological factor that we could name “familial responsibilities.”

Codiagnosis

As expected, the existence of a codiagnosis is clearly associated with longer TD as we can see in Table 2. However, its influence in the model is only moderate (see Figure 3) because there were very few entries with a codiagnosis. That is why there is not a specific figure for this variable.

Occupation

The influence of occupation in the model was smaller than expected. Although the longest median TD corresponds by far with managers and armed forces (see Table 2), we must point out that these two occupations presented also by far the largest percentage of direct payment (59% and 53%, respectively). Thus, we could regard the payment regime as the real source of this influence. We did not find a strong correlation between diagnosis and occupation, except in the case of inguinal hernia and skeletal-muscle pathologies, which are associated with craft, agriculture, forestry, fishery, and plant and machine operators. Some of these statements could be considered similar to conclusions founded in previous studies.⁹

Geographic location

Zone played a more important rule in the numerical model than in the categorical model (see Figure 3 and Table 2). We can also see in Table 2 that Zone 1, corresponding mainly to the west and south of Spain, had a median of 13 days of TD, while Zone 4, composed mainly of the most developed areas from a socioeconomic point of view in north and central Spain, had a median of just 5 days.

Gender

Gender turned out to be relevant in the model due to its moderate interaction with diagnosis and with age. For example, low back pain, which, according to Table 1, is the most common diagnosis, lasted a median of 11 days for men and 22 days for women. Something similar happened with sciatica.

Differences between genders were most observed under 42 years old. Since there was little interaction between gender and children in the ANOVA model, we cannot explain these differences in terms of familial responsibility. Therefore, according to the data, we are inclined to think that young women are more sensitive than young men to lumbar problems.

Some other variables not considered

When we say that the 10 variables studied were the most relevant, we are not saying that the predictors removed at the beginning of the study had no correlation with duration. We just mean that they explained nothing apart from what the 10 main variables explain. In the case of smoking status and drinking status, we found not even a simple correlation with duration. Nevertheless, this conclusion is not reliable since unfortunately these two variables were evaluated only for a small portion of special patients (3.2% and 2.7%, respectively). We feel these variables would have been able to improve the model if they had been evaluated better. The same problem (missing data) to a lesser degree happened with some other variables, such as children, civil status, and occupation (see Table 2).

Unemployment and other irregularities

What do payment and employment regime really depend on? Although these strongly correlated variables were shortly explained in previously, we add 3 remarks about them. First, it must be taken into account that 90.7% of workers were employees, and only 3.8% of them received direct payment (34,483 patients in all). Employees were paid indirectly through their company until the end of TD unless they were fired or their labor contract ended during TD. In such case, payment was assumed directly by Ibermutuamur from that moment.

From this point of view, long disabilities could be understood as a refuge against unemployment. So, while many employees probably tried to hide their sickness or to shorten their TD for fear of being fired,⁴ others probably tried to prolong TD because they were about to get fired or had been already fired (direct payment). This could explain the extremely large skewness of TD duration that has been pointed out in previous studies.⁸ We are even inclined to think about a vicious circle of disability-dismissal: Many employees were probably fired because of disability, and they prolonged disability because of it.

Second, 61.3% of direct payment corresponded to self-employed (managers). According to our experience, it is possible that many of these workers kept charging Ibermutuamur while their companies were really broken, so that TD could be considered to be a refuge against unemployment as well. It is also very likely that other managers kept running their small business in an irregular way while receiving a direct payment from Ibermutuamur for a long time, as has been described in recent studies.⁸

Third, as we can see in Figure 4, these kinds of unexpected long TD are associated with typical diagnoses that are hard to check. We can see that, in these cases, the median of direct payment is about 4 times (or more) the median of indirect payment. Thus, there exists a not strictly medical component in the diagnosis. All these facts help us to understand the huge importance of nonmedical factors in TD.

Limitations of the study

Although temporal disability is a global subject, we found it difficult to compare this study with others in different countries because of the sociological differences and obvious variations between legislations. In addition, we have taken into account a few recent studies in Spain dealing with this issue that were not published in English, so it could be hard to compare our conclusions with the statements of previous studies. Nevertheless, the most conflictive point in the study (unemployment) has been explained according to the professional experience of the authors.

Moreover, despite the large sample size, several circumstances could spoil some conclusions. First, we did not control data gathering. As a result, some important variables, such as smoking and drinking status, were not properly registered. Furthermore, medical diagnosis should have been less specific in order to analyze the whole data. Second, since Ibermutuamur is a widespread company in Spain, we think that our sample is representative, but we are not able to assure this statement. Third, the sample was not stable over the years.

Conclusions

We have obtained an algorithm to distinguish between brief (15 days at most) and long (over 15 days) TD, depending on some biological and medical factors— diagnosis, codiagnosis, age, and gender—and on some socioeconomic factors—unemployment, familial responsibility, geographic location, and others. The influence of the second kind of factor is large and difficult to control. We think this is the reason numerical predictions are not accurate enough. Nevertheless, we have studied them as deeply as possible.

Funding

This work was supported by the Junta de Extremadura (Autonomous Government of Extremadura, Spain) under the project GR15013.