Piecewise Linear regression graph

[Megan Farris]

Hello,
I have this data that I created these dichotomized variables for two separate relationships of the log odds of CHD for chol (cholesterol continuity at 280) with smoke and age as confounder variables. I am not sure who to create this piecewise linear regression graph.


gen chol280= (chol>280)
gen chol2= chol280*(chol-280)
gen chol1= chol-chol2
logit chd chol1 chol2 smoke age
logit chd chol1 smoke age

Can anyone lead me in the right direction?

[Aspen Chen]

The first three lines can be accomplished simply with

mkspline chol1 280 chol2=chol

I am not sure what the second logit model was intended for. For the first, this thread provides an example that should be adaptable to your need.

[Maarten Buis]

Is this what you are looking for?

//load and prepare some example data
sysuse auto, clear
recode rep78 1/2=3

// make splines
mkspline mpg1 20 mpg2 = mpg

// estimate the model
logit foreign mpg1 mpg2 price i.rep78

// change the data in a way you don't want to save
preserve
// fix price at the mean
sum price if e(sample)
replace price = r(mean)

// fix rep78 at acceptable
replace rep78 = 3

// predict the probabilities
predict p

// graph the probabilities
twoway line p mpg, sort

// get the original data back
restore

[Marcos Almeida]

With regards to graphical approaches, you may also take profit of the - margins, at() - commands.

Hopefully it helps.

Best,

Marcos

[Megan Farris]

Thanks for the advice, My graph does not look correct though. I carried out these commands:

mkspline chol1 280 chol2=chol
logit chd chol1 chol2 smoke age
preserve
sum age if e(sample)
replace age = r(mean)
replace smoke = 1
predict p
twoway line p chol, sort

I think the issue is that I need to graph chol1 and chol2 as separate relationships but keep continuity in the line. So chol1 will have one slope and chol2 will have a different slope but they will be joined at a knot at a chol levee of 280.

So I think my twoway command needs to account for the different slopes. I am just not sure how to accomplish this.

[Marcos Almeida]

Please try something like:

. margins chol1, at(age=40(5)80)) predict(xb)
. marginsplot, noci

Hopefully it helps.

Best,

Marcos

[Clyde Schechter]

I think the issue is that I need to graph chol1 and chol2 as separate relationships but keep continuity in the line. So chol1 will have one slope and chol2 will have a different slope but they will be joined at a knot at a chol levee of 280.

Well, we can't see what kind of graph you actually got. But your use of the word "slope" leads me to think that the problem lies in your expectations. The code you used looks correct to me. Remember that you are graphing p, the predicted probability, against chol. That relationship is not linear, not even piecewise linear--so there really is no "slope" except locally at any given point. If chol had been entered as a single variable, rather than as a spline, you would not see a straight line plot: you would see a sigmoid curve (or part of one depending on the range of values in your data). With a spline, you will get something that looks like pieces of two sigmoid curves colliding at a join point. Depending again on the range of values in your data, and where along the sigmoid curves they fall, this could take on a rather large variety of appearances.

If you want something that looks piecewise linear you have to -predict- and plot xb, not p. To get a feel for the possibilities try this:

set obs 100
gen x = _n
gen y = 2*x if x < 50
replace y = 5*(x-50) + 100 if x >= 50
graph twoway line y x, name(piecewise_linear, replace)
gen z = invlogit(y)
graph twoway line z x, name(right_angle_bend, replace)
replace y = y/100
replace z = invlogit(y)
graph twoway line z x, name(almost_linear_then_arcs, replace)

Here y plays the role of your model's xb, and z plays the role of the corresponding p's.

[Megan Farris]

Oh I see, thank you for clarifying.

[Richard Williams]

With regards to graphical approaches, you may also take profit of the - margins, at() - commands.

Hopefully it helps.

Best,

Marcos

I haven't tried it, but I have a feeling margins would get confused because it would not understand that chol, chol1 and chol2 are all related to each other. One of my wish list items for Stata 14 is more powerful factor variables, but spline functions might be especially touch to program.

[Marcos Almeida]

