Back disorders are a major health problem in many occupations. Physical loads at work, such as manual materials handling, frequent bending and twisting, lifting, and forceful movements, have been identified as possible risk factors for low back pain (LBP). Furthermore, psychosocial work factors may have a contribution as well; the results of several studies suggest that Karasek’s model—that is, high work demands in combination with low job control and decision latitude, affect the risk of LBP.^1–3 Many of those studies, however, did not include measures of both physical and psychosocial work factors as well as other potential confounders. Although we know very little about dose-response associations between occupational exposures and LBP risk, some authors have suggested that the relation between physical work load and LBP is U-shaped rather than monotonic.^4,5 Apart from the need to assess quantitatively the broad array of occupational exposures, we must also consider the full spectrum of LBP, including minor episodes of short duration as well as severe and disabling episodes that persist for long periods.^6–8

In order to estimate dose-response relations with sufficient detail, quantitative information on exposures to possible risk factors is necessary. In the past two decades, observational methods have been developed to facilitate a quantitative assessment of physical load at the workplace.^9–13 The large measurement efforts involved with these methods require statistical techniques that make complete use of all data. Furthermore, these statistical techniques should allow for possible non-linearities between exposure and risk, for example, using spline analysis.^14,15

In many occupations, exposure to physical load is not limited to one physical factor but encompasses many factors such as trunk flexion, rotation, lifting, and carrying. Often these physical factors are strongly inter-correlated. In epidemiological studies, conventional approaches for analysing such correlated exposures (for example, deleting certain variables from the model or adjusting the alpha level) may invoke problems of multiple comparisons.¹⁶ Several authors have shown that hierarchical regression models provide an alternative to these conventional approaches and can have important advantages in the adjustment for covariates that are strongly interrelated.^16–19

In the study reported here, dose-response relations of physical load at work, psychosocial work factors, and other factors with the incidence of LBP and LBP with disability are investigated simultaneously.

METHODS

The study population consisted of workers from seven nursing homes and homes for the elderly in the Netherlands. Workers employed for more than 10 hours per week were eligible to participate in the study. Baseline measurements performed between March 1998 and March 1999 included the assessment of physical load for each occupation at the workplace and a self administered questionnaire. Follow up measurements were performed one year after the baseline measurements.

Of the 1208 subjects invited to participate in the study, 769 (64%) agreed to participate. After one year, 523 (68%) of the 769 subjects were observed again. These subjects worked in nine professions: 85 nurses, 197 care givers, 26 kitchen workers, 23 housekeepers, 9 transportation and maintenance workers, 8 laundry workers, 24 physical therapists, 108 office workers, and 43 miscellaneous workers. For the analyses reported in this paper, we selected those workers who reported at baseline no LBP in the previous 12 months.

Data collection

Quantitative assessment of physical load

In a group of 212 workers randomly sampled at baseline, observations were made at the workplace on physical load during work for each of the nine occupations. The proportion of observed subjects per group ranged from 16% for the nurses to 64% for the transportation and maintenance workers, thereby oversampling those occupations with fewer subjects. The primary aim was to have at least 10 workers in each occupation in order to arrive at a reasonable estimate of the average exposure of each group.⁹ In the three largest occupational groups, at least 10 workers were sampled in four different nursing homes in order to increase generalisability across homes. An observational multimoment method was used to describe three measures of physical load: trunk flexion between 20 degrees and 45 degrees, trunk flexion over 45 degrees, and lifting or carrying loads over 10 kg, all expressed in percentages of work time. Observations were made on selected workers every 20 seconds during four periods of 30 minutes each in one working day; thus, we collected 360 observations per worker. For each occupation, the average exposure to each type of physical load was calculated as the mean percentage of time devoted to that activity. An underlying assumption of this strategy is that the average physical load of the observed workers is equal to the average physical load of the total occupational group. Thus, the arithmetic mean for the group was used as a proxy for exposure in all subjects, both observed and unobserved, in that occupation.^9,11–13 To estimate an individual’s cumulative exposure to physical load, the occupation specific exposure expressed in percentage of work time was multiplied by self reports of the number of hours worked per week.

