# Simplifying variable effects in a GLM in R

Apologies, but it looks like my question is off topic for this forum. Thanks for all the excellent replies though. For those who have come across this question if they've been looking for something similar - then the short answer to my question below is very likely 'no'.

Please note this question has been edited in the light of the excellent responses below.

Can anyone recommend a way to easily ‘simplify’ the effects of variables in an R based GLM model? By simplify I mean any or all of the following:

• Group levels of a variable with similar effects.
• Apply a curve to the effects of an ordinal variable.
• Hand-smooth / alter the effects of a variable.
• Band (and re-band) a continuous variable.

Given that one has accepted the pros and cons of taking any of these approaches, are there any packages that might help with this? I’ve drawn a blank so far.

Thanks everyone for your responses.

For clarity (I hope) I have edited the question to remove the (inadvertent) confusing use of terminology , notably removing the use of ‘factor’, which I used in the original wording as a general term for ‘predictor variable’ rather than in the particular R sense of a variable with specific levels (something which as a relatively new user of R I unfortunately overlooked).

Also please note that I’m asking this question from an actuarial rather than a pure research standpoint.

In this field, there are reasonably well understood risk profiles (the chance of having a car accident, illness or death at different age for instance) and also a degree of leeway in applying ‘reasonable’ tweaks to the effects of individual variables (to the extent that most practitioners end up using software specifically designed for this purpose - a package called ‘Emblem’ being pretty much ubiquitous in the insurance industry in the UK).

I’m aware of the costs and benefits of treating any results in this way, but was just wondering if there were methods in R that could aid the process once the pros and cons are accepted.

Original Question:

Can anyone recommend a way to easily ‘simplify’ factors in an R based GLM model? By simplify I mean any or all of the following:

• Combine levels of a factor with similar effects.
• Apply a curve to the effects of a factor. Hand-smooth / alter the effects of a factor.
• Band (and re-band) factors in the model.

Are there any packages that might help with this? I’ve drawn a blank so far.

• What does "band" (&/or "re-band") mean? – gung Oct 13 '14 at 16:05
• Note: although this looks like a question asking for code / packages in R, I think the question is motivated by a lack of understanding of the statistical issues associated w/ categorical variables. As a result, I think this Q should stay here (not be closed or migrated). – gung Oct 13 '14 at 16:07
• What does it mean to "band" (&/or "re-band") a continuous variable? Also, what do you mean by "hand smooth / alter the effects of a variable"? – gung Oct 14 '14 at 12:59
• Note that asking for packages or code is off-topic, both here & on SO. You might try asking on the r-help listserv. As for the broadness of your questions about modeling strategies, it is material that I cover in multiple semesters worth of stats courses. – gung Oct 14 '14 at 13:49
• I appreciate your respect for the rules of this site. I would say that 'hand-smoothing' as you describe it is not valid. Your best bet would be to read a thorough text on regression (I like Kutner et al., 2005) & then read Harrell's Regression Modeling Strategies. – gung Oct 15 '14 at 14:31

It's a simplification to have fewer parameters, but simplifying a model in the light of particular results has a cost too. Statistically-minded people differ in their views on modifying a model in the light of model results. Opposite arguments can both have merit: it's natural to want a model to be parsimonious as well as a good fit to the data (what else do we want?), but it's also risky to capitalise on chance and be over-responsive to what may be idiosyncrasies in a particular dataset. Some would argue that you are spending degrees of freedom in making such modifications; the question is whether the accountancy is honest and explicit. Whether there is a specific significance test that justifies a choice is not quite the issue here.

What is often neglected is the importance of fitting a model in a form that is both more transparent to the reader and easy for other researchers to compare given different data.

So consider this kind of statement:

"Factor $X$ as initially coded had 3 levels. But on an initial fit levels 2 and 3 were similar in their effects, so we merged them and the published model is for a two-level version."

(In some fields, you might be lucky to get this kind of explanation, as the project might be written up as if the final coding was that determined at the outset.)

There are at least two problems with this other than those already hinted at (discussed in various literatures under headings such as "data snooping"):

1. If it was thought worthwhile to use 3 levels in data production (e.g. these were the categories offered to people taking a questionnaire survey) there is information in the pattern of the coefficients, even if some seem similar in magnitude in one particular dataset.

2. This kind of adhockery makes it difficult for other researchers to do similar studies and check whether patterns are similar with other datasets. It's not obvious that the simplification defensible in one dataset will make sense with another dataset.

All of this implies, as a counsel of perfection, being open about initial and final models and why a model was modified. In practice there often need to be compromises and may be severe constraints depending on (e.g.) instructions from supervisors, conventions on reporting in different fields, or pressures from reviewers and editors on allowed space in journals.

This answer deliberately does not focus on what software might be used, within R or within any other statistical software, and why can now be explained. Model simplification does, and should, depend on substantive knowledge on what makes sense in the light of the underlying science, or in the light of practical knowledge about the situation in which the data were produced or in which the model is to be used. It is precisely the kind of decision that is, and should be, difficult to automate. That prejudice is, naturally, not a proof that software purporting to help is impossible or non-existent.

