@andrew I would recommend to build the model matrix with the formula/terms stored in the model itself. This makes sure it is in sync with the model coefficients.
x <- model |>
terms() |>
delete.response() |>
model.matrix(data = data_to_plug_in)
Then you can compute the linear predictor as before:
lin_pred <- drop(x %*% coef(model))
where drop() preserves the row names.
Finally, the inverse link is also stored in the model object: