Extracting a linear predictor formula from GEE model

Skip to first unread message

Pedro Emmanuel Alvarenga Americano do Brasil

Jun 18, 2021, 6:52:25 PM6/18/21
to MedStats-list
Hello StatMastars, 

I'm coming enlightenment. Im sure this question is something very simple that Im missing. I would like to add in a report the linear predictor formula from a gee model. There is a Hmisc::Formula that extract the formula but it does not work for geeglm objects. Additionally, the model has a spline term and after a while looking around on how to extract the formula I was not able to do so.  

fitar1_11 <- geeglm(outcome ~ v1 + v2 + v4_recoded + v5_trunc + bs(weeks, df = 3), family = binomial(), data = M2, id = f1p0,  corstr = "ar1")
geeglm(formula = outcome ~ v1 + v2 + v4_recoded + v5_trunc +
    bs(weeks, df = 3), family = binomial(), data = M2, id = f1p0,
    corstr = "ar1")
                           (Intercept)                                     v1
                         -1.586386e+00                           8.363458e-02
                                    v2 v4_recodedNET, NGT, JF, OET, GF or TPN
                         -3.763166e-02                           4.503599e+15
                              v5_trunc                     bs(weeks, df = 3)1
                          7.261799e-02                           3.469451e-03
                    bs(weeks, df = 3)2                     bs(weeks, df = 3)3
                          3.756985e-01                           4.542696e+00

Degrees of Freedom: 472 Total (i.e. Null);  464 Residual
> nd # newdata to predict at
  v1 v2  v4_recoded v5_trunc weeks
1  8  3 No catheter       36     2
> predict(fitar1_11, type = "response", newdata = nd) # prediction from the model
# Not working prediction from the handmade formula
> plogis(-1.586386 + 0.08363458*nd$v1 + -0.03763166*nd$v2 + 4503599142291968*0 + 0.07261799*nd$v5_trunc + 0.003469451*nd$weeks + 0.3756985*pmax(nd$weeks-1,0)^3 + 4.542696*pmax(nd$weeks-8,0)^3)
[1] 0.8772541

What am I missing to make this formula to make correct predictions? Any further reading would also be nice.

Best regards,  

Pedro Brasil

Martin Holt

Jun 20, 2021, 12:22:55 PM6/20/21
to MedStats
Actually, Stata is more intuitive ...


Mitchell Maltenfort

Jun 20, 2021, 12:27:11 PM6/20/21
to meds...@googlegroups.com
The spline term is the problem.    Spline coefficients aren't easily put in a formula.    Use poly with the raw option instead of bs.

To post a new thread to MedStats, send email to MedS...@googlegroups.com .
MedStats' home page is http://groups.google.com/group/MedStats .
Rules: http://groups.google.com/group/MedStats/web/medstats-rules

You received this message because you are subscribed to the Google Groups "MedStats" group.
To unsubscribe from this group and stop receiving emails from it, send an email to medstats+u...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/medstats/CAFfGvyKnpprHVyLH6z6bRoFQNaOa-Djrvt7nk_zKB8X84jc9OQ%40mail.gmail.com.
Sent from Gmail Mobile

Marc Schwartz

Jun 21, 2021, 8:32:14 AM6/21/21
to meds...@googlegroups.com, Pedro Emmanuel Alvarenga Americano do Brasil

If you want to stay with regression splines, you might want to look at this package on Github and the two SO threads that are pointed to in the "Story" section of the README file:


The functions in that package appear to enable you to extract the piecewise polynomials that are fit on either side of the spline knots. I might also suggest that you consider using ns() rather than bs(), if you want a base R function that more closely replicates what Frank's rcs() function in his rms package creates, where the spline function is linear in the tails.

Here are some additional readings on regression splines in R that you might find of interest:



Also, I think that you wanted to reference Frank's Function() in Hmisc rather than Formula() below, to extract the regression equation.

There is also:


which appears intended to generalize what Frank's Function() does to a more broad class of models. However, it still appears to be in development, and does not yet support GEE models nor splines/poly fits.


Marc Schwartz

Pedro Emmanuel Alvarenga Americano do Brasil wrote on 6/18/21 6:52 PM:

Karl Ove Hufthammer

Jun 21, 2021, 2:33:25 PM6/21/21
to meds...@googlegroups.com
'Marc Schwartz' via MedStats skreiv 21.06.2021 14:32:
> Also, I think that you wanted to reference Frank's Function() in Hmisc
> rather than Formula() below, to extract the regression equation.

In Frank’s rms package, there’s also a nice latex() function which
creates a human-readable LaTeX version of the regression equation (but
which doesn’t support GEE models). Example:

    mod = ols(Sepal.Length ~ rcs(Sepal.Width, 4) +
              pol(Petal.Width, 2), data = iris)

This returns:

\[{\rm E({\rm Sepal.Length}}) = X\beta, {\rm \ \ where} \\ \]
& & 2.482784 \\
& & + 0.6574976 {\rm Sepal.Width}-0.2502618({\rm
Sepal.Width}-2.345)_{+}^{3}  \\
& &  +1.315794  ({\rm Sepal.Width}-2.9)_{+}^{3}-1.366806 ({\rm
Sepal.Width}-3.2)_{+}^{3}  \\
& &   +0.301274  ({\rm Sepal.Width}-3.8)_{+}^{3}  \\
& &  +1.649973  \:{\rm Petal.Width}-0.2768114\:{\rm Petal.Width}^{2} \\
and \((x)_{+}=x\) if \(x > 0\), 0 otherwise\\

which one can paste into a LaTeX document to view (try using
https://latexbase.com/ for a quick preview).

Karl Ove Hufthammer

Reply all
Reply to author
0 new messages