Extracting a linear predictor formula from GEE model

[Pedro Emmanuel Alvarenga Americano do Brasil]

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)
[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]

Actually, Stata is more intuitive ...

Martin

[Mitchell Maltenfort]

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]

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]

'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