Nothing
fit.mult <-
function(y, rates.frame, cov.frame, start)
{
if (missing(start)) {
## Fit model without covariates to get initial rates estimates
glm.out.rates <- fit.baseline(y, rates.frame)
## Initial values for iterative fitting
lambda <- coef(glm.out.rates)
beta <- rep(0, ncol(cov.frame))
}
else {
lambda <- start[1:ncol(rates.frame)]
beta <- start[ncol(rates.frame) + 1:ncol(cov.frame)]
}
niter <- 1
cy <- 1 - y
while(TRUE) {
## covariates model
off <- log(-as.matrix(rates.frame) %*% lambda)
glm.out.cov <- glm(cy ~ -1 + offset(off) + .,
family=binomial(link=cloglog),
data=cov.frame, start=beta, maxit=100)
beta <- coef(glm.out.cov)
## rates model
wgt <- exp(as.matrix(cov.frame) %*% beta)
temp.rates.frame <- wgt * rates.frame
glm.out.rates <- glm(y ~ -1 + ., family=binomial(link=log),
data=temp.rates.frame, start=lambda, maxit=100)
lambda <- coef(glm.out.rates)
## Check convergence
## Convergence <==> deviances are equal
TOL <- max(glm.out.cov$control$epsilon, glm.out.rates$control$epsilon)
dev1 <- glm.out.cov$deviance
dev2 <- glm.out.rates$deviance
if (abs(dev1 - dev2)/(0.1 + abs(dev1)) < TOL)
break
else
niter <- niter + 1
}
return(list( rates=glm.out.rates, cov=glm.out.cov, niter=niter))
}
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.