predict.logit: Predicted Probability with Influence Function

In BuyseTest: Generalized Pairwise Comparisons

Predicted Probability with Influence Function

Description

Compute the predicted probabilities from a logistic regression, with their (robust) standard error, and the corresponding influence function.

Usage

.predict.logit(object, newdata, level = NULL, robust = TRUE)

Arguments

object	a logistic model.
newdata	[data.frame] dataset containing
level	[character] level of the outcome for which the probability should be computed.
robust	[logit] when FALSE uses the individual contribution to the modeled variance-covariance matrix as iid decomposition.

Examples

## Not run:  ## will not run as .predict.logit is not exported
set.seed(10)
n <- 100
df <- data.frame(Y = rbinom(n, prob = 0.5, size = 1), X1 = rnorm(n), X2 = rnorm(n))
e.logit <- glm(Y~X1+X2, data = df, family = binomial(link="logit"))

test1 <- .predict.logit(e.logit, newdata = df[1:5,])
test1["se",] - sqrt(diag(crossprod(attr(test1,"iid")))) ## exact
test1["var.se",] - diag(crossprod(attr(test1,"iid.se"))) ## exact

GS <- predict(e.logit, newdata = df[1:5,], se = TRUE, type =  "response")
test2 <- .predict.logit(e.logit, newdata = df[1:5,], robust = FALSE)
test2["estimate",] - GS$fit ## exact
test2["se",] - GS$se.fit ## exact
test2["se",] - sqrt(diag(crossprod(attr(test2,"iid")))) ## approximate
## since it uses the robust estimator of the correlation
## (and the modeled estimator of the variance)

## Sanity check (fully stratified model)
df <- data.frame(Y = rbinom(n, prob = 0.5, size = 1),
                 X1 = rnorm(n),
                 X2 = as.factor(rbinom(n, size = 1, prob = 0.5)))
newdata <- data.frame(X1=c(0,1),X2=as.factor(0:1))

e.logit <- glm(Y~X1+X2, data = df, family = binomial(link="logit"))
e.predlogit <- .predict.logit(e.logit, newdata = newdata)
cor(attr(e.predlogit,"iid"))

e.logitS <- glm(Y~X1*X2, data = df, family = binomial(link="logit"))
e.predlogitS <- .predict.logit(e.logitS, newdata = newdata)
cor(attr(e.predlogitS,"iid"))

## End(Not run)

BuyseTest documentation built on March 31, 2023, 6:55 p.m.