Questionnaire survey

A questionnaire was used to collect personal data on age, job history, including the number of years employed at that facility, and LBP in the past 12 months.

Information was collected on three psychosocial work factors included in Karasek’s model:²⁰ decision authority, skill discretion, and work demands. The 11 questions on decision authority reflected aspects such as influence on the planning of tasks, influence on the pace at work, brief pauses when needed, and decisions on time spent on given tasks (Cronbach’s alpha = 0.90); the six questions on skill discretion reflected required skills, task variety, learning new things, and amount of repetitive work (Cronbach’s alpha = 0.80); and the 11 questions on work demands reflected working fast, working hard, excessive work, insufficient time to complete the work, and conflicting demands (Cronbach’s alpha = 0.88). The response to each question was scored on a four point ordinal scale. The total score for each factor was computed by summing the item scores for all questions related to that factor and expressing the sum on a 0–100 scale. A score of 100 was defined as least desirable for that factor—that is, low decision authority, low skill discretion, and high work demands.^8,21

Information on the occurrence and nature of LBP was measured at baseline and at follow up with questions derived from the Nordic questionnaire for the analysis of musculoskeletal symptoms.²² This questionnaire has been shown to be valid for collecting information on the nature, duration, and frequency of symptoms. In this study, LBP was defined as any episode of pain that lasted for at least a few hours in the previous 12 months. An incident case of LBP was defined as the report of LBP at follow up among subjects who reported no LBP at baseline. LBP with disability (LBP/D) in the previous 12 months was defined as a report of LBP for which the Von Korff et al disability score was greater than 50, indicating “high disability”.²³ This second outcome event is a subset of the first. An incident case of LBP/D was defined as the report of LBP/D at follow up among subjects who did not report LBP/D at baseline (but who might have reported any LBP at baseline).

Data analysis

Log-linear models with quadratic splines were used to assess dose-response relations for eight possible predictors of LBP and LBP/D incidence: trunk flexion between 20 and 45 degrees, trunk flexion over 45 degrees, lifting and carrying loads over 10 kg, decision authority, skill discretion, psychosocial work demands, years of employment in the facility, and age at baseline. The choice for these predictors and their particular cut-off values were based on a literature review and a cross sectional analysis^1,13 These variables were entered in the multivariate model simultaneously. In order to capture non-linearities all continuous variables were categorised into tertiles of the observed distributions to allow for a quadratic spline function.^14,15 On the basis of likelihood ratio statistics, it appeared that a full quadratic spline (for all categories) did not fit the data better than did a restricted quadratic spline in which a linear spline was used for the lowest category. Since this restricted spline has one linear and two quadratic terms per factor, the eight predictors resulted in 3×8 = 24 parameters and an intercept that had to be estimated. The model was fitted using Proc Genmod (with a log link function) available with SAS statistical software.²⁴

This conventional (one-stage) method may result in large standard errors when the number of parameters in the model is large or these parameters are highly correlated, due to instability of the maximum likelihood estimates.¹⁶ As shown by Greenland and colleagues,^16–18 the maximum likelihood estimates can be improved by fitting a second stage linear model in which a priori similarities among certain parameters are taken into consideration. The variances of the physical factors as well as the psychosocial factors were constrained and assumed to have a similar magnitude, thereby introducing more stable and accurate estimates. In this hierarchical regression approach, estimated parameters from the first stage are treated as outcomes in the second stage. For these analyses, the eight predictors were grouped into three relatively homogeneous domains: the three physical factors, the three psychosocial factors, and the two time related factors (age and years of employment). In the second stage, we assumed that the effect estimates for the factors within a domain are similar. One second stage covariate per domain is not sufficient to describe similarities among the predictors because we have both linear and quadratic terms for each predictor and because we cannot assume that these two terms are similar. Therefore, two second stage covariates per domain were used, which yielded six second stage covariates in the model. The actual values of the second stage covariates (design matrix with elements z_pc, with p = 1, ..., 24 and c = 1, ..., 6) were calculated as follows: z_pc = 1/w for the linear parameters and z_pc = 1/w² for the quadratic parameters, where w is the length (in units of exposure) of the category corresponding to the particular parameter. With these values assigned to the covariates, the second stage model described the prior curve of the spline for each of the eight predictors with the relative risk at a given centile of exposure equal for all factors in a domain.

