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