Predict probabilities from Firth logistic regression in R

Question

I'm using the logistf package in R to perform Firth logistic regression on an unbalanced dataset. I have a logistf object:

fit = logistf(a~b)

Is there a predict() function like on that's used in the lm class to predict probabilities for future data points? Or do I have to manually input the estimated parameters from the Firth regression.

Answer 1

You can probably compute any predictions you want with little algebra. Let consider the example dataset,

data(sex2)
fm <- case ~ age+oc+vic+vicl+vis+dia
fit <- logistf(fm, data=sex2)

A design matrix is the only missing piece to compute predicted probabilities once we get the regression coefficients, given by

betas <- coef(fit)

So, let's try to get prediction for the observed data, first:

X <- model.matrix(fm, data=sex2)       # add a column of 1's to sex2[,-1]
pi.obs <- 1 / (1 + exp(-X %*% betas))  # in case there's an offset, δ, it 
                                       # should be subtracted as exp(-Xβ - δ)

We can check that we get the correct result

> pi.obs[1:5]
[1] 0.3389307 0.9159945 0.9159945 0.9159945 0.9159945
> fit$predict[1:5]
[1] 0.3389307 0.9159945 0.9159945 0.9159945 0.9159945

Now, you can put in the above design matrix, X, values you are interested in. For example, with all covariates set to one

new.x <- c(1, rep(1, 6))
1 / (1 + exp(-new.x %*% betas)) 

we get an individual probability of 0.804, while when all covariates are set to 0 (new.x <- c(1, rep(0, 6))), the estimated probability is 0.530.

Answer 2

An alternative approach is the brglm package. For example, using the same data/model as @chl's Answer

data(sex2)
fm <- case ~ age + oc + vic + vicl + vis + dia
fit <- brglm(fm, data = sex2)

predict(fit, newdata = sex2[1:5, ], type = "response")

That yields the

> predict(fit, newdata = sex2[1:5, ], type = "response")
        1         2         3         4         5 
0.3389307 0.9159945 0.9159945 0.9159945 0.9159945

Note that in the brglm() case, because of the way the function works, what you see above is simply the result of the standard predict.glm() function/method in R.