A two step procedure, described Greenland and Witte,^16,18,25 was used to do the hierarchical regression analyses. The second stage variances τ² could be estimated from the data with empirical Bayes methods; however, the large number of parameters relative to the sample size led us to pre-specify the τ² to upper and lower bounds for the parameters of the splines (that is, a semi-Bayes approach).^16–18 We set τ_p² to 0.5, which corresponds to a prior certainty of 95% that the relative risk (RR) per unit exposure lies in a 16-fold range, that the RR per unit exposure, reflecting departure from linearity comparing the first and second categories, lies in a 16-fold range, and that the difference in the quadratic term between the second and third categories also lies in a 16-fold range.

For both the conventional spline models and the hierarchical spline models, estimated relative risks and corresponding 95% confidence intervals (CI) are presented at the 25th, 50th, 75th, and 90th centiles of the exposure distribution, relative to the 10th centile. Effect estimates are presented for both the conventional and hierarchical models to compare the point estimates, precision, and stability of the estimates.

RESULTS

The estimated prevalence of LBP at baseline among workers in each elderly facility was not related to the response rate in that facility. Of the 769 subjects observed at baseline, 246 (32%) were not available to provide outcome information at follow up. Of those lost to follow up, 39 subjects changed jobs to employers not participating in the study. The estimated prevalences of LBP and low back pain with disability (LBP/D) at baseline among subjects lost to follow up were similar to the prevalences among subjects available for follow up. Among those subjects free of LBP at baseline, no differences in the physical factors were found between workers available for follow up and those lost to follow up. Those lost to follow up, however, were younger and had fewer years of employment at the facility. In addition, among those subjects free of LBP/D at baseline, those who were lost to follow up reported less exposure to lifting and carrying loads over 10 kg and less skill discretion than did those not lost to follow up. No differences were found for decision authority and work demands between these groups.

Descriptive findings

In the baseline cohort of 523 workers, the estimated prevalence of LBP was 57.9%, and the estimated prevalence of LBP/D was 8.4%. The prevalence of LBP one year later was 54.5%, and the prevalence of LBP/D was 11.6%. The estimated one year cumulative incidence (risk) of LBP was 58/220 = 26.4%, and the estimated one year cumulative incidence of LBP/D was 42/479 = 8.8%. Table 1 presents the estimated cumulative incidences of LBP and LBP/D for each occupation. Although the estimated cumulative incidence of LBP and LBP/D did not vary appreciably across occupations (at least for those occupations with more than 10 subjects at risk), the incidence of LBP was relatively low for nurses, and the incidence of LBP/D was relatively low for office workers.

Table 2 presents the means and standard deviations of the eight exposures at baseline for all 523 workers. These results were similar for those subsets of workers free of LBP at baseline (n = 220) and free of LBP/D at baseline (n = 479). The greatest exposure to trunk flexion between 20 and 45 degrees was found among the therapists (5.7 h/wk), followed by the nurses (5.3 h/wk) and office workers (4.7 h/wk). Nurses experienced greater exposure to trunk flexion between 20 and 45 degrees than did caregivers (4.2 h/wk). The greatest exposure to trunk flexion over 45 degrees was found among transportation and maintenance workers (2.6 h/wk), followed by housekeepers (1.7 h/wk) and caregivers (1.3 h/wk). Caregivers experienced greater exposure to this type of trunk flexion than did nurses (1.0 h/wk). Furthermore, the greatest exposure to lifting and carrying loads was found among nurses (38 min/wk). No other appreciable differences among occupations were observed for the other exposures.

