diff --git a/DESCRIPTION b/DESCRIPTION index 177f51b92..ba9dc2e9e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DescTools Type: Package Title: Tools for Descriptive Statistics -Version: 0.99.49.1 -Date: 2023-04-30 +Version: 0.99.50 +Date: 2023-09-06 Authors@R: c( person("Andri", "Signorell", email = "andri@signorell.net", role = c("aut", "cre")), person("Ken" , "Aho", role = c("ctb")), diff --git a/NAMESPACE b/NAMESPACE index 6fed2b0c5..a8c4584d0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,7 +7,7 @@ export( "Clockwise","Closest","Coalesce","CourseData", "CochranArmitageTest", "CochranQTest", -"CoefVar","CohenD","CohenKappa","CollapseTable","ColorLegend", +"CoefVar","CoefVarCI","CohenD","CohenKappa","CollapseTable","ColorLegend", "ColToGray","ColToGrey","ColToHex","ColToHsv","ColToRgb","ColToOpaque","RgbToHex", "CombN","CombPairs","CombSet","ConDisPairs","Conf","ConnLines","ContCoef","Contrasts","CorCI", "CorPolychor","CramerV","CramerVonMisesTest","CronbachAlpha","CutQ","Day","Day<-","dBenf", @@ -89,7 +89,7 @@ importFrom("stats", ".nknots.smspl", "IQR", "acf", "addmargins", "anova", "aov", "uniroot", "var", "weighted.mean", "weights", "wilcox.test", "xtabs", "as.dendrogram", "dist", "hclust", "order.dendrogram", "smooth.spline","var.test", "median", "reorder","AIC", "BIC", "confint", "deviance", "df", "logLik","lm.wfit","profile", "coef", "coefficients", "cov2cor", "df.residual","nobs", "vcov","update","relevel","rbeta", "rexp", "qexp","reshape", - "confint.default", "ecdf", "glm", "reformulate", "aggregate") + "confint.default", "ecdf", "glm", "reformulate", "aggregate", "dchisq") importFrom("data.table", "frankv", "frank") diff --git a/NEWS b/NEWS index 2aab23d36..1dc10f0ae 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,24 @@ -DescTools 0.99.49 (2023-04-28) + +DescTools 0.99.50 (2023-09-02) +------------------------------ + +NEW FUNCTIONS ADDED: + * The confidence intervals for the coefficient of variation are calculated + by the new function CoefVar(). Additional methods have been implemented. + So far these are: "vangel","mckay","nct","naive". + +UPDATED FUNCTIONS: + * Mode() returns NA for a vector of length = 1. + * CoefVar() confidence intervals have been organised in a new function. + * DFLOAT replaced by DBLE everywhere in Fortran code. + + + + + +DescTools 0.99.49 (2023-05-14) ------------------------------ UPDATED FUNCTIONS: @@ -8,7 +26,7 @@ UPDATED FUNCTIONS: (Credits to Michael Chirico in https://github.com/AndriSignorell/DescTools/pull/102) - * BoxedText() uses xy.coords now. + * BoxedText() supports xy.coords now. * Two new methods for MultinomCI() have been implemented (Fitzpatrick-Scott and Quesensberry-Hurst). * BinomCI() no longer reports p_tilde but the usual x/n as @@ -17,6 +35,8 @@ UPDATED FUNCTIONS: which the method specific non-standard point estimator can be returned. (Thanks to julienbio99, MrJerryTAO and Vilmantas for the research and discussion) + * VarTest() gets a more exact calculation for twosided p values + in the case of the one sample test. BUGFIXES: * MultinomCI method goodman uses the correct chisquare quantile now. diff --git a/R/StatsAndCIs.r b/R/StatsAndCIs.r index 661f9feb4..cd0c2daf4 100644 --- a/R/StatsAndCIs.r +++ b/R/StatsAndCIs.r @@ -1186,7 +1186,9 @@ Mode <- function(x, na.rm=FALSE) { if(length(x) == 1L) # only one value in x, x is the mode - return(structure(x, freq = 1L)) + # return(structure(x, freq = 1L)) + # changed to: only one value in x, no mode defined + return(structure(NA_real_, freq = NA_integer_)) # we don't have NAs so far, either there were then we've already stopped # or they've been stripped above @@ -4476,6 +4478,24 @@ CohenD <- function(x, y=NULL, pooled = TRUE, correct = FALSE, conf.level = NA, n } + + +# find non-centrality parameter for the F-distribution +ncparamF <- function(type1, type2, nu1, nu2) { + + # author Ali Baharev + + # Returns the noncentrality parameter of the noncentral F distribution + # if probability of Type I and Type II error, degrees of freedom of the + # numerator and the denominator in the F test statistics are given. + + + .C("fpow", PACKAGE = "DescTools", as.double(type1), as.double(type2), + as.double(nu1), as.double(nu2), lambda=double(1))$lambda +} + + + .nctCI <- function(t, df, conf) { alpha <- 1 - conf @@ -4574,7 +4594,7 @@ CoefVar <- function (x, ...) { -CoefVar.lm <- function (x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) { +CoefVar.lm <- function (x, unbiased = FALSE, na.rm = FALSE, ...) { # source: http://www.ats.ucla.edu/stat/mult_pkg/faq/general/coefficient_of_variation.htm @@ -4592,11 +4612,12 @@ CoefVar.lm <- function (x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ... res <- res * ((1 - (1/(4 * (n - 1))) + (1/n) * res^2) + (1/(2 * (n - 1)^2))) } - if (!is.na(conf.level)) { - ci <- .nctCI(sqrt(n)/res, df = n - 1, conf = conf.level) - res <- c(est = res, low.ci = unname(sqrt(n)/ci["upr.ci"]), - upr.ci = unname(sqrt(n)/ci["lwr.ci"])) - } + + # if (!is.na(conf.level)) { + # ci <- .nctCI(sqrt(n)/res, df = n - 1, conf = conf.level) + # res <- c(est = res, low.ci = unname(sqrt(n)/ci["upr.ci"]), + # upr.ci = unname(sqrt(n)/ci["lwr.ci"])) + # } return(res) @@ -4607,33 +4628,133 @@ CoefVar.lm <- function (x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ... # dv <- unname(nlme::getResponse(x)) -CoefVar.default <- function (x, weights=NULL, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) { +# CoefVar.default <- function (x, weights=NULL, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) { +# +# if(is.null(weights)){ +# if(na.rm) x <- na.omit(x) +# res <- SD(x) / Mean(x) +# n <- length(x) +# +# } +# else { +# res <- SD(x, weights = weights) / Mean(x, weights = weights) +# n <- sum(weights) +# +# } +# +# if(unbiased) { +# res <- res * ((1 - (1/(4*(n-1))) + (1/n) * res^2)+(1/(2*(n-1)^2))) +# } +# +# if(!is.na(conf.level)){ +# ci <- .nctCI(sqrt(n)/res, df = n-1, conf = conf.level) +# res <- c(est=res, low.ci= unname(sqrt(n)/ci["upr.ci"]), upr.ci= unname(sqrt(n)/ci["lwr.ci"])) +# } +# +# return(res) +# +# } - if(is.null(weights)){ - if(na.rm) x <- na.omit(x) - res <- SD(x) / Mean(x) - n <- length(x) + + +CoefVarCI <- function (K, n, conf.level = 0.95, + sides = c("two.sided", "left", "right"), + method = c("nct","vangel","mckay","verrill","naive")) { + + # Description of confidence intervals + # https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/coefvacl.htm + + + .iCoefVarCI <- Vectorize(function(K, n, conf.level=0.95, + sides = c("two.sided", "left", "right"), + method = c("vangel","mckay","verrill","nct","naive")) { - } - else { - res <- SD(x, weights = weights) / Mean(x, weights = weights) - n <- sum(weights) + method <- match.arg(method) + sides <- match.arg(sides, choices = c("two.sided", "left", "right"), + several.ok = FALSE) - } + # double alpha in case of one-sided intervals in order to be able + # to generally calculate twosided intervals and select afterwards.. + if (sides != "two.sided") + conf.level <- 1 - 2 * (1 - conf.level) + + alpha <- 1 - conf.level + + df <- n - 1 + u1 <- qchisq(1-alpha/2, df) + u2 <- qchisq(alpha/2, df) + + switch(method, verrill = { + CI.lower <- 0 + CI.upper <- 1 + + }, vangel = { + CI.lower <- K / sqrt(((u1+2)/n - 1) * K^2 + u1/df) + CI.upper <- K / sqrt(((u2+2)/n - 1) * K^2 + u2/df) + + }, mckay = { + CI.lower <- K / sqrt((u1/n - 1) * K^2 + u1/df) + CI.upper <- K / sqrt((u2/n - 1) * K^2 + u2/df) + + }, nct = { + ci <- .nctCI(sqrt(n)/K, df = df, conf = conf.level) + CI.lower <- unname(sqrt(n)/ci[2]) + CI.upper <- unname(sqrt(n)/ci[1]) + + }, naive = { + CI.lower <- K * sqrt(df / u1) + CI.upper <- K * sqrt(df / u2) + } + ) + + ci <- c(est = K, + lwr.ci = CI.lower, # max(0, CI.lower), + upr.ci = CI.upper) # min(1, CI.upper)) + + if (sides == "left") + ci[3] <- Inf + + else if (sides == "right") + ci[2] <- -Inf + + return(ci) + + }) + - if(unbiased) { - res <- res * ((1 - (1/(4*(n-1))) + (1/n) * res^2)+(1/(2*(n-1)^2))) - } + sides <- match.arg(sides) + method <- match.arg(method) + + res <- t(.iCoefVarCI(K=K, n=n, method=method, sides=sides, conf.level = conf.level)) - if(!is.na(conf.level)){ - ci <- .nctCI(sqrt(n)/res, df = n-1, conf = conf.level) - res <- c(est=res, low.ci= unname(sqrt(n)/ci["upr.ci"]), upr.ci= unname(sqrt(n)/ci["lwr.ci"])) + return(res) + +} + + +CoefVar.default <- function (x, weights = NULL, unbiased = FALSE, + na.rm = FALSE, ...) { + + + if (is.null(weights)) { + if (na.rm) + x <- na.omit(x) + res <- SD(x)/Mean(x) + n <- length(x) + + } else { + res <- SD(x, weights = weights)/Mean(x, weights = weights) + n <- sum(weights) + } + + if (unbiased) { + res <- res * ((1 - (1/(4 * (n - 1))) + (1/n) * res^2) + (1/(2 * (n - 1)^2))) } return(res) - } + # aus agricolae: Variations Koeffizient aus aov objekt # # CoefVar.aov <- function(x){ @@ -4641,7 +4762,7 @@ CoefVar.default <- function (x, weights=NULL, unbiased = FALSE, conf.level = NA, # } -CoefVar.aov <- function (x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) { +CoefVar.aov <- function (x, unbiased = FALSE, na.rm = FALSE, ...) { # source: http://www.ats.ucla.edu/stat/mult_pkg/faq/general/coefficient_of_variation.htm @@ -4659,11 +4780,12 @@ CoefVar.aov <- function (x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, .. res <- res * ((1 - (1/(4 * (n - 1))) + (1/n) * res^2) + (1/(2 * (n - 1)^2))) } - if (!is.na(conf.level)) { - ci <- .nctCI(sqrt(n)/res, df = n - 1, conf = conf.level) - res <- c(est = res, low.ci = unname(sqrt(n)/ci["upr.ci"]), - upr.ci = unname(sqrt(n)/ci["lwr.ci"])) - } + + # if (!is.na(conf.level)) { + # ci <- .nctCI(sqrt(n)/res, df = n - 1, conf = conf.level) + # res <- c(est = res, low.ci = unname(sqrt(n)/ci["upr.ci"]), + # upr.ci = unname(sqrt(n)/ci["lwr.ci"])) + # } return(res) @@ -5236,6 +5358,24 @@ ContCoef <- function(x, y = NULL, correct = FALSE, ...) { } + +.ncchisqCI <- function(chisq, df, conf) { + + alpha <- 1 - conf + probs <- c(alpha/2, 1 - alpha/2) + + ncp <- suppressWarnings(optim(par = 1.1 * rep(chisq, 2), fn = function(x) { + p <- pchisq(q = chisq, df = df, ncp = x) + abs(max(p) - probs[2]) + abs(min(p) - probs[1]) + }, control = list(abstol = 0.000000001))) + + chisq_ncp <- unname(sort(ncp$par)) + + return(chisq_ncp) + +} + + CramerV <- function(x, y = NULL, conf.level = NA, method = c("ncchisq", "ncchisqadj", "fisher", "fisheradj"), correct=FALSE, ...){ @@ -5376,19 +5516,6 @@ CramerV <- function(x, y = NULL, conf.level = NA, -ncparamF <- function(type1, type2, nu1, nu2) { - - # author Ali Baharev - - # Returns the noncentrality parameter of the noncentral F distribution - # if probability of Type I and Type II error, degrees of freedom of the - # numerator and the denominator in the F test statistics are given. - - - .C("fpow", PACKAGE = "DescTools", as.double(type1), as.double(type2), - as.double(nu1), as.double(nu2), lambda=double(1))$lambda -} - YuleQ <- function(x, y = NULL, ...){ @@ -8180,8 +8307,10 @@ print.HoeffD <- function(x, ...) } -# find non-centrality parameter for the F-distribution -ncparamF <- function(type1, type2, nu1, nu2){ - .C("fpow", PACKAGE = "DescTools", as.double(type1), as.double(type2), as.double(nu1), as.double(nu2), lambda=double(1))$lambda -} + + + + + + diff --git a/R/Tests.r b/R/Tests.r index 26e8cfd41..46267171e 100644 --- a/R/Tests.r +++ b/R/Tests.r @@ -723,6 +723,59 @@ VarTest <- function(x, ...) UseMethod("VarTest") VarTest.default <- function (x, y = NULL, alternative = c("two.sided", "less", "greater"), ratio = 1, sigma.squared = 1, conf.level = 0.95, ...) { + + TwoSided_pval <- function(x, DF){ + + # https://stats.stackexchange.com/questions/195469/calculating-p-values-for-two-tail-test-for-population-variance + + # What you are dealing with in this question is a two-sided + # variance test, which is a specific case of a two-sided test + # with an asymmetric null distribution. The p-value is the total + # area under the null density for all values in the lower and + # upper tails of that density that are at least as "extreme" + # (i.e., at least as conducive to the alternative hypothesis) + # as the observed test statistic. Because this test has an + # asymmetric null distribution, we need to specify exactly + # what we mean by "extreme". + # + # Lowest-density p-value calculation: The most sensible thing + # method of two-sided hypothesis testing is to interpret + # "more extreme" as meaning a lower value of the null density. + # This is the interpretation used in a standard likelihood-ratio + # (LR) test. Under this method , the p-value is the probability of + # falling in the "lowest density region", where the density + # cut-off is the density at the observed test statistic. With + # an asymmetric null distribution, this leads you to a p-value + # calculated with unequal tails. + + # example: + # TwoSided_pval(15.35667, DF=17) + # TwoSided_pval(14.6489, DF=17) + + + InvDChisq <- function(x2, DF){ + + fun <- function(x) dchisq(x, df=DF) - dchisq(x2, df=DF) + # the mode of chisq distribution + mod_x2 <- DF-2 + + if(x2 < mod_x2){ + # don't know how to set the right boundary in a sensible way + # the treshold 1e12 is selected randomly + # benchmark show no difference in performance between 1e6 and 1e12 + UnirootAll(fun, interval = c(mod_x2, 1e12)) + } else { + UnirootAll(fun, interval = c(0, mod_x2)) + } + } + + + pchisq(x, df = DF, lower.tail = x < DF-2) + + pchisq(InvDChisq(x, DF), df = DF, lower.tail=!x < DF-2) + + } + + if(is.null(y)){ # perform a one sample variance test @@ -755,11 +808,14 @@ VarTest.default <- function (x, y = NULL, alternative = c("two.sided", "less", " } else if (alternative == "greater") { pval <- pchisq(xstat, df, lower.tail = FALSE) cint <- c(df * vx/qchisq((1 - conf.level), df, lower.tail = FALSE), Inf) + } else { - pval <- 2 * min(pchisq(xstat, df), pchisq(xstat, df, lower.tail = FALSE)) + # this is a "quick-and-nasty" approximation, let's use a better one.. + # pval <- 2 * min(pchisq(xstat, df), + # pchisq(xstat, df, lower.tail = FALSE)) + pval <- TwoSided_pval(xstat, df) + alpha <- 1 - conf.level - cint <- qt(1 - alpha/2, df) - cint <- xstat + c(-cint, cint) cint <- df * vx / c(qchisq((1 - conf.level)/2, df, lower.tail = FALSE), qchisq((1 - conf.level)/2, df)) } diff --git a/man/CoefVar.Rd b/man/CoefVar.Rd index cfb859b61..685594219 100644 --- a/man/CoefVar.Rd +++ b/man/CoefVar.Rd @@ -3,21 +3,27 @@ \alias{CoefVar.lm} \alias{CoefVar.aov} \alias{CoefVar.default} +\alias{CoefVarCI} %- Also NEED an '\alias' for EACH other topic documented here. \title{Coefficient of Variation %% ~~function to do ... ~~ } -\description{Calculates the coefficient of variation and its confidence limits using the noncentral \code{t}-distribution.. +\description{Calculates the coefficient of variation and its confidence limits using various methods. %% ~~ A concise (1-5 lines) description of what the function does. ~~ } \usage{ CoefVar(x, ...) -\method{CoefVar}{lm}(x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) +\method{CoefVar}{lm}(x, unbiased = FALSE, na.rm = FALSE, ...) -\method{CoefVar}{aov}(x, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) +\method{CoefVar}{aov}(x, unbiased = FALSE, na.rm = FALSE, ...) -\method{CoefVar}{default}(x, weights = NULL, unbiased = FALSE, conf.level = NA, na.rm = FALSE, ...) +\method{CoefVar}{default}(x, weights = NULL, unbiased = FALSE, + na.rm = FALSE, ...) + +CoefVarCI(K, n, conf.level = 0.95, + sides = c("two.sided", "left", "right"), + method = c("nct","vangel","mckay","verrill","naive")) } %- maybe also 'usage' for other objects documented here. @@ -31,9 +37,19 @@ CoefVar(x, ...) \item{unbiased}{logical value determining, if a bias correction should be used (see. details). Default is FALSE. %% Reference??? %% ~~Describe \code{unbiased} here~~ } - \item{conf.level}{confidence level of the interval. + \item{K}{the coefficient of variation as calculated by \code{CoefVar()}. +%% ~~Describe \code{x} here~~ +} + \item{n}{the number of observations used for calculating the coefficient of variation.} + + \item{conf.level}{confidence level of the interval. Defaults to 0.95. %% ~~Describe \code{conf.level} here~~ } + \item{sides}{a character string specifying the side of the confidence interval, must be one of \code{"two.sided"} (default), +\code{"left"} or \code{"right"}. You can specify just the initial letter. \code{"left"} would be analogue to a hypothesis of +\code{"greater"} in a \code{t.test}.} +\item{method}{character string specifing the method to use for calculating the confidence intervals, can be one out of: + \code{"nct"} (default), \code{"vangel"}, \code{"mckay"}, \code{"verrill"} (currently not yet implemented) and \code{"naive"}. Abbreviation of method is accepted. See details.} \item{na.rm}{logical. Should missing values be removed? Defaults to FALSE. %% ~~Describe \code{na.rm} here~~ } @@ -42,6 +58,23 @@ CoefVar(x, ...) \details{ In order for the coefficient of variation to be an unbiased estimate of the true population value, the coefficient of variation is corrected as: \deqn{ CV_{korr} = CV \cdot \left( 1 - \frac{1}{4\cdot(n-1)} + \frac{1}{n} \cdot CV^2 + \frac{1}{2 \cdot (n-1)^2} \right) } + +For determining\verb{ }\bold{the confidence intervals}\verb{ } for the coefficient of variation a number of methods have been proposed. \code{CoefVarCI()} currently supports five different methods. +The details for the methods are given in the specific references. + +The \bold{"naive" method} \verb{ } +is based on dividing the standard confidence limit for the standard deviation by the sample mean. + +\bold{McKay's} \verb{ } +approximation is asymptotically exact as n goes to infinity. McKay recommends this approximation only if the coefficient of variation is less than 0.33. Note that if the coefficient of variation is greater than 0.33, either the normality of the data is suspect or the probability of negative values in the data is non-neglible. In this case, McKay's approximation may not be valid. Also, it is generally recommended that the sample size should be at least 10 before using McKay's approximation. + +\bold{Vangel's modified McKay method} \verb{ } +is more accurate than the McKay in most cases, particilarly for small samples.. According to Vangel, the unmodified McKay is only more accurate when both the coefficient of variation and alpha are large. However, if the coefficient of variation is large, then this implies either that the data contains negative values or the data does not follow a normal distribution. In this case, neither the McKay or the modified McKay should be used. +In general, the Vangel's modified McKay method is recommended over the McKay method. It generally provides good approximations as long as the data is approximately normal and the coefficient of variation is less than 0.33. This is the default method. + +See also: https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/coefvacl.htm + +\bold{nct} \verb{ }uses the noncentral t-distribution to calculate the confidence intervals. See Smithson (2003). } \value{if no confidence intervals are requested: the estimate as numeric value (without any name)\cr\cr @@ -52,17 +85,22 @@ else a named numeric vector with 3 elements } \references{ +McKay, A. T. (1932). Distribution of the coefficient of variation and the extended +\emph{t} distribution, \emph{Journal of the Royal Statistical Society}, \emph{95}, 695--698. + Johnson, B. L., Welch, B. L. (1940). Applications of the non-central t-distribution. \emph{Biometrika}, 31, 362--389. +Mark Vangel (1996) Confidence Intervals for a Normal Coefficient of Variation, \emph{American Statistician}, Vol. 15, No. 1, pp. 21-26. + Kelley, K. (2007). Sample size planning for the coefcient of variation from the accuracy in parameter estimation approach. \emph{Behavior Research Methods, 39} (4), 755-766 Kelley, K. (2007). Constructing confidence intervals for standardized effect sizes: Theory, application, and implementation. \emph{Journal of Statistical Software, 20} (8), 1-24 -McKay, A. T. (1932). Distribution of the coefficient of variation and the extended -\emph{t} distribution, \emph{Journal of the Royal Statistical Society}, \emph{95}, 695--698. - Smithson, M.J. (2003) \emph{Confidence Intervals, Quantitative Applications in the Social Sciences Series}, No. 140. Thousand Oaks, CA: Sage. pp. 39-41 +Steve Verrill (2003) Confidence Bounds for Normal and Lognormal Distribution Coefficients of Variation, \emph{Research Paper 609}, USDA Forest Products Laboratory, Madison, Wisconsin. + +Verrill, S. and Johnson, R.A. (2007) Confidence Bounds and Hypothesis Tests for Normal Distribution Coefficients of Variation, \emph{Communications in Statistics Theory and Methods}, Volume 36, No. 12, pp 2187-2206. } \author{Andri Signorell , \cr Michael Smithson (noncentral-t) @@ -85,6 +123,11 @@ CoefVar(x, conf.level=0.95) # Coefficient of variation for a linear model r.lm <- lm(Fertility ~ ., swiss) CoefVar(r.lm) + +# the function is vectorized, so arguments are recyled... +# https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/coefvacl.htm +CoefVarCI(K = 0.00246, n = 195, method="vangel", + sides="two.sided", conf.level = c(.5,.8,.9,.95,.99,.999)) } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. diff --git a/man/DurbinWatsonTest.Rd b/man/DurbinWatsonTest.Rd index 762d37bc0..b09147e30 100644 --- a/man/DurbinWatsonTest.Rd +++ b/man/DurbinWatsonTest.Rd @@ -1,6 +1,5 @@ \name{DurbinWatsonTest} \alias{DurbinWatsonTest} -\encoding{latin1} \title{Durbin-Watson Test} \description{ Performs the Durbin-Watson test for autocorrelation of disturbances. diff --git a/man/GoodmanKruskalTau.Rd b/man/GoodmanKruskalTau.Rd index aaa94030d..09f315d28 100644 --- a/man/GoodmanKruskalTau.Rd +++ b/man/GoodmanKruskalTau.Rd @@ -92,6 +92,28 @@ GoodmanKruskalTau(tab, direction="row", conf.level=0.95) GoodmanKruskalTau(tab, direction="column", conf.level=0.95) # reduce both to: Phi(tab)^2 + + +# example 1 in Liebetrau (1983) + +tt <- matrix(c(549,93,233,119,225,455,402, + 212,124,78,42,41,12,132, + 54,54,33,13,46,7,153), ncol=3, + dimnames=list(rownames=c("Gov", "Mil", "Edu", "Eco", "Intel", "Rel", "For"), + colnames=c("One", "Two", "Multi"))) + +GoodmanKruskalTau(tt, direction = "row", conf.level = 0.95) +GoodmanKruskalTau(tt, direction = "column", conf.level = 0.95) + + +# SPSS +ttt <- matrix(c(225,53,206,3,1,12), nrow=3, + dimnames=list(rownames=c("right","center", "left"), + colnames=c("us","ussr"))) + +round(GoodmanKruskalTau(ttt, direction = "r", con=0.95), d=3) +round(GoodmanKruskalTau(ttt, direction = "c"), d=3) + } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. diff --git a/src/ADinf.o b/src/ADinf.o index 5d93fbc85..40416911a 100644 Binary files a/src/ADinf.o and b/src/ADinf.o differ diff --git a/src/AnDarl.o b/src/AnDarl.o index 3c76d7944..6f11c92f3 100644 Binary files a/src/AnDarl.o and b/src/AnDarl.o differ diff --git a/src/DescTools.dll b/src/DescTools.dll index 6abbf362e..4538f5a53 100644 Binary files a/src/DescTools.dll and b/src/DescTools.dll differ diff --git a/src/RcppExports.o b/src/RcppExports.o index cb659ca30..fbe3c8d20 100644 Binary files a/src/RcppExports.o and b/src/RcppExports.o differ diff --git a/src/aux_fct.o b/src/aux_fct.o index 59df5643b..6e5eb5d8a 100644 Binary files a/src/aux_fct.o and b/src/aux_fct.o differ diff --git a/src/between.o b/src/between.o index f9a6c6421..a1cbad304 100644 Binary files a/src/between.o and b/src/between.o differ diff --git a/src/cidxcn.o b/src/cidxcn.o index d8c3a3c33..d516c514a 100644 Binary files a/src/cidxcn.o and b/src/cidxcn.o differ diff --git a/src/conc.o b/src/conc.o index d2139e7e2..6ecd3662c 100644 Binary files a/src/conc.o and b/src/conc.o differ diff --git a/src/exactsum.o b/src/exactsum.o index 5855ef10e..50c4962df 100644 Binary files a/src/exactsum.o and b/src/exactsum.o differ diff --git a/src/extremes.o b/src/extremes.o index 183ca5989..8f7b69924 100644 Binary files a/src/extremes.o and b/src/extremes.o differ diff --git a/src/fdwilcox.o b/src/fdwilcox.o index daaa0232c..7ea5d1b4f 100644 Binary files a/src/fdwilcox.o and b/src/fdwilcox.o differ diff --git a/src/fpow.o b/src/fpow.o index a3b73288c..3d61ad57f 100644 Binary files a/src/fpow.o and b/src/fpow.o differ diff --git a/src/gompertz.o b/src/gompertz.o index 5b33f3afa..520783b87 100644 Binary files a/src/gompertz.o and b/src/gompertz.o differ diff --git a/src/hoeffd.f b/src/hoeffd.f index 7b46b5347..edde8d898 100644 --- a/src/hoeffd.f +++ b/src/hoeffd.f @@ -61,7 +61,7 @@ subroutine hoeff(x, y, n, d, aad, maxad, rx, ry, rj) s = 0d0 aad = 0d0 maxad = 0d0 - z = dfloat(n) + z = dble(n) do23018 i=1,n rxi = rx(i) ryi = ry(i) diff --git a/src/hoeffd.o b/src/hoeffd.o index 3941dd7a2..fb139b0c3 100644 Binary files a/src/hoeffd.o and b/src/hoeffd.o differ diff --git a/src/init.o b/src/init.o index a5184195a..aa3b0b3e7 100644 Binary files a/src/init.o and b/src/init.o differ diff --git a/src/jtpdf.f b/src/jtpdf.f index 644628ed4..cabcd1682 100644 --- a/src/jtpdf.f +++ b/src/jtpdf.f @@ -18,10 +18,10 @@ subroutine jtpdf(mxsum, pdf, ng, cgsize, pdf0, pdf1) m = cgsize(ng-1) - cgsize(ng) n = cgsize(ng) mn1 = m*n - dm = dfloat(m) - dn = dfloat(n) + dm = dble(m) + dn = dble(n) do 10 i = 0, mn1 - di = dfloat(i) + di = dble(i) pdf(i+1) = fdwilcox(di, dm, dn) 10 continue @@ -35,12 +35,12 @@ subroutine jtpdf(mxsum, pdf, ng, cgsize, pdf0, pdf1) m = cgsize(g) - cgsize(g+1) n = cgsize(g+1) mn0 = m*n - dm = dfloat(m) - dn = dfloat(n) + dm = dble(m) + dn = dble(n) c call intpr(" m",2,m,1) c call intpr(" n",2,n,1) do 30 i = 0, mn0 - di = dfloat(i) + di = dble(i) dmw = fdwilcox(di, dm, dn) c call dblepr(" dmw", 4, dmw, 1) pdf0(i+1) = dmw diff --git a/src/jtpdf.o b/src/jtpdf.o index 26a261146..cd86f4f9c 100644 Binary files a/src/jtpdf.o and b/src/jtpdf.o differ diff --git a/src/ks.o b/src/ks.o index 129a36b8d..a2623c100 100644 Binary files a/src/ks.o and b/src/ks.o differ diff --git a/src/moments.o b/src/moments.o index 95bfa38fd..963673f14 100644 Binary files a/src/moments.o and b/src/moments.o differ diff --git a/src/pan.o b/src/pan.o index 43af3af5a..4522ecfc7 100644 Binary files a/src/pan.o and b/src/pan.o differ diff --git a/src/pointinpolygon.o b/src/pointinpolygon.o index c848e7198..ccb5bb6d1 100644 Binary files a/src/pointinpolygon.o and b/src/pointinpolygon.o differ diff --git a/src/roman2int.o b/src/roman2int.o index cb996f3d7..22a298874 100644 Binary files a/src/roman2int.o and b/src/roman2int.o differ diff --git a/src/tbrm.o b/src/tbrm.o index 7a4fd6350..68a02074b 100644 Binary files a/src/tbrm.o and b/src/tbrm.o differ diff --git a/src/wgt_himed.o b/src/wgt_himed.o index 83752cd39..d9a357676 100644 Binary files a/src/wgt_himed.o and b/src/wgt_himed.o differ diff --git a/tests/misc.R b/tests/misc.R index d469a27e1..ae39e906b 100644 --- a/tests/misc.R +++ b/tests/misc.R @@ -17,7 +17,7 @@ set.seed(45) (z <- as.numeric(names(w <- table(x <- sample(-10:20, size=50, r=TRUE))))) stopifnot(all( - identical(Mode(5), structure(5, freq = 1L)) + identical(Mode(5), structure(NA_real_, freq = NA_integer_)) , identical(Mode(NA), structure(NA_real_, freq = NA_integer_)) , identical(Mode(c(NA, NA)), structure(NA_real_, freq = NA_integer_)) , identical(Mode(c(NA, 0:5)), structure(NA_real_, freq = NA_integer_)) @@ -76,26 +76,28 @@ stopifnot(all( meth <- c("wald","waldcc","hal","jp","mee","mn","score","scorecc","ha","ac","blj") +# use all(IsZero(x - y)) to take into account numerical properties of +# certain operating systems (especially PowerPC) stopifnot(all( - identical(unname(round(BinomDiffCI(56, 70, 48, 80, method = meth), 4)[, -1]), - cbind(c(0.0575, 0.0441, 0.0535, 0.0531, 0.0534, - 0.0528, 0.0524, 0.0428, 0.0494, 0.0525, 0.054), - c(0.3425, 0.3559, 0.3351, 0.3355, 0.3377, - 0.3382, 0.3339, 0.3422, 0.3506, 0.3358, 0.34))), - identical(unname(round(BinomDiffCI(9, 10, 3, 10, method = meth), 4)[, -1]), + all(IsZero(unname(round(BinomDiffCI(56, 70, 48, 80, method = meth), 4)[, -1]) - + cbind(c(0.0575, 0.0441, 0.0535, 0.0531, 0.0534, + 0.0528, 0.0524, 0.0428, 0.0494, 0.0525, 0.054), + c(0.3425, 0.3559, 0.3351, 0.3355, 0.3377, + 0.3382, 0.3339, 0.3422, 0.3506, 0.3358, 0.34)))), + all(IsZero(unname(round(BinomDiffCI(9, 10, 3, 10, method = meth), 4)[, -1]) - cbind(c(0.2605, 0.1605, 0.1777, 0.176, 0.1821, 0.17, 0.1705, 0.1013, 0.1922, 0.16, 0.1869), c(0.9395, 1, 0.8289, 0.8306, 0.837, 0.8406, - 0.809, 0.8387, 1, 0.84, 0.904))), - identical(unname(round(BinomDiffCI(10, 10, 0, 20, method = meth), 4)[, -1]), - cbind(c(1, 0.925, 0.7482, 0.7431, 0.7224, 0.7156, - 0.6791, 0.6014, 0.95, 0.6922, 0.7854), - c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))), - identical(unname(round(BinomDiffCI(84, 101, 89, 105, method = meth), 4)[, -1]), - cbind(c(-0.1162, -0.1259, -0.1152, -0.116, -0.1188, - -0.1191, -0.1177, -0.1245, -0.1216, -0.1168, -0.117), - c(0.0843, 0.094, 0.0834, 0.0843, 0.0858, 0.086, 0.0851, - 0.0918, 0.0898, 0.085, 0.0852))) + 0.809, 0.8387, 1, 0.84, 0.904)))), + all(IsZero(unname(round(BinomDiffCI(10, 10, 0, 20, method = meth), 4)[, -1]) - + cbind(c(1, 0.925, 0.7482, 0.7431, 0.7224, 0.7156, + 0.6791, 0.6014, 0.95, 0.6922, 0.7854), + c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)))), + all(IsZero(unname(round(BinomDiffCI(84, 101, 89, 105, method = meth), 4)[, -1]) - + cbind(c(-0.1162, -0.1259, -0.1152, -0.116, -0.1188, + -0.1191, -0.1177, -0.1245, -0.1216, -0.1168, -0.117), + c(0.0843, 0.094, 0.0834, 0.0843, 0.0858, 0.086, 0.0851, + 0.0918, 0.0898, 0.085, 0.0852)))) ))