# R/2.2_DR_Var_RR.R In brm: Binary Regression Model

#### Defines functions var.rr.dr

```## Sandwich estimator for variance of RR

var.rr.dr = function(y, x, va, vb, vc, alpha.dr, alpha.ml, beta.ml, gamma,
optimal, weights) {
########################################

pscore = as.vector(expit(vc %*% gamma))
n = length(pscore)

### 1. - E[dS/d(alpha.ml,beta.ml)] ############################## Computing
### the Hessian:

Hrr = Hessian2RR(y, x, va, vb, alpha.ml, beta.ml, weights)
hessian = Hrr\$hessian
p0 = Hrr\$p0; p1  = Hrr\$p1; pA = Hrr\$pA
dpsi0.by.dtheta  = Hrr\$dpsi0.by.dtheta
dpsi0.by.dphi    = Hrr\$dpsi0.by.dphi
dtheta.by.dalpha = Hrr\$dtheta.by.dalpha
dphi.by.dbeta    = Hrr\$dphi.by.dbeta
dl.by.dpsi0      = Hrr\$dl.by.dpsi0

############# extra building blocks ##########################

H.alpha = y * exp(-x * (as.vector(va %*% alpha.dr)))

############# Calculation of optimal vector (used in several places below) ##

if (optimal == TRUE) {
theta = as.vector(va %*% alpha.ml)  # avoid n by 1 matrix
dtheta.by.dalpha.beta = cbind(va, matrix(0, n, length(beta.ml)))
wt = 1/(1 - p0 + (1 - pscore) * (exp(-theta) - 1))
} else {
wt = rep(1, n)
}

### 2. -E[dU.by.dalphaml.betaml] ####################################

dU.by.dp0 = -va * wt * (x - pscore)  # n by 2
dp0.by.dpsi0 = p0
dpsi0.by.dalpha.beta = cbind(dpsi0.by.dtheta * dtheta.by.dalpha, dpsi0.by.dphi *
dphi.by.dbeta)  # n by 4
# 4 = 2 (alpha) + 2 (beta)
dp0.by.dalpha.beta = dpsi0.by.dalpha.beta * dp0.by.dpsi0  # n by 4

dU.by.dwt = va * (x - pscore) * (H.alpha - p0)  # n by 2
dwt.by.dwti = -wt^2  # n
# wti is short for wt_inv
dU.by.dwti = dU.by.dwt * dwt.by.dwti  # n by 2
if (optimal == TRUE) {
dwti.by.dalpha.beta = -dp0.by.dalpha.beta - (1 - pscore) * exp(-theta) *
dtheta.by.dalpha.beta  # n by 4
} else {
dwti.by.dalpha.beta = matrix(0, n, ncol(va) + ncol(vb))
}

dU.by.dalpha.ml.beta.ml = t(dU.by.dp0 * weights) %*% dp0.by.dalpha.beta +
t(dU.by.dwti * weights) %*% dwti.by.dalpha.beta

### 3. tau = -E[dU/dalpha.dr] ######################################## (This
### is the bread of the sandwich estimate)

dU.by.dH = va * wt * (x - pscore)  # n by 2
dH.by.dalpha.dr = -va * x * H.alpha  # n by 2

tau = -t(dU.by.dH * weights) %*% dH.by.dalpha.dr/sum(weights)  # 2 by 2

### 4. E[d(prop score score equation)/dgamma]

dpscore.by.dgamma = vc * pscore * (1 - pscore)  # n by 2
part4 = -t(vc * weights) %*% dpscore.by.dgamma  # 2 by 2

### 5. E[dU/dgamma]

dU.by.dpscore = -va * wt * (H.alpha - p0)  # n by 2

if (optimal == TRUE) {
dwti.by.dpscore = 1 - exp(-theta)  # n
dwti.by.dgamma = dpscore.by.dgamma * dwti.by.dpscore  # n by 2
} else {
dwti.by.dgamma = matrix(0, n, ncol(vc))
}

dU.by.dgamma = t(dU.by.dpscore * weights) %*% dpscore.by.dgamma + t(dU.by.dwti *
weights) %*% dwti.by.dgamma  # 2 by 2

############################################################################# Assembling semi-parametric variance matrix

U = va * wt * (x - pscore) * (H.alpha - p0)  # n by 2

S = cbind(dl.by.dpsi0 * (dpsi0.by.dtheta + x) * dtheta.by.dalpha, dl.by.dpsi0 *
dpsi0.by.dphi * dphi.by.dbeta)
pscore.score = vc * (x - pscore)

Utilde = U - t(dU.by.dalpha.ml.beta.ml %*% (-solve(hessian)) %*% t(S)) -
t(dU.by.dgamma %*% (solve(part4)) %*% t(pscore.score))  # n by 2
USigma = t(Utilde * weights) %*% Utilde/sum(weights)

################################### Asymptotic var matrix for alpha.dr

alpha.dr.variance = solve(tau) %*% USigma %*% solve(tau)/sum(weights)

return(alpha.dr.variance)

}
```

## Try the brm package in your browser

Any scripts or data that you put into this service are public.

brm documentation built on July 1, 2020, 10:35 p.m.