Dose-response findings

Table 3 presents adjusted dose-response associations for the estimated effects of eight exposures on the one year risk of LBP. The relative risks (and 95% CIs), comparing selected values of each exposure, were estimated using a conventional (one stage) model and a hierarchical (two stage) model. Table 4 presents similar results for LBP/D. Note that, in general, the 95% confidence intervals are narrower with hierarchical regression than with conventional modelling.

Physical work load

Although the estimated risk of LBP is monotonically associated with trunk flexion over 45 degrees, using hierarchical regression, the association is weak and the relative risks for high exposures are imprecisely estimated. The association with this exposure is much stronger, however, for LBP/D. The relative risk estimated from the hierarchical model was 3.18 (95% CI 1.13 to 9.00) for 1 hour and 45 minutes of bending per week (90th centile) relative to 30 minutes per week. There was little association between trunk flexion between 20 and 45 degrees and either outcome. For lifting and carrying loads over 10 kg, there was little association with LBP/D using either method, but there was an unexpected inverse association with LBP.

The 95% confidence intervals were narrower for effect estimates derived from hierarchical models than for those effect estimates derived from the conventional models. This difference was largest for the estimated effect of trunk flexion over 45 degrees on the incidence of LBP/D. For the effects of trunk flexion between 20 and 45 degrees on LBP and LBP/D and for the effect of trunk flexion over 45 degrees on LBP/D, the point estimates were closer to unity with the hierarchical model. The shape of the dose-response curve was not altered appreciably and remained linear for the hierarchical model. In contrast, for the effect of trunk flexion over 45 degrees on LBP and for the effect of lifting and carrying loads on LBP and LBP/D, the estimated dose-response relations changed from non-monotonic or inconsistent with the conventional model to more linear with the hierarchical model.

Psychosocial work factors

We observed adjusted associations between LBP incidence and low decision authority and high work demands, but little association with skill discretion. In all cases, however, the confidence intervals were wide. No consistent dose-response associations were observed between any of these psychosocial variables and LBP/D. For both outcomes, effect estimates were similar with conventional and hierarchical models.

Time related factors

We observed an inverse monotonic association between age and the incidence of LBP; the estimated relative risk derived from the hierarchical model, comparing 55 year olds with 25 year olds was 0.33 (95% CI 0.12 to 0.88). In contrast, a weak but positive association was observed between age and LBP/D; the corresponding estimate of the relative risk was 1.68 (95% CI 0.57 to 4.90). Years of employment in the facility was inversely associated with LBP/D but relatively unassociated with LBP. Similar findings were obtained for these variables with conventional and hierarchical models.

DISCUSSION

Few differences in the incidence of any LBP were found between occupational groups. Although a few of the work related factors were consistently associated with the incidence of LBP, these associations were not precisely estimated and may have been chance findings. Younger workers were at higher risk of any LBP than were older workers. Findings differed somewhat for the other outcome—low back pain with disability (LBP/D). In particular, younger workers were at lower risk than were older workers, and trunk flexion over 45 degrees was strongly and consistently associated with LBP/D. It appears, therefore, that trunk flexion at work affects the risk of more severe or disabling low back pain but not the risk of less severe or transient low back pain. Since the objective of this study was to investigate the incidence of low back pain among workers who had reported at baseline no low back pain during the previous year, the results cannot be generalised to predicting the recurrence of low back pain among chronic cases.

Methodological issues

