1"float" <-
2function(object, factor, iter.max = 50)
3{
4  float.variance <- function(V, tol=1.0e-3, iter.max = 50)
5    {
6      ## Calculate floated variances for variance matrix V, which is
7      ## assumed to represent a set of treatment contrasts
8
9      m <- nrow(V)
10      if (!is.matrix(V) || ncol(V) != m || m == 1)
11        stop ("V must be a square matrix of size 2 x 2 or more")
12
13      evals <- eigen(V, only.values=TRUE)$values
14      if(any(evals < 0))
15        stop("V not positive definite")
16
17      ## Starting values from Easton et al (1991)
18      R <- V - diag(diag(V))
19      V00 <- sum(R)/(m * (m-1))
20      V10 <- apply(R, 1, sum)/(m-1)
21      fv <- c(V00, V00 - 2 * V10 + diag(V))
22
23      for(iter in 1:iter.max) {
24        w <- 1/fv
25        S <- sum(w)
26        w1 <- w[-1]/S
27        ##Augment data matrix
28        V10 <- as.vector(V %*% w1)
29        V00 <- as.vector(1/S + t(w1) %*% V %*% w1)
30        ##Calculate new estimates
31        fv.old <- fv
32        fv <- c(V00, V00 - 2 * V10 + diag(V))
33        ## Check convergence
34        if(max(abs(fv.old - fv)/fv) < tol)
35          break
36      }
37      if (iter == iter.max)
38        warning("Floated variance estimates did not converge")
39
40      Vmodel.inv <- S * (diag(w1) - w1 %*% t(w1))
41      evals <- 1/(eigen(V %*% Vmodel.inv, only.values=TRUE)$values)
42      divergence <- sum(1/evals - 1 + log(evals))/2
43      return(list(variance=fv, error.limits=sqrt(range(evals)),
44                  divergence=divergence))
45    }
46
47  if (is.null(object$xlevels)) {
48    stop("No factors in model")
49  }
50
51  if (missing(factor)) {
52    i <- 1
53    factor <- names(object$xlevels)[1]
54  }
55  else {
56    i <- pmatch(factor, names(object$xlevels))
57    if (is.na(i)) {
58      stop(paste("Factor",i,"not found in model"))
59    }
60  }
61
62  xcontrasts <- object$contrasts[[i]]
63  xlevels <- object$xlevels[[i]]
64  xname <- names(object$xlevels)[i]
65  nlevels <- length(xlevels)
66
67  ## Extract the coefficients and variance matrix for a single factor
68  ## from object
69
70  if (nlevels <= 2) {
71    stop ("Floated variances undefined for factors with less than 3 levels")
72  }
73
74  ## Get contrast matrix
75  C <- if (is.matrix(xcontrasts)) {
76    xcontrasts
77  }
78  else {
79    get(xcontrasts, mode="function")(xlevels)
80  }
81  if (qr(C)$rank < nlevels - 1) {
82    stop ("Impossible to reconstruct treatment contrasts")
83  }
84
85  ## Get coefficients and variance matrix
86  if(is.null(cnames <- colnames(C)))
87    cnames <- 1:(nlevels-1)
88  contr.names <- paste(xname, cnames, sep="")
89  coef <- coef(object)[contr.names]
90  V <- vcov(object)[contr.names, contr.names]
91
92  ## Convert to treatment contrast parameterization
93  if (identical(xcontrasts, "contr.treatment")) {
94    V.tc <- V
95    coef.tc <- c(0, coef)
96  }
97  else {
98    D.inv <- cbind(rep(-1,nlevels-1), diag(nlevels-1))
99    S <- D.inv %*% cbind(rep(1, nlevels), C)
100    S <- S[,-1]
101    ## coefficients
102    coef.tc <- c(0, S %*% coef)
103    ## If we find a baseline level (implicitly defined
104    ## by having a row of zeros in the contrast matrix)
105    ## then adjust the coefficients
106    is.base <- apply(abs(C), 1, sum) == 0
107    if (any(is.base))
108      coef.tc <- coef.tc - coef.tc[is.base]
109    ## variance matrix
110    V.tc <- S %*% V %*% t(S)
111  }
112  names(coef.tc) <- xlevels
113
114  float.out <- float.variance(V.tc, iter.max = iter.max)
115  var <- float.out$var
116  names(var) <- xlevels
117  ans <- list(coef=coef.tc, var=var, limits=float.out$error.limits,
118              factor=factor)
119  class(ans) <- "floated"
120  return(ans)
121}
122