Please note that this page is incomplete and there may be inconsistencies in the code or explanations. Additionally, this page will not be updated in the near future.
Adjusted predictions are often calculated to predict the response at a given set of predictor values, usually to get an idea of the response value at representative predictor values. (NOT SURE ABOUT THIS) Marginal means are often calculated to isolate the effect of a single predictor by averaging the predicted response over various levels of the predictor of interest while holding the values of other predictors constant. Both conditional and marginal means are functions of the model coefficients.
We can use the delta method to get the variance of a function of random variable. As model coefficients are themselves random variables, we can use the delta method to get the variance of conditional and marginal means, because they are functions of the model ceofficients. For a function G(b) of a random variable vector, b, the variance of G(b) by the delta method is:
var(G(B)) = dG/db’ %*% var(b) %*% dG/db,where dG/db is a vector of partial derivatives of G(b) with respect to
each b coefficient, also known as the Jacobian matrix, and var(b) is the variance-covariance matrix of the
b coefficients. The predict()
function calculates delta-method standard errors for conditional means, but it will not quite work for marginal means.
Example 1: Delta method standard error for conditional mean of Y at mean of X
First let’s make up some data and run a very simple linear regression.
x <- 1:10 y <- c(1,3,3,4,5,7,7,8,9,10) m1 <- lm(y~x)
Let’s get the conditional mean first and delta standard error returned by predict()
.
predict(m1, newdata=(x=5.5), se=T) $fit 1 5.7 $se.fit [1] 0.1365151 $df [1] 8 $residual.scale [1] 0.4316985
Let’s see if we can replicate this standard error manually. The two components we need are the variance-covariance matrix of the coefficients, and the vector of partial derivatives. First, get the variance-covariance matrix of b0 (the intercept) and b1 (the coefficient for x) using vcov()
.
vb <- vcov(m1)
For linear models, the transformation from model coefficients to conditional means is simple: G(b) = b0 + b1*X. We want standard error of G(b), the conditonal mean, at the mean of x, x=5.5. So the transformation equation is G(b) = b0*1 + b1*5.5. The partial derivatives with respect to each coefficient are dG/db0=1 and dG/db2=5.5. The vector of partial derivatives of G(b) with respect to b is thus: [1 5.5]. One can use the deriv() function to get the vector of partial derivatives, but for very simple transformation functions such as these, it is easier to just enter it manually.
*************FIX THIS UP
> dgdb <- deriv(y ~ b0 + 5.5*b1, c(“b0”, “b1”)) > b0 <- 1 > b1 <- 1 > eval(dgdb)
***********ALMOST THERE
dgdb <- c(1,5.5)
Now we apply the delta method: var(G(B)) = dG/db’ %*% var(b) %*% dG/db. We then take the sqrt of the variance to get the standard error of the conditional mean.
vg <- dgdb %*% vb %*% dgdb delta <- sqrt(vg) delta
####EXAMPLE 2 ################################################# # Delta method standard error for predicted value of Y at set values of # several Xs d <- read.csv(“https://stats.idre.ucla.edu/stat/data/hsb2.csv”) d$prog <- factor(d$prog) d$female <- factor(d$female) m2 <- lm(write ~ female + prog + math, data=d)
#get var-covar matrix of all bs vb <- vcov(m2) vb
# G(b) = b0 + b_female1*female1 + b_prog2*prog2 + b_prog3*prog3 + b_math*math # here we want standard error of G(b) for females in prog2 with math = 50 # regression equation is Y = b0*1 + b_female*1 + b_prog2*1 + b_prog3*0 +b_math*50 # so vector of partial derivatives is 1,1,1,0,50 dgdb <- c(1,1,1,0,50) #variance of G(b) vg <- dgdb %*% vb %*% dgdb #sqrt of variance for delta method standard error delta <- sqrt(vg) delta
####EXAMPLE 3 ################################################# # Delta method standard error for predicted value of Y at set values of # several Xs IN LOGISTIC MODEL d <- read.csv(“https://stats.idre.ucla.edu/stat/data/hsb2.csv”) d$prog <- factor(d$prog) d$female <- factor(d$female) m3 <- glm(female ~ prog + math + read, data=d, family=binomial)
#get var-covar matrix of all bs vb <- vcov(m3) vb
# FOR STANDARD ERROR OF XB # G(b) = b0 + b_prog2*prog2 + b_prog3*prog3 + b_math*math + b_read*read # here we want standard error of G(b) for males in prog1 with math = 50 and read = 50 # regression equation is Y = b0*1 + b_prog2*0 + b_prog3*0 +b_math*50 + b_read*50 # so vector of partial derivatives is 1,0,0,50,50 dgdb <- c(1,0,0,50,50) #variance of G(b) vg <- dgdb %*% vb %*% dgdb #sqrt of variance for delta method standard error delta <- sqrt(vg) delta
#FOR STANDARD ERROR ON RESPONSE (PROBABILITY) METRIC #The key here is to apply the delta method twice #First to get the standard error on the linear metric #Second to get the standard error on response metric #because response is function of linear predictor #In logistic regression, to transform linear predictor L to #probability P, the function is: # P = 1/(1+exp(-L))
#Using the chain rule, we see that dP/dL = 1/(1+exp(-L))*exp(-L)/(1+exp(-L)), #so dP/dL = P * (1-P) #Work out derivative: #Let P = (1 + g)^-1, g = exp(h), and h = -1. #Therefore: dP/dg = -(1 + g)^2 # dg/dh = exp(h) # dh/dL = -x #Let dP/dL = dP/dg*dg/dh*dh/dL #so dP/dL = -exp(-x)/-(1+exp(-x))^2 = 1/(1+exp(-L))*exp(-L)/(1+exp(-L))
d <- read.csv(“https://stats.idre.ucla.edu/stat/data/hsb2.csv”) d$prog <- factor(d$prog) d$female <- factor(d$female) m3 <- glm(female ~ prog + math + read, data=d, family=binomial) vb <- vcov(m3) vb #first delta is standard error on log odds scale delta <- sqrt(vg) delta #We also need to get P itself newd <- d newd$female <- 0 newd$prog <- as.factor(1) newd$math <- 50 newd$read <- 50 newd L <- mean(predict(m3, newdata=newd)) P <- 1/(1+exp(-L))
dpdl = P*(1-P) vp = dpdl * delta^2 *dpdl deltap = sqrt(vp) deltap