The participation rate among the nine professions and seven nursing homes varied considerably, but was not related to the prevalence or severity of LBP at baseline. Hence, selective participation does not seem to pose a problem in this study, although the lack of information on individual characteristics of non-participants prohibits a further underpinning of this important assumption. Subjects lost to follow up were younger, had fewer years of employment, reported less exposure to lifting and carrying loads, and reported less skill discretion than did subjects not lost to follow up. Thus, the presence of bias due to selective loss to follow up cannot be ruled out. Loss to follow up is determined by non-response and job change (to an employer other than those participating in the study) during follow up. Of those lost to follow up, 39 subjects changed jobs. Furthermore, we know that there is a high turnover rate in the source population during the first years of employment and that changing jobs was associated with age and years of employment among our subjects. These factors might explain the difference in age and years of employment between subjects lost and not lost to follow up.

An observational method was used to measure physical load at the workplace. The advantage of this approach is the ability to collect detailed quantitative information. In earlier analyses of these data, we found that this approach was subject to a certain degree of random measurement error, especially for lifting and carrying loads, where the within subject variance captured up to 80% of the total variance. In general, random or non-differential measurement error results in underestimation of effect.²⁶ In our study, the direct measurement of physical load was not made at the individual level but at the occupational group level. To measure an individual worker’s cumulative exposure to physical load, the average percentage of work time devoted to that activity in the worker’s occupation was multiplied by that individual’s self reported number of hours worked per week. It is difficult to predict whether, and to what extent, this approach resulted in biased estimates of physical load effects in our study and in what direction the bias is likely to have occurred.²⁷

In many previous epidemiological studies of LBP, the study population has been limited to an occupational group with a more or less fixed ratio of high levels of exposure to different physical factors and a reference group with low levels of exposure. Furthermore, when estimating the effects of physical factors that were strongly interrelated, those researchers typically did not fully adjust each effect estimate for all other factors. Those results, therefore, might have been confounded. In this study, correlation coefficients among physical factors ranged from 0.29 to 0.59. We estimated the effect of each occupational exposure by controlling for all other exposures—first, by including several occupational groups with different physical exposure distributions; and second, by using a hierarchical regression approach that enhanced effect estimation for multiple interrelated exposures.^16–19

Issues in the hierarchical modelling of the dose-response splines

When results from the conventional and hierarchical models were compared, clear differences were found only for the physical factors, especially for the effect of trunk flexion over 45 degrees on LBP/D. Estimated confidence intervals tended to be narrower with hierarchical modelling than with conventional modelling. These findings are in accordance with the results of earlier research on log-linear hierarchical models.^16–18 In general, with hierarchical regression, large and unstable effect estimates are pulled (regressed) towards a common mean. With a linear model, the common mean or prior distribution describes a monotonic linear function; in practice, this often means that extremely unstable effect estimates are reduced towards unity. In contrast, with a spline model, the prior function is non-monotonic (when the values of the second stage design matrix are not zero). Essentially, this prior function is the mean of all dose-response curves found with the conventional model in one covariate domain. Unstable effect estimates are regressed towards this mean spline. The implication, in contrast to linear models, is that hierarchical regression affects not only the point estimates of effect, but also the shape of the dose-response relations.

In order to minimise bias in the estimates caused by regression to an inappropriate common mean or prior function, we must assume that parameters are exchangeable, and quantitative exposures must be scaled so that regression estimates for one unit of exposure can be compared.¹⁶ Therefore, in this study, we grouped the variables into three domains (physical factors, psychosocial factors, and time related factors), and we identified exchangeable linear and quadratic parameters as described by positive second stage covariates. A comparison of regression estimates was possible since the same unit of exposure was used for all variables within a domain.

When τ² was preset to 0.25 (instead of 0.5), which implies larger regression towards the prior curve,^16–18 results were different only for estimates of less stable physical factors. Estimated relative risks for the effects of trunk flexion over 45 degrees were about 20% smaller and 95% confidence intervals were narrower when τ² was set to 0.25, whereas estimated relative risks for the effects of lifting and carrying loads were somewhat larger. When τ² was set to 1, estimated relative risks obtained from the hierarchical model were closer to the maximum likelihood estimates and the 95% confidence intervals were wider than when τ² was set to 0.5.

Comparison with other studies and evaluation of results

