1fit.mult <-
2function(y, rates.frame, cov.frame, start)
3{
4  if (missing(start)) {
5    ## Fit model without covariates to get initial rates estimates
6    glm.out.rates <- fit.baseline(y, rates.frame)
7    ## Initial values for iterative fitting
8    lambda <- coef(glm.out.rates)
9    beta <- rep(0, ncol(cov.frame))
10  }
11  else {
12    lambda <- start[1:ncol(rates.frame)]
13    beta <- start[ncol(rates.frame) + 1:ncol(cov.frame)]
14  }
15
16  niter <- 1
17  cy <- 1 - y
18  while(TRUE) {
19    ## covariates model
20    off <- log(-as.matrix(rates.frame) %*% lambda)
21    glm.out.cov <- glm(cy ~ -1 + offset(off) + .,
22                       family=binomial(link=cloglog),
23                       data=cov.frame, start=beta, maxit=100)
24    beta <- coef(glm.out.cov)
25
26    ## rates model
27    wgt <- exp(as.matrix(cov.frame) %*% beta)
28    temp.rates.frame <- wgt * rates.frame
29    glm.out.rates <- glm(y ~ -1 + ., family=binomial(link=log),
30                         data=temp.rates.frame, start=lambda, maxit=100)
31    lambda <- coef(glm.out.rates)
32
33    ## Check convergence
34    ## Convergence <==> deviances are equal
35    TOL <- max(glm.out.cov$control$epsilon, glm.out.rates$control$epsilon)
36    dev1 <- glm.out.cov$deviance
37    dev2 <- glm.out.rates$deviance
38    if (abs(dev1 - dev2)/(0.1 + abs(dev1)) < TOL)
39      break
40    else
41      niter <- niter + 1
42  }
43
44  return(list( rates=glm.out.rates, cov=glm.out.cov, niter=niter))
45}
46