Revealing hidden curvilinear relations between work engagement and its predictors: Demonstrating the added value of generalized additive model (GAM)

Previous studies measuring different aspects of the quality of life have, as a rule, presumed linear relationships between a dependent variable and its predictors. This article utilizes non-parametric statistical methodology to explore curvilinear relations between work engagement and its main predictors: job demands, job control and social support. Firstly, the study examines what additional information nonlinear modeling can reveal regarding the relationship between work engagement and the three predictors in question. Secondly, the article compares the explanatory power of non-linear and linear modeling with regard to work engagement. The Generalized Additive Model (GAM), that makes possible non-linear modeling, is compared with the widely used simply linear Generalized Linear Model (GML) procedure. Based on the survey data (N=7867) collected in eight European countries in 2007, the article presents the following main results. GAM clearly fitted the data better than GLM. All investigated job characteristics had curvilinear relationships with work engagement, although job demands and job control relationships were almost linear. Social support had a clear U-shaped curvilinear connection to work engagement. Interactions between the three job characteristics were also found. Interaction between job demands and social support was curvilinear in shape. Finally, GAM proved to be a more practical and efficient tool of analysis than GLM in situations where there are reasons to assume curvilinear relationships, complex interactions effects between predictors.


Introduction
Previous studies measuring different aspects of the quality of life, and more specifically quality of work, have almost, as a rule, presumed linear relationships between a dependent variable and its predictors (e.g. job demands, social support and job control). In this article we will make good use of non-parametric statistical methodology to study curvilinear relations between work engagement and its main predictors. Work engagement is a rather novel concept that has been applied to measure the well-being of working populations. In contrast to work stress and burnout studies, it adopts a positive view regarding an employee's well-being at work (e.g. Mauno, Kinnunen & Ruokolainen 2007;Bakker et al. 2008). However, time after time, work engagement studies have been based on the idea of linear relationships between independent and dependent variables and presented the relationships as regression models (e.g. Hakanen, Bakker & Demerouti 2005;Mauno et al. 2007), or as structural equation models (e.g. Hakanen, Bakker & Schaufeli 2005;Schaufeli & Bakker 2004). This article has two aims which deal with empirical and methodological questions. The first aim of the study is to examine what additional information non-linear modeling can bring out regarding the relationship between work engagement and its three main predictors: job demands, job control and social support. With the help of this aim we aspire to obtain a more detailed understanding of variations in work engagement among working populations. The second aim of the paper is to compare the explanatory power of non-linear and linear modeling with regard to work engagement. By doing this we aim to find out what non-linear models have to offer for employee well-being and particularly work engagement studies compared with the traditional approaches based on linear modeling. Both aims are premised on the comparison of two statistical tools. The results provided by the non-parametric Generalized Additive Models (GAM) are compared with those based on the General Linear Model (GML) procedure.
In terms of data, we use a structured survey (N=7867) collected from Finland, Sweden, the United Kingdom, the Netherlands, Germany, Portugal, Hungary and Bulgaria in 2007. The survey was targeted towards the service sector, and the collected data represents: the retail trade, finance, banking, telecommunications and public hospitals.
The structure of this paper is as follows. First, we will take a look at previous studies dealing with work engagement in order to map out its main predictors and to determine what have been the dominant statistical methods in this field of research. Second, we will present Karasek's wellknown model, which is applied in order to predict work well-being according to job characteristics. We will use Karasek's extended model of job demands, control and social support (Karasek 1979;Karasek & Theorell 1990). Karasek's model and its successors, such as Job Demands-Resources (JD-R) model (Demerouti, Bakker, Nachrainer & Schaufeli 2001;, are widely exploited in work engagement research. In contrast to many previous studies, which have presumed linear relations between Karasek's predictors, we will investigate how non-linear modeling works with respect to this model and work engagement. The curvilinear modeling of job characteristics is theoretically motivated by Warr's vitamin model (Warr 1987), which presumes curvilinear relationships between job characteristics and well-being. After the introduction to the previous research, we will specify our main research hypotheses, and we will describe the data, measures and statistical tools used in the study.