Several reviewers have presented clear evidence that physical work load contributes to the occurrence of LBP.^1,2 Since our study population included many nurses and caregivers, it is informative to compare our results with the results of other LBP studies of nurses. Nursing work is physically strenuous due to patient handling, and a major part of the work involves trunk flexion and lifting. The results of our study do not confirm the association between lifting and LBP that has been reported in other studies.^28,29 One explanation for our negative finding might be random error in measuring this exposure (as indicated above). Another explanation might be selection factors operating before baseline, for example, if workers previously exposed to lifting and carrying heavy loads and at high risk of LBP were more likely to leave their jobs or reduce their strenuous activities at work.

Other epidemiological studies of LBP that incorporated observational techniques to measure physical load at work were conducted by Punnett and colleagues³⁰ and Hoogendoorn and colleagues.³¹ Punnet et al used the same cut-off points for trunk flexion that we used in our study: 20–45 degrees, and over 45 degrees. They also found a positive association between trunk flexion and the risk of LBP, but their effect estimates were larger than ours. In their prospective cohort study of low-back pain, Hoogendoorn et al found similar estimates for the effect of trunk flexion that we report in this paper. They also found a positive association between lifting loads of at least 25 kg more than 15 times per work day and the risk of LBP, but the association was weaker and inverse for lighter loads and lower frequencies of lifting. This latter finding is similar to the results of our study.

Although there is evidence that psychosocial factors at work affect the risk of LBP, the specific factors responsible for such effects are not well understood.^3,28 Many studies of psychosocial factors were cross sectional or did not involve adequate adjustment for physical factors and other potential confounders. In both the case-control study of nursing personnel by Josephson and colleagues³² and the cohort study of different worker groups by Hoogendoorn and colleagues,²¹ the investigators found weak associations between psychosocial work factors and low back pain, controlling for physical factors and other potential confounders. The evidence from our study does little to refute or confirm the results of those studies.

Results from several previous studies suggest that the occurrence of LBP increases with age.^1,33 We also found a positive association with LBP/D but an inverse association with LBP. In contrast, we found an inverse association between years of employment and LBP/D but no association between years of employment and LBP. These seemingly contradictory findings might be due to selection factors that operated before the start of follow up. It is possible that the observed risk of LBP was lower in older workers because many workers of that age had developed LBP before follow up and thus were excluded from this incidence study. Indeed, the estimated prevalence of LBP at baseline in the original worker population was about 70% in nurses and caregivers. It is also possible that the observed risk of LBP/D was lower in workers with more years of employment because workers with more years of employment and free of LBP/D at baseline represent a healthy (low risk) group of “survivors” who had not developed LBP/D before follow up.

There is a growing body of evidence suggesting that different types of back pain are affected by different factors.^6–8 In this study, we found that while LBP/D was more strongly associated with exposure to physical factors at work, LBP was more strongly associated with psychosocial factors. Although these associations may reflect true differences in effect, another possible explanation is differential reporting (that is, misclassification) of low back pain. It is possible, for example, that workers with less decision authority and more work demands might be more likely than workers without these stressors to report minor episodes of low back pain. Another possible explanation might be the difference in study populations for predicting our two outcomes: the prediction of LBP was restricted to only those 220 workers who were free of LBP at baseline, while the prediction of LBP/D was conducted in 479 workers who were free of LBP/D at baseline (56% of whom had LBP). Thus, the effect of occupational exposures to physical factors on LBP/D might be greater for workers with existing (non-disabling) LBP than in workers without any LBP. Unfortunately, we could not test that hypothesis in our study because there were too few incident cases of LBP/D among workers without LBP at baseline.

This study has several methodological strengths including the prospective cohort design, the detailed quantitative measures of occupational exposures, the use of hierarchical regression to adjust for confounders, and the use of spline techniques to assess dose-response relations. We conclude that occupational exposure to trunk flexion over 45 degrees appears to be a risk factor for low back pain with disability among persons employed in nursing homes or homes for the elderly in the Netherlands.