A little follow up to my testosterone comment (written in German):
When three variables are correlated to each other, and two of the three correlations are known, the range for the third correlation is restricted according to this formula (Olkin, 1981):



Now comes the new part: here's the graphical representation of that range restriction:




As one can see, one, or both of the two given correlations have to be fairly high to imply a positive third correlation.


Olkin, I. (1981). Range restrictions for product-moment correlation matrices. Psychometrika, 46, 469-472. doi:10.1007/BF02293804

Tags: R, statistics

Recently I started to build a new package for R containing Wilcox’ collection of functions for robust statistics.

Wilcox provides 700+ functions for robust statistics, including:

• robust correlations (e.g. percentage bend correlation)
• robust measures of location and mean differences (e.g. Yuen’s test, trimmed mean test with bootstrapping)
• bootstrapping routines for many indices. Using bootstrapping, asymmetric and non-parametric confidence intervals can be calculated
• etc ...

I’ve simply put Wilcox’ functions into a package, so that the environment isn’t cluttered with 700+ functions. I also started to add some help functions (but until now there are really few).

The package is hosted on R-Forge. You can install the package directly in R by typing:
install.packages("WRS", repos="http://R-Forge.R-project.org")

Tags: R, statistics

Although there is a weighted.mean function in R, so far I couldn't find a implementation of weighted.var and weighted.t.test - here they are (the weighted variance is from Gavin Simpson, found on the R malining list):


# weighted variance, inspired by a function from Gavin Simpson on R-Help
var.wt <- function(x, w, na.rm = FALSE) {
if (na.rm) {
w <- w[i <- !is.na(x)]
x <- x[i]
}
sum.w <- sum(w)
return((sum(w*x^2) * sum.w - sum(w*x)^2) / (sum.w^2 - sum(w^2)))
}



weighted.t.test <- function(x, w, mu, conf.level = 0.95, alternative="two.sided", na.rm=TRUE) {

if(!missing(conf.level) &&
(length(conf.level) != 1 || !is.finite(conf.level) ||
conf.level < 0 || conf.level > 1))
stop("'conf.level' must be a single number between 0 and 1")

if (na.rm) {
w <- w[i <- !is.na(x)]
x <- x[i]
}

# to achieve consistent behavior in loops, return NA-structure in case of complete missings
if (sum(is.na(x)) == length(x)) return(list(estimate=NA, se=NA, conf.int=NA, statistic=NA, df=NA, p.value=NA))

# if only one value is present: this is the best estimate, no significance test provided
if (sum(!is.na(x)) == 1) {
warning("Warning weighted.t.test: only one value provided; this value is returned without test of significance!", call.=FALSE)
return(list(estimate=x[which(!is.na(x))], se=NA, conf.int=NA, statistic=NA, df=NA, p.value=NA))
}

x.w <- weighted.mean(x,w, na.rm=na.rm)
var.w <- var.wt(x,w, na.rm=na.rm)
df <- length(x)-1
t.value <- sqrt(length(x))*((x.w-mu)/sqrt(var.w))
se <- sqrt(var.w)/sqrt(length(x))

if (alternative == "less") {
pval <- pt(t.value, df)
cint <- c(-Inf, x.w + se*qt(conf.level, df) )
}
else if (alternative == "greater") {
pval <- pt(t.value, df, lower.tail = FALSE)
cint <- c(x.w - se * qt(conf.level, df), Inf)
}
else {
pval <- 2 * pt(-abs(t.value), df)
alpha <- 1 - conf.level
cint <- x.w + se*qt(1 - alpha/2, df)*c(-1,1)
}

names(t.value) <- "t"
return(list(estimate=x.w, se=se, conf.int=cint, statistic=t.value, df=df, p.value=pval))
}


Tags: R, statistics