This is an example of – margins, at() – in logistic regression under piecewise models, taken from the date set gss_ivrm, according to the book: Mitchel, Michael N. Interpreting and Visualizing Regression Models Using Stata, p. 494.

. codebook smoke

----------------------------------------------------------------------------------------------------------------------------------
smoke                                                                                         Does R smoke (yes=1 no=0)? (recoded)
----------------------------------------------------------------------------------------------------------------------------------

                  type:  numeric (byte)
                 label:  yesno

                 range:  [0,1]                        units:  1
         unique values:  2                        missing .:  0/55087
       unique mv codes:  2                       missing .*:  38714/55087

            tabulation:  Freq.   Numeric  Label
                         10648         0  0. No
                          5725         1  1. Yes
                         38680        .i  
                            34        .n  

. mkspline edprehsm 12 edhsm 16 edcom = educ, marginal

. logit smoke c.edprehsm c. edhsm c.edcom hsgrad cograd age, nolog

Logistic regression                               Number of obs   =      16274
                                                  LR chi2(6)      =     806.95
                                                  Prob > chi2     =     0.0000
Log likelihood = -10136.225                       Pseudo R2       =     0.0383

------------------------------------------------------------------------------
       smoke |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    edprehsm |   .0540074   .0146957     3.68   0.000     .0252044    .0828105
       edhsm |  -.1571181   .0269648    -5.83   0.000    -.2099681   -.1042682
       edcom |   .0589168   .0420088     1.40   0.161    -.0234189    .1412526
      hsgrad |   -.599279    .062409    -9.60   0.000    -.7215984   -.4769597
      cograd |  -.3026172   .0953092    -3.18   0.001    -.4894197   -.1158147
         age |  -.0201919    .001057   -19.10   0.000    -.0222636   -.0181202
       _cons |   .2709492   .1535817     1.76   0.078    -.0300655    .5719638
------------------------------------------------------------------------------

. margins, at(edprehsm=0 edhsm=0 hsgrad=0 cograd=0) at(edprehsm=12 edhsm=0 edcom=0 hsgrad=0 cograd=0) at(edprehsm=16 edhsm=4 edcom
> =0 hsgrad=1 cograd=0) at(edprehsm=16 edhsm=4 edcom=0 hsgrad=1 cograd=1) at(edprehsm=20 edhsm=8 edcom=4 hsgrad=1 cograd=1) predic
> t(xb) noatlegend

Predictive margins                                Number of obs   =      16274
Model VCE    : OIM

Expression   : Linear prediction (log odds), predict(xb)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         _at |
          1  |  -.6257571   .1342834    -4.66   0.000    -.8889478   -.3625665
          2  |   .0119962   .0560904     0.21   0.831     -.097939    .1219315
          3  |  -.9997256   .0792405   -12.62   0.000    -1.155034    -.844417
          4  |  -1.302343   .0542245   -24.02   0.000    -1.408621   -1.196065
          5  |  -1.479118   .1183285   -12.50   0.000    -1.711038   -1.247199
------------------------------------------------------------------------------

[Marcos Almeida]

Below, I'm trying to display the graphic after - marginsplot -. according to the model above-mentioned.

I just excluded the labels on the x axis and recoded them into numbers for a better visualization.

[Richard Williams]

Thanks Marcos. So you can use margins, but you have to plug in the correct values for the spline variables yourself. Which is kind of a nuisance, but I imagine it would be non-trivial to extend factor variables to this kind of situation.

[Michael N Mitchell]

Thanks to Marcos for suggesting the example from my book! Unfortunately, the results from the -margins- command need to be massaged a little bit to create a more sensible graph. Following the -margins- command, try adding these commands:

mat yhat = r(b)'
mat educ = (0 \ 12 \ 12 \ 16 \ 16 \ 20)
svmat yhat
svmat educ
graph twoway line yhat1 educ1, xline(12 16) title("Piecewise Model")

graph export pw.png, replace

This yields the following graph:

Chapter 4 (namely, section 4.9) of my book explains how this code works to create the graph above. For more information, see http://www.stata.com/bookstore/inter...ession-models/

Warmly,

Michael N. Mitchell