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