• Thanks for your comment it's an excellent summary of the pros and cons of the kind of 'simplification' I was asking about. I've edited the question to include a little more about the standpoint from which the question was asked. I very much agree with the point that any such 'simplification' should be well documented and accounted for. – AdrianD Oct 14 '14 at 9:42
• You're still pressing for R software to do this. I only use R occasionally so I can only advise searching the usual R places. Also, many users here have R as their main or indeed only statistical software, and although absence of evidence is not evidence of absence, no one has mentioned any such. – Nick Cox Oct 14 '14 at 10:29

(@Nick Cox has done a good job talking about the issues involved in combining levels of a factor; I thought I would explain a little bit about categorical variables.)

It is not possible to 'fit a curve' to a categorical (nominal) variable. Nor is it possible to 'smooth' the fit. There are many ways to represent / code a factor (also here). Nonetheless, in one way or another, all methods use $k$ degrees of freedom to represent $k$ groups. In the end, the groups are fit by finding $k$ group means. By contrast, it is possible to have $k$ groups at fixed locations along a single, continuous variable (such as $0$ mg, $10$ mg, $20$ mg, etc.). In that case, a single degree of freedom can be expended to fit a straight line, or more for curved functions (e.g., $2$ for a curve with a single 'bend'). The implication of using $k$ degrees of freedom for $k$ groups is that each group is placed on its own dimension in the model's space. With, at most, $2$ groups on any dimension, curved functions are not meaningful. Thus, the reason you cannot find any R packages to do this is because it cannot be done. To illustrate, consider the following:

cont.x = c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3)
cat.x  = as.factor(c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3))
y      = c(1,2,3,4,5, 3,4,5,6,7, 5,6,7,8,9)


In both the continuous and categorical cases, what is entered is 1,2,3, but let's inspect the model matrices (a model matrix is how the variables are represented 'behind the scenes') and the resulting model fits:

model.matrix(lm(y~cont.x))
#    (Intercept) cont.x
# 1            1      1
# 2            1      1
# 3            1      1
# 4            1      1
# 5            1      1
# 6            1      2
# 7            1      2
# 8            1      2
# 9            1      2
# 10           1      2
# 11           1      3
# 12           1      3
# 13           1      3
# 14           1      3
# 15           1      3
# attr(,"assign")
# [1] 0 1
summary(lm(y~cont.x))
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   1.0000     1.0377   0.964  0.35282
# cont.x        2.0000     0.4804   4.163  0.00111 **
#  ...


The continuous x is represented by a single variable. In the model fit, the slope of its line is given and the intercept is just the y value where the line crosses x=0. For the categorical x, it's different:

model.matrix(lm(y~cat.x))
#    (Intercept) cat.x2 cat.x3
# 1            1      0      0
# 2            1      0      0
# 3            1      0      0
# 4            1      0      0
# 5            1      0      0
# 6            1      1      0
# 7            1      1      0
# 8            1      1      0
# 9            1      1      0
# 10           1      1      0
# 11           1      0      1
# 12           1      0      1
# 13           1      0      1
# 14           1      0      1
# 15           1      0      1
# attr(,"assign")
# [1] 0 1 1
# attr(,"contrasts")
# attr(,"contrasts")\$cat.x
# [1] "contr.treatment"
summary(lm(y~cat.x))
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)   3.0000     0.7071   4.243  0.00114 **
# cat.x2        2.0000     1.0000   2.000  0.06866 .
# cat.x3        4.0000     1.0000   4.000  0.00176 **
# ...


A separate variable is created for x=2 and for x=3, x=1 is indicated by 0's in both of those variables. In the model fit, the intercept now represents the mean of the first group, and the following two variables give the difference between the latter groups (each individually) and the first group.

Here are some pictures of what they look like:

• I agree with the main point, naturally, but in practice many researchers are broad-minded (read: cavalier) about measurement scales and flip at whim between e.g. regarding a 5-point scale as ordinal and as interval. So, your point is perhaps "If you really regard a factor as such, treating it as measured is out of the question". – Nick Cox Oct 14 '14 at 1:26
• @NickCox, you raise a reasonable point, but I suspect you are thinking of the factor as ordinal; I am referring to nominal / categorical factors. I can also explain ordinal factors, but it's much more involved. – gung Oct 14 '14 at 1:29
• Indeed; to treat nominal as interval would be so crazy I did not even imagine it. A researcher wanting to treat nominal scales as measured needs some considerable advice! It's the ordinal/interval divide that is often jumped and that I was indeed thinking of. NB I don't feel wedded to exact R terminology of what is/isn't a factor; "factors" to me are categorical variables declared to software as such. – Nick Cox Oct 14 '14 at 1:33
• @NickCox, I think that I (& R) mean the same thing by 'factor' that you do. I think that the OP is referring to nominal variables, though. If the factor were ordinal, R would return several variables (linear, quadratic, etc) based on orthogonal polynomials by default, so there's no need to ask about this. (Note that nominal groups are often listed as 1, 1, 1, ..., 2, 2, ....) – gung Oct 14 '14 at 1:39
• Thanks for the comments. When I wrote the original question I too did not consider that it might be assumed that I was considering apply a curve to a nominal variable, but can see that my slapdash use of 'factor' might have given that impression. – AdrianD Oct 14 '14 at 10:04