It should be a straightforward regression model of the form Y = b0 + b1X1 ... +bkXk + e

I'm guessing that you have identified the problem as automatic coding of the categorical predictors using dummy variables. You can force a more sensible coding by fitting linear contrasts.

e.g.,

hsb2

It should be a straightforward regression model of the form Y = b0 + b1X1 ... +bkXk + e

I'm guessing that you have identified the problem as automatic coding of the categorical predictors using dummy variables. You can force a more sensible coding by fitting linear contrasts.

e.g.,

hsb2

I think that in R ordinal factors are turned into orthogonal polynomials - see

http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm

Yes, a GLM with normal probability model (and identity link) is nothing but a LM.

Yes, factors (R's representation of categorical variables) in models are dummy-(reference-)coded (unsell stated differently). The first level of the factor is used as reference level (-> goes to the intercept), and the other coefficients (=betas) model differences to this reference.

You get the estimates of all betas (and the intercept) when you call

summary(your_model)

or only the estimates (=beta-hats) without further statistics with

coefficients(your_model)

If the model was specified including, for example

... ~ ... + TASTE + ...

where TASTE is a factor with the levels "bad","good","so.so" (in this order), the level "bad" is used as reference, so the intercept refers to the (expected) response for TASET=="bad".  The names of the coefficients will include

TASTEgood

TASTEso.so

and their values represent the (expected) differences in the response between TASTE=="good" vs. TASTE=="bad" and of TASTE=="so.so" vs. TASTE=="bad".

You can see the design matrix including the dummy variables created by R with

model.matrix(your_model)

Thom, thanks for that.  I don't have a copy of the code on my PC so unfortunately I can't post.  If it's equivalent I think I'll have him use the LM code since I'm much more familar with it and know how to specific contrasts etc.

Kelvyn, I think you're probably right.  I seem to recall the coefficients being labelled as L, Q and C which confused me at the time.  This will make manually displaying the regression equation problematic since it is no longer simply the X score multiplied by the beta coefficient, but rather a trend fitted to the X score.  Do you know how we would write out the regression equation in full if it contains a linear polynomial?

Jochen, I didn't know about the model.matrix() command, that will be useful thanks.

Garvin, you can use R to reate a character string containing the fitted regression equation as a text. Here is an example based on the data "Toothgrowths" (comes with R, is a data.frame with a metric and a categorical predictor):

model

Here are some orthogonal polynomials for 3 , 4 and 6 categories 

Categorical Variables Codings for 3 categories

                  Linear Quadratic

3 groups 1 -.707 .408

               2 .000 -.816

               3 .707 .408

Categorical Variables Codings for 4 categories

                   Lin   Quad  Cubic

4 groups 1 -.671 .500 -.224

               2 -.224 -.500 .671

               3 .224 -.500 -.671

                4 .671 .500 .224

Categorical Variables Codings for 6 categories

Linear Quadratic Cubic 4th order 5th order

6 groups 1 -.598 .546 -.373 .189 -.063

               2 -.359 -.109 .522 -.567 .315

               3 -.120 -.436 .298 .378 -.630

               4 .120 -.436 -.298 .378 .630

               5 .359 -.109 -.522 -.567 -.315

               6 .598 .546 .373 .189 .063

Thank you both for your responses.  I apologize but I'm afraid I'm still not getting it.  I do not have any experience of applying polynomial terms to predictors in a regression model in this way.  I think those are the contrast weightings you've provided me with Kelvyn which (I believe) don't appear directly in the regression equation themselves.

The standard multiple regression model is Y = B0 + B1*X1 + B2*X2..... +Bi*Xi

What happens to this equation when there are polynomial terms applied to the predictors?

My understanding is that B1*X1 IS a linear polynomial term so if you've modeled that it doesn't change how you deal with things.  To turn that into a quadratic term it would become B1*X1^2, to turn it into a cubic it becomes B1*X1^3 and so on - am I getting this right?

So as an example a model which contained 1 predictor which had linear, quadratic and polynomial terms applied to it, and 1 "standard" continuous predictor, and no interaction term

Y = B0 + B1*X1 + B1*X1^2 + B1*X1^3 + B2*X2

Is this right?

If you fit the model in R using

x1quad

Add-on:

Fitting the model instead as

model = lm(Y ~ poly(x1,degree=3) + X2)

will fit othogonal polynomials. The coefficients refer to the transformed values of X1, X1², and X1³ so that the covariance of the transformed data is 0. So the coefficients you get from this fit are referring to transformed X1-values, and you can not simply use X1-values on their original with these coefficients to predict Y.

There are two options (I know) to get is anyway:

(1) use poly(X1, degree=3, raw=TRUE)

to use the poly-function but without transforming the X1-values. Now the coefficients are interpretable on the original scale of x.

(2) back-transform the coefficients from the model to match the original scale. This is a bit more complicated. See the attached link for the solution in R.

Here is an example:

attach(cars)

m1 = lm(dist ~ speed + I(speed^2))

m2 = lm(dist ~ poly(speed,2,raw=TRUE))

m3 = lm(dist ~ poly(speed,2))

# comparison of the coefficients:

cbind(manual=coef(m1),raw.poly=coef(m2),orth.poly=coef(m3))

# all three models give the same predictions:

newSpeed = data.frame(speed=10:15)

cbind( newSpeed,

manual=predict(m1, new=newSpeed),

raw.poly=predict(m2, new=newSpeed),

orth.poly=predict(m3, new=newSpeed))

# equation is dist = 2.47 + 0.91*speed + 0.01*speed^2

# using the function poly2raw from the attached link:

poly2raw( poly(speed,2), m3$coefficients )

# --> 2.47, 0.91, 0.01

https://dennisphdblog.wordpress.com/2010/02/01/extracting-the-mean-polynomial-from-an-r-poly-regression/

Apologies for not being clear

If you specify a 4 category ordinal variable in a regression,  R turns this into a constant ( a set of 1's ie X0) which will  gives the intercept term (B0) and X1 which is the linear orthogonal polynomial and X2 which is the quadratic and X3 which is the cubic.

So say your factor variable is Class and an extract is

Person  Class

1            1st

2            2nd

3            4th

4             2nd

5             3rd

500       1st

These data will become the following when factor(Class) is used

Person  Class       X0        X1            X2             X3  

1            1st           1         -.671          .500        -.224 

2            2nd          1          -.224         -.500        .671

3            4th           1           .671           .500         .224

4             2nd         1          -.224         -.500        .671

5             3rd          1          .224          -.500        -.671

500         1st           1           -.671          .500        -.224 

So it fits  Yhat = B0X0 +B1X1 + B2X2 + B3X3

You can then possibly use a more parsimonious form if for example B1 is a  large effect while B1 and B2 are small then this suggests an underlying linear trend across Classes   and you could drop X2 and X3 to have fewer terms in a simplified model. You do not have to refit the model as the X's are designed to be orthogonal so the linear trend will be given by B1.

An example is given in this manual

https://www.researchgate.net/publication/234015221_Manual_supplement_for_MLwiN_Version_214

Data Manual supplement for MLwiN Version 2.14

Thank you for your responses.  The impression I'm getting is that making predictions from a model which contains both elements of polynomial regression and standard linear regression is not a trivial problem.  For our purposes it would probably be simpler to simply avoid using the polynomial transformations on the predictors and see if something simpler will do the job.

Possibly - but you can also use the predict() function to get the predictions or perhaps even easier use the simulate() function to simulate what is going on. The predict() function is a bit fiddly, though there are some packages that make it easier ISTR.

OK I met with my student today.  His finished model contained a Linear, Quadratic and Cubic function applied to one of his predictors.  We input some test values using predict() and as expected we couldn't figure out how R was arriving at its predictions.

So I computed 3 new variables - the ordinal variable converted to an integer (i.e. the linear), this value squared (quadratic) and this value cubed (cubic).  I removed the original predictor and inserted these 3 new ones in its place.

The predictions generated by these two models agreed with one another down to about 9 decimal places so I think we can say they are equivalent.

The new model gave us coefficients which can be used with the equation of the form:

Y = B0 + B1*X1 + B1*X1^2 + B1*X1^3 + B2*X2 .....

So not only can we now make predictions, we know exactly where those predictions are coming from.

Thanks

Gavin, could you show the R call how the model was fit?

This could help to understand why you "couldn't figure out how R was arriving at its predictions" and how to get the equation right from that model (without using new variables).

Jochen, we did the work on my students laptop so I do not have a copy.

As I discussed above I was struggling to understand what to do with the slope coefficients for a variable to which polynomial terms have been applied.

The approach I've adopted here (which I got from your 3rd post, thanks) produces essentially identical predictions so it meets our needs.  I believe it is doing the same thing as R does automatically by a slightly different route.

OK I think I have figured it out.  I apologize if you have all been trying to explain this and I've just not been getting it.

I created some mock data with 1 continuous predictor and 1 categorical predictor which has polynomial contrasts applied and ran an LM on it.

If the categorical predictor has 4 levels then the polynomial contrast weightings are:

    .L                          .Q          .C

1 -0.670820393     0.5     -0.223606798

2 -0.223606798     -0.5     0.670820393

3 0.223606798     -0.5     -0.670820393

4 0.670820393     0.5      0.223606798

The regression coefficients provided by my model are:

(Intercept)   60.40375308

Group.L       -31.45423223

Group.Q      -5.025481752

Group.C      28.71857561

ContPred   -3.071113661

Therefore, when the group is equal to 1, the regression equation is:

Y = Intercept + (ContPred*-3.07) + (-0.67*-31.45(this is the linear component)) + (0.5*-5.03 (quadratic component)) + (-0.22*28.71 (cubic component))

When the groups are 2, 3 or 4 the first part of the linear, quadratic and cubic components are replaced with the appropriate value from the contrasts weighting table.  I've checked my working on this and it returns predictions which are within 0.00000001 of what R provides for  a given value which is well within rounding error.