Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down at the end of May, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: Calculation of cubic splines


From   Roger Newson <r.newson@imperial.ac.uk>
To   "statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu>
Subject   Re: st: Calculation of cubic splines
Date   Thu, 19 May 2011 14:45:38 +0100

Another possibility might be to use the -flexcurv- option of the -bspline- package, downloadable from SSC. The -flexcurv- module creates a spline basis whose corresponding parameters are the values of the spline at reference points on the X-axis or (alternatively) differences between the values of the spline at each reference point and the value of the spline at a baseline reference point. For differences, of course, you can substitute ratios, if your model is a generalized linear model with a log link, or a linear regression model on the logs. The -bspline- package comes with a .pdf manual, as well as on-line help.

I hope this helps.

Best wishes

Roger


Roger B Newson BSc MSc DPhil
Lecturer in Medical Statistics
Respiratory Epidemiology and Public Health Group
National Heart and Lung Institute
Imperial College London
Royal Brompton Campus
Room 33, Emmanuel Kaye Building
1B Manresa Road
London SW3 6LR
UNITED KINGDOM
Tel: +44 (0)20 7352 8121 ext 3381
Fax: +44 (0)20 7351 8322
Email: r.newson@imperial.ac.uk
Web page: http://www.imperial.ac.uk/nhli/r.newson/
Departmental Web page:
http://www1.imperial.ac.uk/medicine/about/divisions/nhli/respiration/popgenetics/reph/

Opinions expressed are those of the author, not of the institution.

On 19/05/2011 14:08, Maarten Buis wrote:
You already have the predicted probabilities, that is what you got
when you typed -predict p_spline-. The analytic representation is not
so easy, so that I cannot meaningfully start trying to type it here in
plain text. I would just look at the methods and formulas section of
the manual entry on -mkspline-(*). Alternatively I often like linear
splines as a good compromise between interpretability and flexibility
of the curve.

Hope this helps,
Maarten

(*) Note however that when you are using Stata 10, the manual entry
contains a typo in the formula. I believe a minus sign on the wrong
side of one of the brackets.

On Thu, May 19, 2011 at 2:57 PM, Mikkel Brabrand<mikkel@brabrand.net>  wrote:
All.

I am trying to use cubic splines to assess risk of in-hospital mortality using some vital signs. I would like to calculate the predicted mortality using cubic splines manually.

I have defined the following knots:
. mkspline _Ssbt = sbt, cubic nknots(5) displayknots

             |     knot1      knot2      knot3      knot4      knot5
-------------+-------------------------------------------------------
         sbt |        97        119        132        146        178

. mat sbt_knots = r(knots)

. mkspline _Stemp = temp, cubic nknots(5) displayknots

             |     knot1      knot2      knot3      knot4      knot5
-------------+-------------------------------------------------------
        temp |      35.9       36.6         37       37.3       38.8

. mat temp_knots = r(knots)

. mkspline _Salder = alder, cubic nknots(5) displayknots

             |     knot1      knot2      knot3      knot4      knot5
-------------+-------------------------------------------------------
       alder |        23         53         66         76         88

. mat alder_knots = r(knots)

And have run the logistic regression as follows:

. xi: logit in_hosp_mort _Ssbt* _Stemp* _Salder*, or

Iteration 0:   log likelihood = -364.70483
Iteration 1:   log likelihood = -322.05257
Iteration 2:   log likelihood = -301.68143
Iteration 3:   log likelihood = -299.60534
Iteration 4:   log likelihood =  -299.5029
Iteration 5:   log likelihood = -299.50192
Iteration 6:   log likelihood = -299.50192

Logistic regression                               Number of obs   =       2979
                                                  LR chi2(12)     =     130.41
                                                  Prob>  chi2     =     0.0000
Log likelihood = -299.50192                       Pseudo R2       =     0.1788

------------------------------------------------------------------------------
in_hosp_mort | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      _Ssbt1 |   .9557737   .0141873    -3.05   0.002     .9283678    .9839888
      _Ssbt2 |    1.19366   .1618851     1.31   0.192     .9150398    1.557117
      _Ssbt3 |   .4041954   .3164691    -1.16   0.247     .0871232    1.875205
      _Ssbt4 |   3.502106   4.366206     1.01   0.315     .3041618    40.32311
     _Stemp1 |   .2456296   .0807298    -4.27   0.000     .1289796    .4677788
     _Stemp2 |   6.086799   31.49901     0.35   0.727     .0002396    154643.8
     _Stemp3 |   1.01e+11   3.39e+12     0.75   0.452     2.19e-18    4.62e+39
     _Stemp4 |   3.34e-36   2.06e-34    -1.33   0.185     1.19e-88    9.42e+16
    _Salder1 |   1.078022   .0906853     0.89   0.372     .9141613    1.271254
    _Salder2 |   1.004445   .1631946     0.03   0.978     .7305158    1.381093
    _Salder3 |   .7612758   .8178836    -0.25   0.800     .0926928    6.252275
    _Salder4 |    2.45316    5.18604     0.42   0.671     .0389283    154.5918
------------------------------------------------------------------------------

. predict p_spline if e(sample)
(option pr assumed; Pr(in_hosp_mort))
.
. roctab in_hosp_mort p_spline, summary

                      ROC                    -Asymptotic Normal--
           Obs       Area     Std. Err.      [95% Conf. Interval]
         --------------------------------------------------------
          2979     0.8207       0.0240        0.77361     0.86786

.

My question is: What is the formula I should use to calculate the predicted mortality? I have spend a great deal of time on this and have not been able to figure it out.

Thanks.

Mikkel
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/




*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index