# R/1.2_MLE_Var.R In brm: Binary Regression Model

#### Defines functions var.mle.rdvar.mle.rr

```### variance calculation

var.mle.rr = function(x, alpha.ml, beta.ml, va, vb, weights) {

p0p1 = getProbRR(va %*% alpha.ml, vb %*% beta.ml)
n = dim(va)
pA = rep(NA, n)     # P(Y=1|A,V); here A = X
pA[x == 0] = p0p1[x == 0, 1]
pA[x == 1] = p0p1[x == 1, 2]

expect.dl.by.dpsi0.squared = (pA)/(1 - pA)
dpsi0.by.dphi = (1 - p0p1[, 1]) * (1 - p0p1[, 2])/((1 - p0p1[, 1]) + (1 -
p0p1[, 2]))
dpsi0.by.dtheta = -(1 - p0p1[, 1])/((1 - p0p1[, 1]) + (1 - p0p1[, 2]))
tmp = cbind((dpsi0.by.dtheta + x) * va, dpsi0.by.dphi * vb)
## since dtheta.by.dalpha = va, and dphi.by.dbeta = vb
fisher.info = (t(expect.dl.by.dpsi0.squared * weights * tmp) %*% tmp)
return(solve(fisher.info))
}

### variance calculation

var.mle.rd = function(x, alpha.ml, beta.ml, va, vb, weights) {

p0p1 = getProbRD(va %*% alpha.ml, vb %*% beta.ml)
# p0p1 = cbind(p0, p1): n * 2 matrix
p0 = p0p1[, 1]
p1 = p0p1[, 2]
n = nrow(va)
pA = p0             # P(Y=1|A,V); here A = X
pA[x == 1] = p1[x == 1]
s0 = p0 * (1 - p0)
s1 = p1 * (1 - p1)
sA = pA * (1 - pA)

rho = as.vector(tanh(va %*% alpha.ml))  #estimated risk differences

expect.dl.by.dpA.squared = 1/sA
dp0.by.dphi = s0 * s1/(s0 + s1)
dp0.by.drho = -s0/(s0 + s1)
drho.by.dalpha = (1 - rho^2) * va
dphi.by.dbeta = vb

tmp = cbind((dp0.by.drho + x) * drho.by.dalpha, dp0.by.dphi * dphi.by.dbeta)
fisher.info = (t(expect.dl.by.dpA.squared * weights * tmp) %*% tmp)
return(solve(fisher.info))
}
```

## 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.