Inverse Probability Tilting (IPT) Weights from Graham et al. (2012)

IPT.R

# IPT_ATT() and IPT_ATE() compute IPT weights for the ATT or ATE, respectively.
# `formula` should be a model formula for the treatment, and `data` should be the
# dataset (as in `glm()`). Code for the ATT weights comes form the `drdid` package.

IPT_ATT <- function(formula, data) {
    
    mf <- model.frame(formula, data = data)
    
    int.cov <- model.matrix(formula, data = mf)
    int.cov[,-1] <- scale(int.cov[,-1])
    D <- model.response(mf)
    
    pslogit <- glm.fit(x = int.cov, y = D, family = binomial())
    
    init.gamma <- pslogit$coefficients
    
    loss.ps.IPT <- function(gamma1, n, D, int.cov){
        #Coefficients for quadratic extrapolation
        cn <- -(n - 1)
        bn <- -n + (n - 1) * log(n - 1)
        an <- -(n - 1) * (1 - log(n - 1) + 0.5 * (log(n - 1))^2)
        vstar <- log(n - 1)
        
        v <- gamma1 %*% t(int.cov)
        phi <- ifelse(v < vstar, - v - exp(v), an + bn * v + 0.5 * cn * (v^2))
        phi1 <- ifelse(v < vstar, - 1 - exp(v), bn + cn * v)
        phi2 <- ifelse(v < vstar, - exp(v), cn)
        
        # Minus is because nlm minimizes functions, and we aim to maximize!
        res <- - sum((1 - D) * phi + v)
        
        attr(res, "gradient") <- - t(int.cov) %*% as.vector(((1-D) * phi1 + 1))
        attr(res, "hessian")  <-  - t(as.vector((1-D) *phi2) * int.cov) %*% int.cov
        return(res)
    }
    
    pscore.IPT <- suppressWarnings(stats::nlm(loss.ps.IPT, init.gamma, iterlim = 10000, gradtol = 1e-06,
                                              check.analyticals = F,
                                              D = D, int.cov = int.cov,
                                              n = length(D)))
    gamma.cal <- try(pscore.IPT$estimate)
    
    #Compute fitted pscore and weights for regression
    pscore.index <- tcrossprod(gamma.cal, int.cov)
    p <- as.numeric(stats::plogis(pscore.index))
    
    ifelse(D == 1, 1, p/1-p)
}

IPT_ATE <- function(formula, data) {
    
    mf <- model.frame(formula, data = data)
    
    int.cov <- model.matrix(formula, data = mf)
    int.cov[,-1] <- scale(int.cov[,-1])
    D <- model.response(mf)
    
    pslogit <- glm.fit(x = int.cov, y = D, family = binomial())
    
    init.gamma <- pslogit$coefficients

    f <- function(b, x, D, t) {
        do.call("cbind", lapply(1:ncol(x), function(i) {
            if (t == 1)
                (D / plogis(x %*% b) - 1) * x[,i]
            else
                ((1-D) / (1-plogis(x %*% b)) - 1) * x[,i]
        })) |>
            colMeans()
    }
    
    sol1 <- rootSolve::multiroot(f, x = int.cov, D = D, t = 1, start = init.gamma)
    sol0 <- rootSolve::multiroot(f, x = int.cov, D = D, t = 0, start = init.gamma)
    
    p <- ifelse(D == 1,
                drop(plogis(int.cov %*% sol1$root)),
                drop(plogis(int.cov %*% sol0$root)))
    
    ifelse(D == 1, 1/p, 1/1-p)
}