1## $Id: estimable.R 2060 2015-07-19 03:22:30Z warnes $ 2estimable <- function (obj, cm, beta0, conf.int=NULL, show.beta0, ...) 3 { 4 UseMethod("estimable") 5 } 6 7estimable.default <- function (obj, cm, beta0, conf.int=NULL, 8 show.beta0, joint.test=FALSE, ...) 9{ 10 if (is.matrix(cm) || is.data.frame(cm)) 11 { 12 cm <- t(apply(cm, 1, .to.est, obj=obj)) 13 } 14 else if(is.list(cm)) 15 { 16 cm <- matrix(.to.est(obj, cm), nrow=1) 17 } 18 else if(is.vector(cm)) 19 { 20 cm <- matrix(.to.est(obj, cm), nrow=1) 21 } 22 else 23 { 24 stop("`cm' argument must be of type vector, list, or matrix.") 25 } 26 27 if(missing(show.beta0)) 28 { 29 if(!missing(beta0)) 30 show.beta0=TRUE 31 else 32 show.beta0=FALSE 33 } 34 35 36 if (missing(beta0)) 37 { 38 beta0 = rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm))) 39 40 } 41 42 43 if (joint.test==TRUE) 44 { 45 .wald(obj, cm, beta0) 46 } 47 else 48 { 49 if ("lme" %in% class(obj)) { 50 stat.name <- "t.stat" 51 cf <- summary(obj)$tTable 52 rho <- summary(obj)$cor 53 vcv <- rho * outer(cf[, 2], cf[, 2]) 54 tmp <- cm 55 tmp[tmp==0] <- NA 56 df.all <- t(abs(t(tmp) * obj$fixDF$X)) 57 df <- apply(df.all, 1, min, na.rm=TRUE) 58 problem <- apply(df.all !=df, 1, any, na.rm=TRUE) 59 if (any(problem)) 60 warning(paste("Degrees of freedom vary among parameters used to ", 61 "construct linear contrast(s): ", 62 paste((1:nrow(tmp))[problem], collapse=", "), 63 ". Using the smallest df among the set of parameters.", 64 sep="")) 65 } 66 else if ("lm" %in% class(obj)) 67 { 68 stat.name <- "t.stat" 69 cf <- summary.lm(obj)$coefficients 70 vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2 71 df <- obj$df.residual 72 if ("glm" %in% class(obj)) 73 { 74 vcv <- summary(obj)$cov.scaled 75 if(family(obj)[1] %in% c("poisson", "binomial")) 76 { 77 stat.name <- "X2.stat" 78 df <- 1 79 } 80 else 81 { 82 stat.name <- "t.stat" 83 df <- obj$df.residual 84 } 85 } 86 } 87 else if ("geese" %in% class(obj)) 88 { 89 stat.name <- "X2.stat" 90 cf <- summary(obj)$mean 91 vcv <- obj$vbeta 92 df <- 1 93 } 94 else if ("gee" %in% class(obj)) 95 { 96 stat.name <- "X2.stat" 97 cf <- summary(obj)$coef 98 vcv <- obj$robust.variance 99 df <- 1 100 } 101 else 102 { 103 stop("obj must be of class 'lm', 'glm', 'aov', 'lme', 'gee', 'geese' or 'nlme'") 104 } 105 if (is.null(cm)) 106 cm <- diag(dim(cf)[1]) 107 if (!dim(cm)[2]==dim(cf)[1]) 108 stop(paste("\n Dimension of ", deparse(substitute(cm)), 109 ": ", paste(dim(cm), collapse="x"), 110 ", not compatible with no of parameters in ", 111 deparse(substitute(obj)), ": ", dim(cf)[1], sep="")) 112 ct <- cm %*% cf[, 1] 113 ct.diff <- cm %*% cf[, 1] - beta0 114 115 vc <- sqrt(diag(cm %*% vcv %*% t(cm))) 116 if (is.null(rownames(cm))) 117 rn <- paste("(", apply(cm, 1, paste, collapse=" "), 118 ")", sep="") 119 else rn <- rownames(cm) 120 switch(stat.name, 121 t.stat={ 122 prob <- 2 * (1 - pt(abs(ct.diff/vc), df)) 123 }, 124 X2.stat={ 125 prob <- 1 - pchisq((ct.diff/vc)^2, df=1) 126 }) 127 128 if (stat.name=="X2.stat") 129 { 130 retval <- cbind(hyp=beta0, est=ct, stderr=vc, 131 "X^2 value"=(ct.diff/vc)^2, 132 df=df, prob=1 - pchisq((ct.diff/vc)^2, df=1)) 133 dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error", 134 "X^2 value", "DF", "Pr(>|X^2|)")) 135 } 136 else if (stat.name=="t.stat") 137 { 138 retval <- cbind(hyp=beta0, est=ct, stderr=vc, "t value"=ct.diff/vc, 139 df=df, prob=2 * (1 - pt(abs(ct.diff/vc), df))) 140 dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error", 141 "t value", "DF", "Pr(>|t|)")) 142 } 143 144 if (!is.null(conf.int)) 145 { 146 if (conf.int <=0 || conf.int >=1) 147 stop("conf.int should be between 0 and 1. Usual values are 0.95, 0.90") 148 alpha <- 1 - conf.int 149 switch(stat.name, 150 t.stat={ 151 quant <- qt(1 - alpha/2, df) 152 }, 153 X2.stat={ 154 quant <- qt(1 - alpha/2, 100) 155 }) 156 nm <- c(colnames(retval), "Lower.CI", "Upper.CI") 157 retval <- cbind(retval, lower=ct.diff - vc * quant, upper=ct.diff + 158 vc * quant) 159 colnames(retval) <- nm 160 } 161 rownames(retval) <- make.unique(rownames(retval)) 162 retval <- as.data.frame(retval) 163 if(!show.beta0) retval$beta0 <- NULL 164 165 class(retval) <- c("estimable", class(retval)) 166 167 return(retval) 168 } 169} 170 171.wald <- function (obj, cm, 172 beta0=rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm)))) 173{ 174 if (!is.matrix(cm) && !is.data.frame(cm)) 175 cm <- matrix(cm, nrow=1) 176 df <- nrow(cm) 177 if ("geese" %in% class(obj)) 178 { 179 cf <- obj$beta 180 vcv <- obj$vbeta 181 } 182 else if ("gee" %in% class(obj)) 183 { 184 cf <- summary(obj)$coef 185 vcv <- obj$robust.variance 186 } 187 else if ("lm" %in% class(obj)) 188 { 189 cf <- summary.lm(obj)$coefficients[, 1] 190 if ("glm" %in% class(obj)) 191 vcv <- summary(obj)$cov.scaled 192 else 193 vcv <- summary.lm(obj)$cov.unscaled * summary.lm(obj)$sigma^2 194 } 195 else if ("lme" %in% class(obj)) 196 { 197 s.o <- summary(obj) 198 cf <- s.o$tTable[,1] 199 se <- s.o$tTable[, 2] 200 rho <- s.o$cor 201 vcv <- rho * outer(se, se) 202 } 203 else 204 stop("obj must be of class 'lm', 'glm', 'aov', 'gee', 'geese', or 'lme'.") 205 u <- (cm %*% cf)-beta0 206 vcv.u <- cm %*% vcv %*% t(cm) 207 W <- t(u) %*% solve(vcv.u) %*% u 208 prob <- 1 - pchisq(W, df=df) 209 retval <- as.data.frame(cbind(W, df, prob)) 210 names(retval) <- c("X2.stat", "DF", "Pr(>|X^2|)") 211 print(as.data.frame(retval)) 212} 213 214## estimable.mer <- function (obj, cm, beta0, conf.int=NULL, show.beta0, 215## sim.mer=TRUE, n.sim=1000, ...) 216## { 217## if (is.matrix(cm) || is.data.frame(cm)) 218## { 219## cm <- t(apply(cm, 1, .to.est, obj=obj)) 220## } 221## else if(is.list(cm)) 222## { 223## cm <- matrix(.to.est(obj, cm), nrow=1) 224## } 225## else if(is.vector(cm)) 226## { 227## cm <- matrix(.to.est(obj, cm), nrow=1) 228## } 229## else 230## { 231## stop("'cm' argument must be of type vector, list, or matrix.") 232## } 233 234## if(missing(show.beta0)) 235## { 236## if(!missing(beta0)) 237## show.beta0=TRUE 238## else 239## show.beta0=FALSE 240## } 241 242 243## if (missing(beta0)) 244## { 245## beta0 = rep(0, ifelse(is.null(nrow(cm)), 1, nrow(cm))) 246 247## } 248 249## if ("mer" %in% class(obj)) { 250## if(sim.mer) 251## return(est.mer(obj=obj, cm=cm, beta0=beta0, conf.int=conf.int, 252## show.beta0=show.beta0, n.sim=n.sim)) 253 254## stat.name <- "mer" 255## cf <- as.matrix(fixef(obj)) 256## vcv <- as.matrix(vcov(obj)) 257## df <- NA 258## } 259## else { 260## stop("obj is not of class mer") 261## } 262 263## if (is.null(rownames(cm))) 264## rn <- paste("(", apply(cm, 1, paste, collapse=" "), 265## ")", sep="") 266## else rn <- rownames(cm) 267 268## ct <- cm %*% cf[, 1] 269## ct.diff <- cm %*% cf[, 1] - beta0 270## vc <- sqrt(diag(cm %*% vcv %*% t(cm))) 271 272## retval <- cbind(hyp=beta0, est=ct, stderr=vc, "t value"=ct.diff/vc) 273## dimnames(retval) <- list(rn, c("beta0", "Estimate", "Std. Error", 274## "t value")) 275 276## rownames(retval) <- make.unique(rownames(retval)) 277## retval <- as.data.frame(retval) 278## if(!show.beta0) retval$beta0 <- NULL 279 280## class(retval) <- c("estimable", class(retval)) 281 282## return(retval) 283 284## } 285