set.seed(1)
fn <- function(beta, y, x, w) {
# binomial deviance using double log link
mu <- exp(x %*% beta)
logLik <- - y * mu + (w - y) * log(1 - exp(-mu))
-2 * sum(logLik)
}
n <- 1000
beta <- c(-1, 0.5)
w <- rpois(n, 100)
x <- rep(c("A", "B"), length = n)
X <- model.matrix(~ x, data.frame(x))
y <- rbinom(n, w, exp(-exp(X %*% beta)))
x1 <- powell(beta, fn, y = y, x = X, w = w)
x2 <- optim(beta, fn, y = y, x = X, w = w, hessian = TRUE)
x3 <- glm(1 - y/w ~ x, data = data.frame(x, y, w),
family = binomial("cloglog"), weights = w)
# compare coefficients from 3 fits
rbind(powell = x1$par, optim = x2$par, glm = coef(x3))
# compare standard errors from 3 fits
rbind(powell = sqrt(diag(2 * solve(x1$hessian))),
optim = sqrt(diag(2 * solve(x2$hessian))),
glm = sqrt(diag(vcov(x3))))
Run the code above in your browser using DataLab