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.

Thanks!

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.