Nothing
make.control <- function(max.retries = 10, method = "BFGS", ...) {
if (max.retries <= 0) stop("maximum number of retries must be > 0")
list(max.retries = max.retries, method = method, ...)
}
fit.algorithm <- function(par, f, g, Y, X, P, setM,
control, silent = FALSE) {
method <- control$method
max.retries <- control$max.retries
control$method <- NULL
control$max.retries <- NULL
retry <- TRUE
retry_count <- 0
while (retry) {
if (method == "BFGS2") {
erg <- ucminf(par = par, fn = f, gr = g, Y = Y, X = X, W = P, SetM = setM,
control = control, hessian = TRUE)
} else {
erg <- optim(par = par, fn = f, gr = g, Y = Y, X = X, W = P, SetM = setM,
method = method, control = control, hessian = TRUE)
}
# perform check
erg$hessian <- (erg$hessian + t(erg$hessian)) / 2
retry <- inherits(try(chol(erg$hessian)), "try-error")
if (retry) {
if (!silent) cat("wrong hessian...resetting parameter ")
#PD <- Matrix::nearPD(erg$hessian)
#print(PD$normF)
#stop()
retry_count <- retry_count + 1
if (retry_count > max.retries) stop("too many retries... hessian is not positive definite!")
#print(erg$par)
i <- which.max(abs(erg$par))
if (!silent) cat(i)
if (!silent) cat(" (")
if (!silent) cat(erg$par[i])
if (!silent) cat(")...")
par <- erg$par
par[i] <- 0
}
}
return(erg)
}
bootstrapping <- function(ret, Y, X, Pmodel, PX, f, g,
boot.fraction = 1, repetitions = 1000) {
if (is.matrix(Y)) {
k <- ncol(Y) - 1
if (k == 0) k <- 1
} else {
k <- 1
}
simuliere.predict <- function(Pmodel, newdata, u, prep, s) {
Pmodel$coefficients <- prep
P <- predict(Pmodel, newdata = newdata[s,], type = "response")
if (u == 2) {
P0 <- 1 - P
P <- cbind(P0, P)
dimnames(P)[[2]] <- c("M0", "M1")
}
return(P)
}
u <- nrow(ret$setM)
n <- nrow(X)
m <- floor(boot.fraction * n)
par <- ret$optim$par
boot.beta <- matrix(0.0, nrow = repetitions, ncol = length(par))
cat("bootstrapping")
s <- matrix(as.integer(0), nrow = m, ncol = repetitions)
prep <- matrix(0.0, nrow = repetitions,
ncol = length(Pmodel$coefficients))
for (i in 1:repetitions) {
# ziehen
s[, i] <- sample.int(n, m, replace = TRUE)
prep[i,] <- mvrnorm(n = 1, mu = coef(Pmodel), Sigma = vcov(Pmodel))
}
`%op%` <- ifelse(getDoParWorkers() > 1, `%dopar%`, `%do%`)
tmp <- foreach(i = 1:repetitions, .inorder = FALSE) %op% {
si <- s[,i]
# P simulieren
P <- simuliere.predict(Pmodel, newdata = PX, u, prep[i,, drop = TRUE], si)
# single step
if (k > 1) {
tmp2 <- fit.algorithm(par, f, g, Y[si,, drop = FALSE], X[si,, drop = FALSE], P, ret$setM,
ret$control, silent = TRUE)$par
} else {
tmp2 <- fit.algorithm(par, f, g, Y[si], X[si,, drop = FALSE], P, ret$setM,
ret$control, silent = TRUE)$par
}
cat(".")
rm(si)
rm(P)
gc()
tmp2
}
for (i in 1:repetitions) {
boot.beta[i,] <- tmp[[i]]
}
rm(s)
rm(prep)
cat("done.\n")
ret$boot.beta <- boot.beta
ret$CoVar <- m / n * cov(boot.beta)
ret$boot.fraction <- boot.fraction
ret$SEtype <- "bootstrap"
return(ret)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.