|
| 1 | + |
| 2 | +#' Levene's Test for Homogeneity of Variance |
| 3 | +#' |
| 4 | +#' Computes Levene's test for homogeneity of variance across groups. |
| 5 | +#' |
| 6 | +#' Let \eqn{X_ij} be the jth observation of X for the ith group. |
| 7 | +#' Let \eqn{Z_ij = |X_ij - X_i|}, where \eqn{X_i} is the mean of X in the ith group. |
| 8 | +#' Levene’s test statistic is |
| 9 | +#' \deqn{ W_0 = \frac{ \sum_i n_i (\bar{Z}_i - \bar{Z})^2 / (g - 1) }{ \sum_i |
| 10 | +#' \sum_j (Z_{ij} - \bar{Z}_i)^2 / \sum_i (n_i - 1) } } |
| 11 | +#' where \eqn{n_i} is the number of observations in group i and g is the number of |
| 12 | +#' groups. |
| 13 | + |
| 14 | +#' @aliases LeveneTest LeveneTest.formula LeveneTest.default |
| 15 | + |
| 16 | +#' @param x response variable for the default method, or a \code{lm} or |
| 17 | +#' \code{formula} object. If \code{y} is a linear-model object or a formula, |
| 18 | +#' the variables on the right-hand-side of the model must all be factors and |
| 19 | +#' must be completely crossed. |
| 20 | + |
| 21 | +#' @param g factor defining groups. |
| 22 | + |
| 23 | +#' @param center The name of a function to compute the center of each group; |
| 24 | +#' \code{mean} gives the original Levene's test; the default, \code{median}, |
| 25 | +#' provides a more robust test (Brown-Forsythe-Test). |
| 26 | + |
| 27 | +#' @param formula a formula of the form \code{lhs ~ rhs} where \code{lhs} gives |
| 28 | +#' the data values and \code{rhs} the corresponding groups. |
| 29 | + |
| 30 | +#' @param data an optional matrix or data frame (or similar: see |
| 31 | +#' \code{\link{model.frame}}) containing the variables in the formula |
| 32 | +#' \code{formula}. By default the variables are taken from |
| 33 | +#' \code{environment(formula)}. |
| 34 | + |
| 35 | +#' @param subset an optional vector specifying a subset of observations to be used. |
| 36 | +#' @param na.action a function which indicates what should happen |
| 37 | +#' when the data contain NAs. Defaults to \code{getOption("na.action")}. |
| 38 | + |
| 39 | +#' @param ... arguments to be passed down, e.g., \code{data} for the |
| 40 | +#' \code{formula} and \code{lm} methods; can also be used to pass arguments to |
| 41 | +#' the function given by \code{center} (e.g., \code{center=mean} and |
| 42 | +#' \code{trim=0.1} specify the 10% trimmed mean). |
| 43 | + |
| 44 | +#' @return returns an object meant to be printed showing the results of the |
| 45 | +#' test. |
| 46 | + |
| 47 | +#' @note This function is rewritten using common R standards based on |
| 48 | +#' car::leveneTest() using the same calculation logic. |
| 49 | + |
| 50 | +#' @author andri.signorell \email{andri@signorell.net}; original version |
| 51 | +#' written by John Fox \email{jfox@@mcmaster.ca} based on a generic version |
| 52 | +#' contributed by Derek Ogle\cr adapted from a response posted by Brian Ripley |
| 53 | +#' to the r-help email list. |
| 54 | + |
| 55 | +#' @seealso \code{\link{fligner.test}} for a rank-based (nonparametric) |
| 56 | +#' \eqn{k}-sample test for homogeneity of variances; \code{\link{mood.test}} |
| 57 | +#' for another rank-based two-sample test for a difference in scale parameters; |
| 58 | +#' \code{\link{var.test}} and \code{\link{bartlett.test}} for parametric tests |
| 59 | +#' for the homogeneity in variance. |
| 60 | +#' |
| 61 | +#' \code{\link[coin:ScaleTests]{ansari_test}} in package \pkg{coin} for exact |
| 62 | +#' and approximate \emph{conditional} p-values for the Ansari-Bradley test, as |
| 63 | +#' well as different methods for handling ties. |
| 64 | +#' @references Fox, J. (2008) \emph{Applied Regression Analysis and Generalized |
| 65 | +#' Linear Models}, Second Edition. Sage. |
| 66 | +#' |
| 67 | +#' Fox, J. and Weisberg, S. (2011) \emph{An R Companion to Applied Regression}, |
| 68 | +#' Second Edition, Sage. |
| 69 | + |
| 70 | +#' @keywords htest |
| 71 | + |
| 72 | +#' @examples |
| 73 | +#' |
| 74 | +#' ## example from ansari.test: |
| 75 | +#' ## Hollander & Wolfe (1973, p. 86f): |
| 76 | +#' ## Serum iron determination using Hyland control sera |
| 77 | +#' serum <- ToLong(data.frame( |
| 78 | +#' ramsay=c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99, |
| 79 | +#' 101, 96, 97, 102, 107, 113, 116, 113, 110, 98), |
| 80 | +#' jung.parekh=c(107, 108, 106, 98, 105, 103, 110, 105, 104, |
| 81 | +#' 100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99) |
| 82 | +#' )) |
| 83 | +#' |
| 84 | +#' LeveneTest(x ~ grp, data=serum) |
| 85 | +#' LeveneTest(x ~ grp, data=serum, center=mean) |
| 86 | +#' LeveneTest(x ~ grp, data=serum, center=mean, trim=0.1) |
| 87 | +#' |
| 88 | +#' LeveneTest( c(rnorm(10), rnorm(10, 0, 2)), factor(rep(c("A","B"), each=10)) ) |
| 89 | +#' |
| 90 | +#' LeveneTest(Ozone ~ Month, data = airquality) |
| 91 | +#' |
| 92 | +#' LeveneTest(count ~ spray, data = InsectSprays) |
| 93 | +#' # Compare this to fligner.test() and bartlett.test() |
| 94 | +#' |
| 95 | + |
| 96 | + |
| 97 | + |
| 98 | +#' @export |
| 99 | +LeveneTest<-function (x,...) |
| 100 | + UseMethod("LeveneTest") |
| 101 | + |
| 102 | +#' @export |
| 103 | +#' @method LeveneTest formula |
| 104 | +#' @name LeveneTest |
| 105 | +LeveneTest.formula<-function (formula,data,subset,na.action,...) { |
| 106 | + |
| 107 | +# kruskal.test |
| 108 | +# if (missing(formula) || (length(formula) != 3L)) |
| 109 | +# stop("'formula' missing or incorrect") |
| 110 | +# m <- match.call(expand.dots = FALSE) |
| 111 | +# if (is.matrix(eval(m$data, parent.frame()))) |
| 112 | +# m$data <- as.data.frame(data) |
| 113 | +# m[[1L]] <- quote(stats::model.frame) |
| 114 | +# mf <- eval(m, parent.frame()) |
| 115 | +# if (length(mf) > 2L) |
| 116 | +# stop("'formula' should be of the form response ~ group") |
| 117 | +#- DNAME <- paste(names(mf), collapse = " by ") |
| 118 | +# y <- LeveneTest(x = mf[[1L]], g = mf[[2L]]) |
| 119 | +# y$data.name <- DNAME |
| 120 | +# y |
| 121 | + |
| 122 | + |
| 123 | +if (missing(formula)|| (length(formula)!=3L)|| (length(attr(terms(formula[-2L]), |
| 124 | +"term.labels"))!=1L)) |
| 125 | + stop("'formula' missing or incorrect") |
| 126 | +m<- match.call(expand.dots=FALSE) |
| 127 | +if (is.matrix(eval(m$data, parent.frame()))) |
| 128 | +m$data<- as.data.frame(data) |
| 129 | +m[[1L]]<- quote(stats::model.frame) |
| 130 | +m$...<-NULL |
| 131 | +mf<- eval(m, parent.frame()) |
| 132 | +DNAME<- paste(names(mf),collapse=" by") |
| 133 | + names(mf)<-NULL |
| 134 | +if (length(mf)>2L) |
| 135 | + stop("'formula' should be of the form response ~ group") |
| 136 | +y<- LeveneTest(x=mf[[1L]],g=mf[[2L]],...) |
| 137 | +y$data.name<-DNAME |
| 138 | +y |
| 139 | + |
| 140 | +} |
| 141 | + |
| 142 | +# example |
| 143 | +# d.set <- data.frame( |
| 144 | +# ramsay = c(111, 107, 100, 99, 102, 106, 109, 108, 104, 99, |
| 145 | +# 101, 96, 97, 102, 107, 113, 116, 113, 110, 98) |
| 146 | +# |
| 147 | +# , jung.parekh = c(107, 108, 106, 98, 105, 103, 110, 105, 104, |
| 148 | +# 100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99) |
| 149 | +# ) |
| 150 | +# |
| 151 | +# |
| 152 | +# ToLong(d.set) |
| 153 | +# |
| 154 | +# LeveneTest(x~grp, ToLong(d.set)) |
| 155 | +# |
| 156 | + |
| 157 | +# LeveneTest.default(x ~ grp, ToLong(d.set)) |
| 158 | +# |
| 159 | +# with(ToLong(d.set), LeveneTest.default(x, grp)) |
| 160 | +# |
| 161 | +# DescTools::LeveneTest(x ~ grp, ToLong(d.set)) |
| 162 | +# |
| 163 | +# debug(LeveneTest.default) |
| 164 | +# |
| 165 | +# var.test(x ~ grp, ToLong(d.set)) |
| 166 | + |
| 167 | +#' @export |
| 168 | +#' @method LeveneTest default |
| 169 | +#' @name LeveneTest |
| 170 | +LeveneTest.default<-function (x,g,center=median,...) { |
| 171 | + |
| 172 | +if (is.list(x)) { |
| 173 | +if (length(x)<2L) |
| 174 | + stop("'x' must be a list with at least 2 elements") |
| 175 | +if (!missing(g)) |
| 176 | + warning("'x' is a list, so ignoring argument 'g'") |
| 177 | +DNAME<- deparse1(substitute(x)) |
| 178 | +x<- lapply(x,function(u)u<-u[complete.cases(u)]) |
| 179 | +if (!all(sapply(x,is.numeric))) |
| 180 | + warning("some elements of 'x' are not numeric and will be coerced to numeric") |
| 181 | +k<- length(x) |
| 182 | +l<- lengths(x) |
| 183 | +if (any(l==0L)) |
| 184 | + stop("all groups must contain data") |
| 185 | +g<-factor(rep.int(seq_len(k),l)) |
| 186 | +x<- unlist(x) |
| 187 | + } |
| 188 | +else { |
| 189 | +if (length(x)!= length(g)) |
| 190 | + stop("'x' and 'g' must have the same length") |
| 191 | +DNAME<- paste(deparse1(substitute(x)),"and", deparse1(substitute(g))) |
| 192 | +OK<- complete.cases(x,g) |
| 193 | +x<-x[OK] |
| 194 | +g<-g[OK] |
| 195 | +g<-factor(g) |
| 196 | +k<- nlevels(g) |
| 197 | +if (k<2L) |
| 198 | + stop("all observations are in the same group") |
| 199 | + } |
| 200 | +n<- length(x) |
| 201 | +if (n<2L) |
| 202 | + stop("not enough observations") |
| 203 | + |
| 204 | +meds<- tapply(x[OK],g[OK],center,...) |
| 205 | +resp<- abs(x-meds[g]) |
| 206 | +ANOVA_TAB<- anova(lm(resp~g)) |
| 207 | + |
| 208 | + rownames(ANOVA_TAB)[2]<-"" |
| 209 | +dots<- deparse(substitute(...)) |
| 210 | + |
| 211 | +dots<- unlist(match.call(expand.dots=FALSE)$...) |
| 212 | +center_x<- deparse(substitute(center)) |
| 213 | +if(!is.null(dots)) |
| 214 | +center_x<- paste0(center_x, |
| 215 | + gettextf("(%s)", |
| 216 | + paste(gettextf("%s=%s", |
| 217 | + names(dots),dots), |
| 218 | +collapse=","))) |
| 219 | + |
| 220 | +STATISTIC<-ANOVA_TAB$`F value`[1] |
| 221 | +PARAMETER<-ANOVA_TAB$Df |
| 222 | +PVAL<-ANOVA_TAB$`Pr(>F)`[1] |
| 223 | + |
| 224 | + names(STATISTIC)<-"F" |
| 225 | + |
| 226 | + names(PARAMETER)<- c("num df","denom df") |
| 227 | + |
| 228 | +RVAL<-list(statistic=STATISTIC,parameter=PARAMETER, |
| 229 | +p.value=PVAL, |
| 230 | +method= gettextf( |
| 231 | +"Levene's Test for Homogeneity of Variance (center = %s)", |
| 232 | +center_x), |
| 233 | +data.name=DNAME,anova_tab=ANOVA_TAB) |
| 234 | + |
| 235 | + class(RVAL)<-"htest" |
| 236 | +return(RVAL) |
| 237 | + |
| 238 | +} |
| 239 | + |
| 240 | + |
| 241 | + |
| 242 | +# |
| 243 | +# with(carData::Moore, DescTools::LeveneTest(conformity, fcategory)) |
| 244 | +# |
| 245 | +# with(carData::Moore, |
| 246 | +# LeveneTest(conformity, interaction(fcategory, partner.status))) |
| 247 | +# |
| 248 | +# LeveneTest(lm(conformity ~ fcategory*partner.status, data = carData::Moore)) |
| 249 | +# |
| 250 | +# LeveneTest(conformity ~ fcategory * partner.status, data = carData::Moore) |
| 251 | +# LeveneTest(conformity ~ fcategory * partner.status, data = Moore, center = mean) |
| 252 | +# LeveneTest(conformity ~ fcategory * partner.status, data = Moore, center = mean, trim = 0.1) |
| 253 | +# |
| 254 | +# |
| 255 | +# |
| 256 | +# |
| 257 | +# |
| 258 | +# LeveneTest(lm(conformity ~ fcategory*partner.status, data = Moore)) |
| 259 | +# |
| 260 | + |
| 261 | + |
| 262 | +# Original |
| 263 | + |
| 264 | + |
| 265 | +# # moved from Rcmdr 13 July 2004 |
| 266 | +# |
| 267 | +# # levene.test.default function slightly modified and generalized from Brian Ripley via R-help |
| 268 | +# # the original generic version was contributed by Derek Ogle |
| 269 | +# # last modified 28 January 2010 by J. Fox |
| 270 | +# |
| 271 | +# LeveneTest <- function (y, ...) { |
| 272 | +# UseMethod("LeveneTest") |
| 273 | +# } |
| 274 | +# |
| 275 | +# LeveneTest.default <- function (y, group, center=median, ...) { # original levene.test |
| 276 | +# |
| 277 | +# if (!is.numeric(y)) |
| 278 | +# stop(deparse(substitute(y)), " is not a numeric variable") |
| 279 | +# |
| 280 | +# if (!is.factor(group)) { |
| 281 | +# warning(deparse(substitute(group)), " coerced to factor.") |
| 282 | +# group <- as.factor(group) |
| 283 | +# } |
| 284 | +# |
| 285 | +# valid <- complete.cases(y, group) |
| 286 | +# meds <- tapply(y[valid], group[valid], center, ...) |
| 287 | +# resp <- abs(y - meds[group]) |
| 288 | +# table <- anova(lm(resp ~ group))[, c(1, 4, 5)] |
| 289 | +# rownames(table)[2] <- " " |
| 290 | +# dots <- deparse(substitute(...)) |
| 291 | +# |
| 292 | +# attr(table, "heading") <- paste("Levene's Test for Homogeneity of Variance (center = ", |
| 293 | +# deparse(substitute(center)), if(!(dots == "NULL")) paste(":", dots), ")", sep="") |
| 294 | +# table |
| 295 | +# } |
| 296 | +# |
| 297 | +# |
| 298 | +# LeveneTest.formula <- function(formula, data, ...) { |
| 299 | +# form <- formula |
| 300 | +# mf <- if (missing(data)) model.frame(form) else model.frame(form, data) |
| 301 | +# if (any(sapply(2:dim(mf)[2], function(j) is.numeric(mf[[j]])))) |
| 302 | +# stop("Levene's test is not appropriate with quantitative explanatory variables.") |
| 303 | +# y <- mf[,1] |
| 304 | +# if(dim(mf)[2]==2) group <- mf[,2] |
| 305 | +# else { |
| 306 | +# if (length(grep("\\+ | \\| | \\^ | \\:",form))>0) stop("Model must be completely crossed formula only.") |
| 307 | +# group <- interaction(mf[,2:dim(mf)[2]]) |
| 308 | +# } |
| 309 | +# LeveneTest.default(y = y, group=group, ...) |
| 310 | +# } |
| 311 | +# |
| 312 | +# |
| 313 | +# |
| 314 | +# LeveneTest.lm <- function(y, ...) { |
| 315 | +# LeveneTest.formula(formula(y), data=model.frame(y), ...) |
| 316 | +# } |
| 317 | +# |
| 318 | +# |
| 319 | +# |
| 320 | + |
| 321 | + |
| 322 | + |