# Estimating risk differences for marginal models

I need to estimate risk differences for marginal models. I have worked on the Probratio function to create an equivalent function for RDs that works.

Please provide feedback on accuracy or comments to improve efficiency where necessary.

 probdiff <- function(object, parm, subset, method=c('ML', 'delta', 'bootstrap'), scale='linear', level=0.95, seed, NREPS=1000, ...) {
if (length(match('glm', class(object))) < 0)
stop('Non GLM input to "object"')
if (family(object)$family != 'binomial') stop('object not a logistic regression model') nc <- length(cf <- coef(object)) if (missing(subset)) subset <- T if (missing(parm)) parm <- seq(2, nc) cf <- cf[parm] method <- match.arg(method, c('ML', 'delta', 'bootstrap')) scale <- 'linear' if (is.na(scale)) stop('scale cannot take values outside of linear') if (scale == 'linear') { f <- function(x) x[2] - x[1] name <- c('Risk difference') null <- 0 } else{ (is.na(scale)) stop('scale cannot take values outside of linear') } cilevel <- c({1-level}/2, 1-{1-level}/2) ciname <- paste0(c('Lower', 'Upper'), ' ', formatC(100*cilevel, format='f', digits=1), '% CI') if (method == 'ML') { newfit <- glm(object, family=binomial(link="identity"), subset=subset, ...) out <- coef(summary(newfit))[parm, , drop=F] if (scale == 'linear') { out <- cbind(out[, 1], out[, 2], out[, 3], out[, 4], out[, 1] + qnorm((1-level)/2)*out[, 2], out[, 1] + qnorm(1-(1-level)/2)*out[, 2]) } else { newfit <- glm(object, family=gaussian(link = "identity"), subset=subset, ...) if (family(object)$family != 'binomial')
out <- coef(summary(newfit))[parm, , drop=F]
out = lmtest::coeftest(out, vcov.=sandwich::vcovHC(out, type="HC3"))
val <- out[, 1]
se <- val * out[, 2]
out <- cbind(val, se, z <- abs(val-1)/se, pnorm(z, lower.tail = F)*2,
val + qnorm((1-level)/2)*se, val + qnorm(1-(1-level)/2)*se)
}
colnames(out) <- c(name, 'Std. Error', 'Z-value', 'p-value', ciname)
return(out)
}## endif MLe code here

Mod1 <- Mod0 <- model.matrix(object)[subset, ]
n <- nrow(Mod0)
Nvec <- matrix(rep(c(1/n,0,0,1/n),each=n), n*2, 2)

if (method == 'delta') {
if (scale == 'linear') {
df <- deriv( ~y/x, c('x', 'y')) # y depends on x
}

out <- sapply(parm, function(p) {
Mod0[, p] <- 0
Mod1[, p] <- 1
Mod <- rbind(Mod0, Mod1)
allpreds <- family(object)$linkinv(Mod %*% coef(object)) avgpreds <- t(Nvec) %*% allpreds val <- f(avgpreds) V <- sweep(chol(vcov(object)) %*% t(Mod), allpreds*(1-allpreds), '*', MARGIN = 2) %*% Nvec V <- t(V) %*% V dxdy <- matrix(attr(eval(df, list('x'=avgpreds[1], 'y'=avgpreds[2])), 'gradient')) se <- sqrt(t(dxdy) %*% V %*% dxdy) out <- c(val, se, z <- abs({val-null}/se), 2*pnorm(abs(val/se), lower.tail=FALSE), val + qnorm(cilevel[1])*se, val + qnorm(cilevel[2])*se) names(out) <- c(name, 'Std. Error', 'Z-value', 'p-value', ciname) out }) # / sqrt( n ) out <- t(out) rownames(out) <- names(cf) return(out) } ## endif deltaMod1 <- Mod0 <- model.matrix(object)[subset, ] n <- nrow(Mod0) Nvec <- matrix(rep(c(1/n,0,0,1/n),each=n), n*2, 2) if (method == 'delta') { if (scale == 'linear') { df <- deriv( ~y/x, c('x', 'y')) # y depends on x } out <- sapply(parm, function(p) { Mod0[, p] <- 0 Mod1[, p] <- 1 Mod <- rbind(Mod0, Mod1) allpreds <- family(object)$linkinv(Mod %*% coef(object))
avgpreds <- t(Nvec) %*% allpreds
val <- f(avgpreds)
V <- sweep(chol(vcov(object)) %*% t(Mod), allpreds*(1-allpreds), '*', MARGIN = 2) %*% Nvec
V <- t(V) %*% V
dxdy <- matrix(attr(eval(df, list('x'=avgpreds[1], 'y'=avgpreds[2])), 'gradient'))
se <- sqrt(t(dxdy) %*% V %*% dxdy)
out <- c(val, se, z <- abs({val-null}/se), 2*pnorm(abs(val/se), lower.tail=FALSE), val + qnorm(cilevel[1])*se, val + qnorm(cilevel[2])*se)
names(out) <- c(name, 'Std. Error', 'Z-value', 'p-value', ciname)
out
}) # / sqrt( n )
out <- t(out)
rownames(out) <- names(cf)
return(out)
} ## endif delta code here

if (method == 'bootstrap') {
if (missing(seed))
stop('seed must be supplied by the user when obtaining results from random number generation')
set.seed(seed)
out <- replicate(NREPS, {

index <- sample(1:n, n, replace=T)
Mod <- model.matrix(object)[subset, ][index, ]
newbeta <- glm.fit(Mod, object$y[index], family=binomial())$coef
out <- sapply(parm, function(p) {
Mod1 <- Mod0 <- Mod
Mod1[, p] <- 1
Mod0[, p] <- 0
Mod <- rbind(Mod0, Mod1)
f(t(Nvec) %*% newpreds)
})
out
})

if (length(parm) == 1) {
out <- c(val <- mean(out), se <- sd(out), z <- abs({val - null}/se), 2*pnorm(abs(val/se), lower.tail=FALSE),
val + qnorm((1-level)/2)*se, val + qnorm(1-(1-level)/2)*se)
names(out) <- c(name, 'Std. Error', 'Z-value', 'p-value', ciname)
} else {
out <- cbind(val <- rowMeans(out), se <- apply(out, 1, sd), z <- abs({val - null}/se), 2*pnorm(abs(val/se), lower.tail=FALSE),
val + qnorm(cilevel[1])*se, val + qnorm(cilevel[2])*se)
colnames(out) <- c(name, 'Std. Error', 'Z-value', 'p-value', ciname)
rownames(out) <- names(cf)
}

return(out)
}
}


• Welcome to Code Review! I changed the title so that it describes what the code does per site goals: "State what your code does in your title, not your main concerns about it.". Feel free to edit and give it a different title if there is something more appropriate. Nov 1 at 17:06