# Extracting a linear predictor formula from GEE model

17 views

### 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")
fitar1_11
geeglm(formula = outcome ~ v1 + v2 + v4_recoded + v5_trunc +
bs(weeks, df = 3), family = binomial(), data = M2, id = f1p0,
corstr = "ar1")
Coefficients:
(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
1
0.8345324
# 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)
 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 ...

Martin

### Mitchell Maltenfort

Jun 20, 2021, 12:27:11 PM6/20/21
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 .

---
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
Hi,

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:

https://github.com/ZheyuanLi/SplinesUtils

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:

https://www.clayford.net/statistics/using-natural-splines-in-linear-modeling/

https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-019-0666-3

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:

https://github.com/datalorax/equatiomatic

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.

Regards,

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
'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:

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

This returns:

${\rm E({\rm Sepal.Length}}) = X\beta, {\rm \ \ where} \\$
\begin{eqnarray*}
\lefteqn{X\hat{\beta}=}\\
& & 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} \\
\end{eqnarray*}
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