Previous literature 2.1 Work engagement and its predictors
In the field of psychology, work engagement is considered as the positive opposite of burnout and stress (Mauno et al. 2007;Bakker et al. 2008). Work engagement is defined as a positive, fulfilling work-related state of mind, which is manifested in three dimensions: absorption, vigor and dedication (Schaufeli, Salanova, Gonzalez-Romá & Bakker 2002). Absorption refers to a state of being where one is completely occupied by one's own work and has difficulties detaching oneself from it. Flow is a closely related to the concept of absorption with the difference that flow refers to a short-term "peak" and absorption to a more pervasive and persistent state of mind (Schaufeli & Bakker, 2004). Vigor implies a high level of energy and mental resilience while working, especially when confronting difficulties at work. Dedication is a dimension of work engagement that manifests itself in a person's genuine involvement in work, including feelings of pride, significance and inspiration (Schaufeli et al. 2002;Spreitzer, Lam & Fritz 2010).
In prior studies, it has been usual to search for predictors for work engagement from workrelated demands and resources, which is a strategy consistent with a variety of theoretical approaches (Halbesleben 2010). For instance, the job demands-resources model (JD-R) developed by Demerouti et al. (2001) assumes that the main antecedents of work engagement are job demands and especially job resources. In the JD-R model job demands refer to the physical, psychological, social, or organizational aspects of the job that require constant physical or psychological effort. Job resources are the physical, psychological, social, or organizational aspects of a job that are functional in achieving work goals, reducing job demands and their associated physiological and psychological costs, or stimulate personal growth, learning, and development .
Prior studies have consistently indicated a positive relationship between work engagement and job resources such as social support, performance feedback and particularly job autonomy (Halbesleben 2010;Schaufeli & Bakker 2004). For example, the JD-R model predicts that job demands have a negative effect on work well-being only if demands are high and constant.
Although studies have indicated that job resources predict work engagement better than job demands (e.g. Mauno et al. 2007), negative relations between work engagement and job demands, such as work-family conflict, family-work conflict and work overload, have been found. The JD-R model also assumes that job resources buffer the negative effects of job demands on work engagement and that a positive relationship between job resources and work engagement is stronger when job demands are high ).
Regarding socio-demographic variables, the literature indicates the following types of relationships to work engagement. In relation to gender, the studies conducted have produced inconsistent results, although the inconsistencies are likely to be related to the special characteristics of the data used. Various studies have found that women are slightly more engaged in their work than men (e.g. Mauno et al. 2007), but also the opposite connections have been found (Schaufeli & Bakker 2003, Schaufeli, Bakker & Salanova 2006. Studies also reveal that work engagement increases slightly with the age of the respondents, but this relationship has turned out to be relatively weak (Schaufeli & Bakker 2003;Schaufeli et al. 2006;Van den Broeck et al. 2008). Also marital status seems to be linked to job engagement, but once again the connection is weak (Hakanen 2009). Compared with unmarried people, married people and cohabitants, as well as widowed and divorced people are more engaged in work (cf. Mauno et al. 2007). There are also signs that people with children living at home tend to be more engaged in work than their counterparts (Mauno et al. 2007). Prior studies also indicate that one's educational level could be positively connected to some aspects of work engagement, especially absorption (Mauno et al. 2007).
Some work-related structural factors also act as predictors of work engagement. With regards to shift work, Mauno et al. (2007) report that working shifts is not associated with any dimension of work engagement -at least in the health care sector, though it can be assumed that the connection may be different depending on the economic sector (cf. Demerouti et al. 2001). Cuyper et al. (2010) have also pointed out that the type of employment contract one has is associated with work engagement; temporary workers seem to be more involved in their work than employees with a permanent contract.
Sector and country comparisons on work engagement have been much rarer at least up to the time of writing. Schaufeli and Bakker (2003) have found some country differences in work engagement, yet they have not paid attention to differences between occupational groups. Many studies have discovered occupational differences as well (Schaufeli et al. 2006;Hakanen 2009).
To sum up, previous studies on relationships between job characteristics (job demands and resources) and work engagement are numerous (e.g. Bakker, Demerouti & Euwema 2005;Mauno et al. 2007;Schaufeli & Bakker 2004) but most of them utilize linear methods such as hierarchical multiple regression analyses and structural equation models. The only exception in this respect is a recent study by Sawang, Brough and Barbour (2009) where the curvilinear relationship between job characteristics and work engagement was analysed by using quadratic terms in a hierarchical regression analysis. Sawang et al. discovered curvilinear relationships between supervisor support and work engagement as well as between the control of time and work engagement. Curvilinear relationships between supervisor support and work engagement took the form of an inverted U-shaped curve. The low and high levels of supervisor support were connected to a low level of work engagement while an average level of supervisor support was connected to a high level of work engagement. Unfortunately, the form of the curvilinear relationship regarding the effect of control of time on work engagement remained rather vague.
Although research on curvilinear relationships between work engagement and job demands, job control/autonomy or social support has been rare, curvilinear relationships regarding other measures of work related well-being have been studied by utilizing quadratic terms in a hierarchical regression approach. De Jonge, Reuvers, Houtman, Bongers and Kompier (2000) have found curvilinear relationships between: job demands and emotional exhaustion (J-shape); social support and emotional exhaustion (U-shape) and between social support and depression (U-shape). De Jonge and Schaufeli (1998) found four curvilinear relationships between: job demands and anxiety (U-shape); job autonomy and emotional exhaustion (inverted U-shape); social support and job satisfaction (inverted U-shape) and social support and emotional exhaustion (U-shape). Also, a study by Rydstedt, Ferrie and Head (2006) indicates small curvilinear relationships between: social support and job satisfaction (U-shape) and between job control and job satisfaction (U-shape).

Karasek's extended model: could interactions be curvilinear?
In this study we take advantage of Karasek's extended Demands-Control-Support-model to examine the consequences of paid work on an employee's well-being measured as work engagement. Karasek's model (1979;Karasek & Theorell 1990) is one of the best known psychosocial balance models, along with the effort-reward imbalance (ERI) model (Siegrist 1996), and they are founded on the presumption that an employee's well-being is dependent on finding a balance between the demands employees are subjected to and the resources they have at their disposal ).
In Karasek's model, job demands refer to stress factors which influence how employees cope with workload, unexpected tasks and conflicts. Job control (or job decision latitude) refers to employees' possibilities to influence their own work settings. Further, the model separates the two components of job control: skill discretion, which refers to possibilities to be creative, use skills and learn new things at work; and job decision authority, which refers to the possibility to choose the way work is done and take part in decisions which affect one's own work (Karasek 1979;Karasek & Theorell 1990). This basic model was later supplemented with the dimension of social support (Karasek & Theorell 1990). In the extended model social support is understood as social interaction with colleagues and supervisors and has a positive impact on well-being at work (Johnson & Hall 1988;Karasek & Theorell 1990).
In Karasek's model, as in many other models drawing on it, the interactions between job demands and resources (both job control/autonomy and social support) have played a pivotal role.
Karasek's original model includes two main hypotheses related to this. On one hand, the strain hypothesis proposes that the risk of stress increases when work demands increase but job control decreases. On the other hand, active learning hypotheses claims that motivation and learning at work rise when both demands and control increase. Karasek's extended model includes the idea that the risk of stress becomes higher when work demands are high but job control is low and there is a lack of social support (Karasek & Theorell 1990). There is empirical evidence showing that high job demands and low job control are important predictors of wellbeing and motivation. Support for the buffer hypothesis i.e. that control can moderate the negative effects of high demands on wellbeing is less consistent ( Karasek's model has been criticized for overlooking the curvilinearity of the relationships between the main variables of the model. Warr (1987) suggests that Karasek's results in fact imply the relationships are curvilinear, and proposes that his vitamin model of well-being would predict non-linear relationships between job characteristics and work-related well-being. The vitamin model presupposes that job characteristics have either a constant or an additional decrement type of effect on well-being. Social support has a constant effect on work related wellbeing which means that at first the increase in social support has a positive effect on work related well-being but at some point an additional amount of social support does not have either a positive or a negative but a constant effect on well-being at work. Job demands and job autonomy/control have an additional decrement type of effect on work-related well-being. Job demands and job autonomy first have a positive effect on well-being at work, but when the levels of the demands and autonomy increase, the effect first becomes constant but then turns negative. Fletcher and Jones (1993) have also tested the curvilinearity of Karasek's model by using quadratic terms entered in a hierarchical regression analysis. Except for job and life satisfaction for men, they did not find curvilinear relationships between the model and outcome variables. In this study we will apply new, non-parametric and non-linear, methodology to find out whether Karasek's model could be used to predict variations in work engagement.

Analytic Strategy and Research Hypotheses
To meet the aims of the article we will present a set of statistical models based on linear (GLM) and non-linear (GAM) methods side by side. The main reason for this is to demonstrate what is missing when the relationship between work engagements and its predicators are presumed to be linear, and to show how curvilinear modeling enriches our understanding of how work engagement spreads within our data.
The chosen analytic strategy will also help us to test our research hypotheses, which are based on our literature review on empirical findings and require comparisons between GLM and GAM procedures. The hypotheses set for the study are: H1. The main effects of job demands, job control, and social support on work engagement are linear, and GAM does not provide much additional information when compared with GML.
H2. An interaction effect between job demands and job control on work engagement cannot be found, regardless of the applied statistical method.
H3. An interaction effect between job demand and social support on work engagement cannot be found, regardless of the applied statistical method.

DATA
This article is based on survey data (N=7867) collected from Finland, Sweden, the United The questionnaires were primarily delivered to the respondents in electronic form. An electronic survey was used to facilitate responses and to minimize errors related to the manual input of data. When an online survey was not feasible (e.g. due to limited or no internet connection at a workplace), paper-based questionnaires were used. The only country in which the web version of the survey was not used at all was Bulgaria. In other countries, the web-based version of the questionnaire was widely used for data collection. In organizations with a high degree of Internet access the web-based survey was typically preferred, while the paper-based questionnaire was more used, for example, in the retail sector (van der Lippe et al., 2011, 59-60).

Work engagement
Work engagement is a six-item composite variable based on the Utrecht Work Engagement Scale by Schaufeli and Bakker (2003)

Predictors of work engagement
The measurement of job demands, job control and social support is based on the Job Content Questionnaire (JCQ) which has been widely used in previous studies and based on Karasek's Demand-Control-Support model (Karasek and Theorell 1990).
Job demands (pressures). To measure job demands we used a shortened version of the scale developed by Karasek and Theorell (1990). In our study the job demand composite variable includes the following five components: "Does your job require you to work fast?", "Does your job require you to work very hard?", "Do you feel that your job requires too much input from you?", "Do you have enough time to complete your job? (item reversed)", and "Does your job often make conflicting demands on you?" High scores on this scale indicate that an employee experiences high pressure in terms of time as well as physical and mental effort (range 1-4, M=2.54, SD=.54, =.733).
Job control. Job control is a composite variable comprised of eight statements, such as, "Are you free to decide what your job involves?" and "Does your job require you to invent your own tasks?" (range 1-4, M=2.22, SD=.50, =.787). This scale is also adopted from Karasek and Theorell (1990). The high scores on this variable indicate a higher degree of independence in organizing work, which may also be reflected by improved possibilities for combining work and family life.

Controlled variables
Several personal and work-related controlling variables were included in the analyses based on the review of the prior studies presented above. Control variables related to personal factors were gender, age (categorised: age under 35 years, age 35-49 years, age over 50 years), married/living together (yes/no), children at home (yes/no), education (basic level, secondary level, lower third level, third level), and country (Finland, Sweden, UK, Netherlands, Germany, Portugal, Hungary, Bulgaria). Work related controlling variables were industry (Retail, Telecom, Hospital, Bank), supervisory position (yes/no), shift work (yes/no) and type of contract (permanent/temporary).

Validation of measurements
Prior to testing the hypothesis a confirmatory factor analysis using Mplus 6.12 (Muthén & Muthén, 1998-2012 was performed to validate the job demands, job control and social support measurements. Given the big sample size, instead of chi-square test the goodness of model was examined using other fit indices: CFI and TLI (> .95 indicates a good fit), RMSEA (< .06 indicates a good fit) and SRMR (< .08 indicates a good fit) (e.g., Brown, 2006;Hu & Bentler, 1999).
The measurement model had a good fit ( 2 (130) = 3173.089, p < .001, CFI = .921, TLI = .907, RMSEA = .055, SRMR = .049) when two correlations between items (Does your job require you to work fast and Does your job require you to work very hard; Do you get to learn new things in your job and Does your job require creativity) were allowed to correlate.

Statistical methods
In this study two different methods are used to examine relationships between work engagement and the three predictors. Generalized Additive Models (GAMs) are applied to estimate flexible curvilinear relations and General/-ized Linear Models (GLMs) are utilized to estimate linear relations, as has often been done in previous studies. Whereas previous studies have usually forced a relationship between work engagement and its main predictors by using linear forms, in this article we use GAM which lets us estimate the shape of the relationship from the data.
Whereas GLM specifies a designated functional form (such as a straight line or quadratic curve) to describe the relationship between response and predictor variables, GAM allows the data to determine the functional form of the relationship, which can be very flexible and curvilinear. In GAMs, smooth functions are represented by using penalized regression splines. Simon Wood (2006), who works as a professor of statistics at University of Bath, states that the estimation of GAM is straightforward because GAM is in fact simply a penalized GLM with a linear predictor involving a sum of the smooth functions of the predictor variable (Wood 2006).
In a GLM procedure, curvilinear relations can also be estimated. This can be done by adding quadratic terms to the model but this often leads to problems with multicollinearity and the curvilinear effects estimated this way are global when in some cases the curvilinearity is rather a local phenomenon, meaning that curvilinearity might occur only on a certain range of the scale (i.e. only at the high end of the scale or at the middle of the scale) (Beck & Jackman, 1998).
GLM is also able to estimate interaction effects between predictors, though interactions are usually forced into a linear form for the sake of simplification. In contrast, GAM can estimate simple multidimensional smooths which can be interpreted as interactions with possible curvilinearity.
In this study, the data analysis was performed by using statistical program R 2.10.1 (R Development Core Team, 2011) and a special mgcv package created by Simon Wood. In GAMs thin plate regression splines were used and an appropriate degree of spline smoothness was estimated from the data by using Generalized Cross Validation (GCV). To avoid overfitting, the model degrees of freedom were inflated with a constant multiplier (gamma = 1.4) -as recommended by Wood (2006). Additionally, the upper limit on the effective degrees of freedom associated with the smooths was set to four (cf. degree-four polynomial).
The differences between the GAM and GLM models were compared by using several statistical tools. Goodness of model fit was observed by using the measures of the proportion of the explained deviance (cf. R 2 ) and adjusted R 2 , AIC-score (Akaike Information Criterion) and the generalized likelihood ratio test (Deviance test). Large values of adjusted R 2 and the proportion of explained deviance indicate good model fit. AIC-scores of multiple models can be used in the model comparison. The model with the smallest value of AIC suggests a best fit to the data. The AIC and GCV-scores were used with a Deviance test to assess the significance of the interaction terms and the linearity of the predictors. The GCV-score is a good measure of model fit in GAM models. The interpretation of the GCV-score is similar to the AIC-score, a smaller score indicates a better fit (Wood 2006).
In total, we estimated seven different models (three GAMs and four GLMs). All models included the main effects of the control variables. One GAM and one GLM were used to estimate the main effects of job characteristics. Two GAM interaction models were built by adding an interaction term to the main effect model. The significance of the interaction effects was tested by comparing the main effect GAM to the interaction GAMs. In GLM the interaction terms were different to those in GAM. Interaction terms were estimated and interpreted with the method proposing by Aiken and West (1991). The simple slopes strategy was used in interpreting and visualizing the results. The effect of job demands on work engagement was investigated at three levels of job control and social support: one standard deviation above the average level, average level and one standard deviation below the average level.

Main effects of the predicting variables
Several statistical indicators show that GAM has a better fit to our data than GLM when estimating main effect models where work engagement was explained by control variables and with the main effects of job demands, job control and social support (see Table 1). GAM was superior to GLM in the proportion of explained deviance (+1.02), adjusted R 2 (+.010) and the AIC-score (-96). Also the deviance test indicated that GAM fits significantly (p < .001) better to our data than GLM. Therefore, the estimates of GAM are prioritized and used hereafter in this study when reporting the effects of job characteristics on work engagement. Table 1 here The main effects of job demands, job control and social support turned out to be curvilinear at a statically extremely significant level (p < .001). However, by looking at Figure 1, it can be noticed that the effects of job demands and job control were almost linear in practice. In relation to social support, GAM provides clear additional information when compared to the linear GML approach. Non-linear modeling reveals that there is a U-shape relation between social support and work engagement, which remains undiscovered when using the linear GLM procedure. It follows from this that Hypothesis 1 (H1) is clearly rejected especially regarding social support. To better understand the curvilinear relationship between social support and work engagement we conducted some further investigations. We noticed that only one item of the applied social support scale echoed the found curvilinear connection to work engagement while others had positive linear effects. That one item, which dominated the whole form of connection between social support and work engagement, was: "I get on well with my colleagues". This unexpected curvilinear connection requires more careful research as prior studies have come up with conflicting results. Rydstedt et al. (2006) have noticed a similar connection to ours, but between social support and job satisfaction. On the other hand, inverted U-shaped relationships between supervisor support and work engagement (Sawang et al. 2009) and social support and job satisfaction (De Jonge & Schaufeli 1998) have been discovered as well. Also De Jonge et al. (2000) found U-shaped relationship between social support and burnout. Table 2 shows that regression coefficients for control variables do not vary much between GLM and GAM procedures. The model intercepts differ in GLM and GAM because of different parameterization. Table 2 points out that men experience lower levels (-0.208) of work engagement than women. The level of work engagement is also higher for respondents over 50 years of age compared to younger respondents. Also respondents without children living at home (+0.073) register higher levels of engagement than those who have children at home. Supervisory positions (+0.151) are also associated with more engagement than non-supervisors.
Highly educated respondents seem to suffer from a lower level (-0.170)

Linearity and non-linearity in interactions
Interaction between job demands and job control was statistically significant in both GLM (see Table 3) and GAM (see Table 4) when the interaction between job demands and social support was statistically significant only in GAM. In GLM the regression coefficients for interaction effects were 0.205 (p<.001) for the interaction between job demands and control and -0.038 (p=.214) for the interaction between job demands and social support.
To sum up, our results suggest the rejection of H2 and H3 proposing that the interaction effects between job demands and control or between job demands and social support cannot be found by means of GML and GAM. Both interaction effects were found by means of GAM and the interaction between job demands and control also by GLM. However, the results make it evident that non-linear modeling paints a more rigorous picture of the interaction effects in question. In particular, GAM revealed that the interaction between job demands and social support is much more complex than is suggested by linear GLM.

Discussion
Prior studies have mostly indicated that there are linear relationships between job characteristics and work-related well-being. On the other hand, only a few studies (e.g. De Jonge 2000; Rydstedt et al. 2006) have so far assumed that these relationships might not be linear but are curvilinear in nature. We used the generalized additive model (GAM) to go beyond linear modeling and to obtain more profound information on the relationship between job characteristics and work engagement.

Theoretical implications
There is lot of empirical evidence to support the assumptions of the JD-R model regarding job demands and resources. Job resources (e.g. job control and social support) have a positive linear effect on work engagement (e.g. Mauno et al., 2007;Schaufeli & Bakker, 2004;Schaufeli et al., 2009;Xanthopoulou et al., 2009). Prior research also states that job resources are particularly important when worker faces high job demands (e.g. Kühnel et al., 2012).
In this study relationships between work engagement, job demands, job control and especially social support turned out to be more complex than other studies -premised on the generalized linear model (GLM) or other linear methods -have suggested. In our data, there were curvilinear relationships between work engagement and all investigated job characteristics. While connections between work engagement and job demands as well as work engagement and job control were only weakly non-linear (i.e. almost linear), there was a U-shape relationship between social support and work engagement. Previous studies have discovered similar or opposite (inverted U-shape) with social support and well-being (De Jonge et al., 2000;De Jonge & Schaufeli, 1998;Rydstedt et al., 2006). Very low levels of social support might indicate a completely isolated work (e.g. teleworking) where the role of social support on well-being is understandably different. Furthermore, the curvilinear relationships discovered do not support Warr's vitamin model, but the curvilinearity of job demands could be explained with the JD-R model. The non-linear effect of job demands might explain why prior studies based on linear methods have not confirmed the negative effect between job demands and work engagement.
Our study pointed out that the interaction effect of job demands and job control on work engagement was very similar in linear GLM and GAM, but GAM indicated the strong curvilinear interaction effect of job demands and social support on work engagement. Karasek's model and the JD-R model gain support regarding the buffer effect of job control and social support on job demands. Job resources are beneficial for workers and the results also suggest that job demands do not have a negative effect on well-being of the worker if the demands are not particularly high or they cannot be matched with high job resources.

Practical implications
This study encourages organizations develop work related well-being at the organizational and task related levels. For instance, employees with demanding work tasks might gain higher levels of well-being not through the reduction of their job demands but through the improvements in their resources, such as job control and social support, which could help converting demands into challenges. It also seems that certain groups of employees favour completely isolated work tasks to the work with only a moderate level of social support. This is the case especially if the work is not very demanding. The findings of the study encourage cultivating and fostering job resources in organizations to thrive work engagement.
This study also has a clear practical implication for researchers by pointing out that GAM is a good alternative to the commonly used GLM. Even though GLM is easy to use and it provides uncomplicated results. This study showed that the use of linear GLM can lead to a too straightforward and therefore skewed interpretation. While it is possible to examine curvilinear relationships with GLM simply by adding quadratic terms into the model, yet adding curvilinearity to the interaction effects is not that straightforward in GLMs (see Aiken & West, 1991). Therefore, GAM is a feasible tool with which to examine especially multifaceted interaction effects. GAM interactions are easy to visualize with three dimensional interaction plots that contain more information than traditional linear GLM interaction plots. GLM is a calculation-effective method which has probably been one reason for its popularity. The rapid development of computer technology has made GAM and other calculation-intensive nonparametric regression methods more available to researchers. GAMs can be estimated using for example the R-package mgcv, which is relatively easy to use and well documented.

Study limitations
Both our data and some strict presuppositions concerning the relationships between job characteristics and work engagement impose a few limitations on the results of this study. The data used in the study was cross-sectional and therefore conclusion about the causality between work engagement and its predictors could not be made. The data was also collected only from the service sector and therefore the results cannot be generalized to other sectors. Although GAM enabled the estimation of very flexible relationships, we made some simplifications. We entered the controlled variables as main effects into the model. By doing this, we were not able to take into account, for example, whether men and women might have different kinds of relationships between job characteristics and work engagement.
One limitation of the study related to the study variables. Even though the data was collected from many different organizations the responses regarding job autonomy, job demands and social support were mostly situated at the middle of the scales. The lack of respondents at the both high and low ends of the scales shows as wide confidence intervals in the Figure 1, which leaves the possible curvilinearity unclear. In future research it is important to collect data from different organizations to ensure sufficient variation of responses.

Opportunities for future research
The results of this study encourage researchers to study curvilinear relationships between work engagement and its predictors. This study focused only to job autonomy, job demands and social support while there are many more aspects that are related to well-being in work. Besides GAM future research can apply other non-linear methods in work engagement research. For example cusp catastrophe modeling can be utilized to estimate "catastrophe points", where minor changes in predictor values are connected to major changes in dependent variable.
GAM can be used not only in work engagement research but in various fields of study also.
According to this study GAM is most practicable in situations where there are reasons to assume curvilinear relationships, complex interactions or other moderator effects between predictors.
GAM can also be used in data-driven exploratory research where there is little or no prior knowledge about the phenomenon concerned. GAM is useful tool for collecting precise information, even in complex research settings. Yet, the other side of the coin is, GAM can overfit the data and make simple relationships look even more complex than they are.    (1) p=.214 Note: ***p < .001, **p < .01; *p < .05.
The interaction effects were investigated by adding interaction terms separately to the main effect model of work engagement.