################################################################
#          TTESTS ETC
################################################################

tt <- function (x,m,nl=0) {
	tmp <- t.test (x, mu=m, alternative="two.sided")
	t <- tmp$statistic
	ndt <- numberOfDigits (t)
	if (nl == 0){
		r <- paste(format (tmp$estimate, digits=3), "%+-%",format(sd(x),digits=numberOfDigits (sd(x))),"%, t(",tmp$parameter,")=",format (t, digits=ndt),", p=",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste (r, ", ns", sep= "")
		}
	} else if (nl > 0) {
		r <- paste(format (tmp$estimate, digits=3), "%+-%",format(sd(x),digits=numberOfDigits (sd(x))),"%\nt(",tmp$parameter,")=",format (t, digits=ndt),", p=",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste(r, ", ns", sep="")
		}
	}
	r
}

# one \n more

tt2 <- function (x,m,nl=0) {
	tmp <- t.test (x, mu=m, alternative="two.sided")
	t <- tmp$statistic
	ndt <- numberOfDigits (t)
	if (nl == 0){
		r <- paste(format (tmp$estimate, digits=3), "%+-%",format(sd(x),digits=numberOfDigits (sd(x))),"%, t(",tmp$parameter,")=",format (t, digits=ndt),"\np=",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste (r, ",ns", sep= "")
		}
	} else if (nl > 0) {
		r <- paste(format (tmp$estimate, digits=3), "%+-%",format(sd(x),digits=numberOfDigits (sd(x))),"%\nt(",tmp$parameter,")=",format (t, digits=ndt),"\np=",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste(r, ", ns", sep="")
		}
	}
	r
}

tt3 <- function (x,m,nl=0) {
	tmp <- t.test (x, mu=m, alternative="two.sided")
	t <- tmp$statistic
	ndt <- 1 + numberOfDigits (t) # get 2 decimals instead of 1
	if (nl == 0){
		r <- paste("(M = ", format (tmp$estimate, digits=4), ", SD = ",format(sd(x),digits=numberOfDigits (sd(x))),"), T(",tmp$parameter,") = ",format (t, digits=ndt),", p = ",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste (r, ", ns", sep= "")
		}
	} else if (nl > 0) {
		r <- paste("(M = ", format (tmp$estimate, digits=4), ", SD = ", format(sd(x),digits=numberOfDigits (sd(x))),"), T(",tmp$parameter,") = ",format (t, digits=ndt),", p = ",format (tmp$p.value, digits=3), sep="")
		if (tmp$p.value>0.05){
			r <- paste(r, ", ns", sep="")
		}
	}
	r
}

tt4 <- function (x,m,nl=0) {
	tmp <- t.test (x, mu=m, alternative="two.sided")
	t <- tmp$statistic
        d <- abs (mean (x) - m) / sd (x)
        ci1 <- tmp$"conf.int"[1]
        ci2 <- tmp$"conf.int"[2]
	ndt <- 1 + numberOfDigits (t) # get 2 decimals instead of 1
	if (nl == 0){
		r <- paste("(M = ", format (tmp$estimate, digits=4),
                           ", SD = ", format(sd(x),digits=numberOfDigits (sd(x))),
                           "), T(", tmp$parameter, ") = ", format (t, digits=ndt),
                           ", p = ",format (tmp$p.value, digits=3),
                           ", D = ",format (d, digits=2),
                           ", CI = ",format (ci1, digits=4), ", ", format (ci2, digits=4),
                           sep="")
		if (tmp$p.value>0.05){
			r <- paste (r, ", ns", sep= "")
		}
	} else if (nl > 0) {
		r <- paste("(M = ", format (tmp$estimate, digits=4),
                           ", SD = ", format(sd(x),digits=numberOfDigits (sd(x))),
                           "), T(",tmp$parameter,") = ", format (t, digits=ndt),
                           ", p = ", format (tmp$p.value, digits=3),
                           ", D = ",format (d, digits=2),
                           ", CI = ",format (ci1, digits=4), ", ", format (ci2, digits=4),
                           sep="")
		if (tmp$p.value>0.05){
			r <- paste(r, ", ns", sep="")
		}
	}
	r
}


ff <- function (v1,v2) {
	y <- c (v1, v2)
	#x <- gl (2, length (v1), length (v1) + length (v2))
        x <- factor (c(rep(1, length (v1)), rep (2, length (v2))))
	tmp <- summary (aov (y ~ x))

	f <- tmp[[1]]$"F value"[1]
	ndt <- numberOfDigits (f)
	f <- format (f, digit=ndt)
	p <- format (tmp[[1]]$"Pr(>F)"[1], digit=3)


	r <- paste("F(", tmp[[1]]$Df[1], ",", tmp[[1]]$Df[2], ")=", f, ", p=", p, sep="")
	if (tmp[[1]]$"Pr(>F)"[1] > 0.05) {
		r <- paste (r, ", ns", sep="")
	}
	r
}

ff2 <- function (v1,v2) {
        et <- format (eta_sq (v1, v2), digit=3)
	y <- c (v1, v2)
	#x <- gl (2, length (v1), length (v1) + length (v2))
        x <- factor (c(rep(1, length (v1)), rep (2, length (v2))))
	tmp <- summary (aov (y ~ x))

	f <- tmp[[1]]$"F value"[1]
	ndt <- numberOfDigits (f)
	f <- format (f, digit=ndt)
	p <- format (tmp[[1]]$"Pr(>F)"[1], digit=3)


	r <- paste("F(", tmp[[1]]$Df[1], ",", tmp[[1]]$Df[2], ") = ", f, ", p = ", p, ", et = ", et, sep="")
	if (tmp[[1]]$"Pr(>F)"[1] > 0.05) {
		r <- paste (r, ", ns", sep="")
	}
	r
}

ff2rank <- function (v1,v2) {
    	y <- rank (c (v1, v2))
        et <- format (eta_sq (y[1:(length(v1))], y[(1+length(v1)):length(y)]),
                              digit=3)

	#x <- gl (2, length (v1), length (v1) + length (v2))
        x <- factor (c(rep(1, length (v1)), rep (2, length (v2))))
	tmp <- summary (aov (y ~ x))

	f <- tmp[[1]]$"F value"[1]
	ndt <- numberOfDigits (f)
	f <- format (f, digit=ndt)
	p <- format (tmp[[1]]$"Pr(>F)"[1], digit=3)


	r <- paste("F(", tmp[[1]]$Df[1], ",", tmp[[1]]$Df[2], ") = ", f, ", p = ", p, ", et = ", et, sep="")
	if (tmp[[1]]$"Pr(>F)"[1] > 0.05) {
		r <- paste (r, ", ns", sep="")
	}
	r
}


ffr <- function (d, nBetwSubjCond = 1) {
	# Mixed analysis of variance
	# Arguments data vector d, number of BETWEEN subject conditions
	# Works for one between-subject factor with two levels,
	# and zero or one within subject factors

	nSubj <- length (d)/2

	subj <- factor (rep(1:nSubj, 2))
	if (nBetwSubjCond > 1){
		bs <- c()
		for (i in c(1:nBetwSubjCond)) {
			bs <- c(bs, rep (i, nSubj/nBetwSubjCond) )
		}
		bs <- rep (bs, 2) # since two within conditions
		bs <- factor (bs)
	}
	ws <- factor (c(rep (1,length(d)/2), rep (2, length (d)/2)))

	if (nBetwSubjCond > 1){
		res <- summary (aov (d~bs*ws+Error(subj/ws)))
	} else {
		res <- summary (aov (d~ws+Error(subj/ws)))
	}
	res
}


numberOfDigits <- function (x) {
	# number of digits for t- or F-value
	if (abs(x) < 1){
		ndt <- 1
	} else {
		ndt <- 2 + floor (log10 (abs(x)))
	}
	ndt
}

strip2 <- function (x, y, xName, yName, printTests="y", yLabel="% Correct", xLabel=""){
	if (printTests == "y"){
		xT <- tt (100*x,50,1)
		xT <- paste (xName, "\n", xT)
		yT <- tt (100*y,50,1)
		yT <- paste (yName, "\n", yT)
	} else {
		xT <- xName
		yT <- yName
	}
	f <- ff (x,y)
	f <- paste (xLabel, f)

	par(mar=c(5.1,6,4,7))
	par(mgp=c(4,1.7,0))
	par (cex.lab=1.1, font.lab=8)
	par (cex.axis=.8, font.axis=6)
	par (cex.main=1.5, font.main=7)
	stripchart (c(100*as.data.frame(x), 100*as.data.frame(y)), method="stack", jitter=0.1, offset=0.7, vertical=TRUE,
	group.names=c(xT,yT),
	xlim=c(0.5, 2.5), ylim=c(0,100), ylab=yLabel, xlab=f, pch=16)
	abline (h=50, lty=3)
}


strip3 <- function (y, x=1:dim(y)[2], ylab="% Correct", main="")
{
	if (max (y, na.rm=T) <= 1){
		y <- 100*y
	}
	par(mar=c(4.5,4.5,2.5,0.5))
	par (cex.lab=1.5, font.lab=6, cex.axis=1.2)
	par (bty="o")
	par (yaxt="s")
	par (xaxt="n")
	par (las=1)
	stripchart (y, method="stack", jitter=0.1, offset=0.7, vertical=TRUE,
		xlim=c(0.5, dim(y)[2]+ 0.5),
		at = x,
		ylim=c(0,100),
		ylab=ylab,
		pch=16,
		main = main)
	abline (h=50, lty=3)
	points (-.35+x, mean (y, na.rm=T), pch=5)
	text (-.35+x, -4+mean (y, na.rm=T),
		labels= paste (format (round (mean (y, na.rm=T), digits=3), digits=3), "%", sep = ""),
		cex=1)
	par (xaxt="s")
}

strip4 <- function (dat, x=1:dim(dat)[2],
			main="", ylab="",
			xlab_big_at=c(), xlab_big_line= 2, xlab_big=c(),
			xlab_sma_at=c(), xlab_sma_line= 0, xlab_sma=c(),
			xlab_exp_at=c(), xlab_exp_line=-2, xlab_exp=c()){

	if (max (dat, na.rm=T) <= 1){
		dat <- 100*dat
	}

	m <- format (mean(dat, na.rm=T), digits=numberOfDigits(mean (dat, na.rm=T)))

	par(mar=c(4.5,4.5,2.5,0.5))
	par (cex.lab=1.5, font.lab=6, cex.axis=1.2)
	par (bty="o")
	par (yaxt="s")
	par (xaxt="n")
	par (las=1)
	stripchart (dat, method="stack", jitter=0.1, offset=0.7, vertical=TRUE,
		xlim=c(0.5, x[length(x)] + 0.5), at = x,
		ylim=c(0,100),
		ylab=ylab, pch=16,
		main = main)
	abline (h=50, lty=3)
	par (xaxt="s")

	if (length (xlab_big) > 0){
		axis (side=1, at=xlab_big_at, line= xlab_big_line, lty=0, labels=xlab_big)
	}

	if (length (xlab_sma) > 0){
		par (cex.axis=.9)
		axis (side=1, at=xlab_sma_at, line=xlab_sma_line, lty=0, labels=xlab_sma)
	}
	if (length (xlab_exp) > 0){
		par (cex.axis=.7)
		axis (side=1, at=xlab_exp_at, line=xlab_exp_line, lty=0, labels=xlab_exp)
	}

	points (-.35+x, m, pch=5)
	text (-.35+x,-4+as.numeric (m),
		labels=paste(m, "%"), cex=1)
}


strip4b <- function (dat, x=1:dim(dat)[2],
			main="", ylab="",
                        margins = c(4.5,4.5,2.5,2.5),
			xlab_big_at=c(), xlab_big_line= 2, xlab_big=c(),
			xlab_sma_at=c(), xlab_sma_line= 0, xlab_sma=c(),
			xlab_exp_at=c(), xlab_exp_line=-2, xlab_exp=c()){

	if (max (dat, na.rm=T) <= 1){
		dat <- 100*dat
	}

	m <- format (mean(dat, na.rm=T), digits=numberOfDigits(mean (dat, na.rm=T)))

	#par(mar=c(4.5,4.5,2.5,0.5))
	par(mar=margins)
	par (cex.main=1.5, cex.lab=1.5, font.lab=6, cex.axis=1.5)
	par (bty="o")
	par (yaxt="s")
	par (xaxt="n")
	par (las=1)
	stripchart (dat, method="stack", jitter=0.1, offset=0.7, vertical=TRUE,
		xlim=c(0.5, x[length(x)] + 0.5), at = x,
		ylim=c(0,100),
		ylab=ylab, pch=16,
		main = main)
	abline (h=50, lty=3)
	par (xaxt="s")

	if (length (xlab_big) > 0){
		par (cex.axis=1.5)
		axis (side=1, at=xlab_big_at, line= xlab_big_line, lty=0, labels=xlab_big)
	}

	if (length (xlab_sma) > 0){
		par (cex.axis=.9)
		axis (side=1, at=xlab_sma_at, line=xlab_sma_line, lty=0, labels=xlab_sma)
	}
	if (length (xlab_exp) > 0){
		par (cex.axis=1.5)
		axis (side=1, at=xlab_exp_at, line=xlab_exp_line, lty=0, labels=xlab_exp)
	}

	points (-.35+x, m, pch=5)
	text (-.35+x,-4+as.numeric (m),
		labels=paste(m, "%"), cex=1.2)
}



################################################################
#          CORRELATIONS
################################################################


cor2 <- function (x, y){
# get correlation with p-value
# t = r*sqrt(n - 2)/(1 - r^2)
# F = t^2 = r^2*(n-2)/(1 - r^2) ~ F(1, n-2)
	r <- cor (x,y)
	r2 <- r^2
#	dfr <- length (x)
# 	only approximation
#	F <- r2 * dfr / (1 - r2)
#	p <- pf(F, 1, dfr, lower.tail=FALSE)
#	paste("r = ", r,"; p = ", p, sep="")
#
	dfr <- length (x) - 2
	t <- r / sqrt((1-r^2)/dfr)
	p <- pt(abs(t), dfr, lower.tail=FALSE)
	paste("r = ", r,", t(", dfr, ") = ", t, ", p = ", p, sep="")
}

v2z <- function(v, s="n"){
	if (s == "y"){
		sort ((v - mean (v))/sd(v))
	} else {
		(v - mean (v))/sd(v)
	}
}

################################################################
#          RUN
################################################################

pRun <- function (s, n, p=0.5)
{
# Computes the number of s consecutive successes
# in n trials with success probability p
# Formula from Wolfram Math
if (n < s){
	pR = 0
} else {
	if (n == s)
    		pR = p^s
	else
    		pR = pRun (s, n-1, p) + ((p ^(s+1)) * (1 - pRun (s, n-s, p)))
}
pR
}

################################################################
#          SIGNAL DETECTION THEORY
################################################################

aprime <-function(hit,fa) {
        a<-1/2+((hit-fa)*(1+hit-fa) / (4*hit*(1-fa)))
        b<-1/2-((fa-hit)*(1+fa-hit) / (4*fa*(1-hit)))
        a[fa>hit]<-b[fa>hit]
        a[fa==hit]<-.5
        a
}

bppd <-function(hit,fa) {
        ((1-hit)*(1-fa)-hit*fa) / ((1-hit)*(1-fa)+hit*fa)
}

dprime <- function(hit,fa) {
	qnorm(hit) - qnorm(fa)
}

beta <- function(hit,fa) {
	zhr <- qnorm(hit)
	zfar <- qnorm(fa)
	exp(-zhr*zhr/2+zfar*zfar/2)
}

################################################################
#          CONTINGENCY TABLES
################################################################

ct <- function (t1o1,t1o2,t2o1,t2o2) {
	# 2x2 Contingency table
	# Arguments: (Treatment 1, Outcome 1), (Treatment 1, Outcome 2), (Treatment 2, Outcome 1), (Treatment 2, Outcome 2)

	t1s = t1o1 + t1o2
	t2s = t2o1 + t2o2
	o1s = t1o1 + t2o1
	o2s = t1o2 + t2o2
	total = t1s + t2s

	# Expected outcomes
	t1o1e = o1s * t1s / total
	t1o2e = o2s * t1s / total
	t2o1e = o1s * t2s / total
	t2o2e = o2s * t2s / total

	c = (((t1o1-t1o1e)**2) / t1o1e) + (((t1o2-t1o2e)**2) / t1o2e) + (((t2o1-t2o1e)**2) / t2o1e) + (((t2o2-t2o2e)**2) / t2o2e)
	p = pchisq (c, 1, lower.tail=FALSE)

	r <- paste ("\chi^2(1)=", c, ", p=", p)
	if (p>0.05){
		r <- paste (r, ", ns")
	}
	r
}

################################################################
#          remove outliers
################################################################

remOL <- function (dat, ind, thre)
{
	# Arguments : data, index for selecting cases, threshold in terms of sd
	dat[abs(dat[[ind]]-mean(dat[,ind]))<thre*sd(dat[,ind]),]
}

################################################################
#          add subjects from a gaussian distribution
################################################################

addSubjects <- function (d, n)
{
	c(d, rnorm (n, mean (d), sd(d)))
}

################################################################
#          effect sizes
################################################################

d_cohen <- function (x, y)
{
	d <- sqrt (2) * abs (mean(x-y))/sqrt (var (x) + var (y))
	d
}

eta_sq <- function (x, y)
{
	ss <- anova (lm (c(x, y) ~ factor (rep(c(1,2), c(length(x), length (y))))))$"Sum Sq"
	eta <- ss[1]/sum(ss)
	eta
}

noetherCI <- function (x, alpha=0.05) {
    # confidence intervals for the median according
    # to \citeA{Gibbons1992}, p. 165-166, based on \cite{Noether1967}

    x <- sort (x)
    # compute walsh differences
    w <- c()
    for (i in 1:length (x)){
        for (j in 1:i){
            w <-  c(w, (x[i]+x[j])/2)
        }
    }
    w <- sort (w)

    t <- qsignrank (alpha/2, n=length (x))
    ci <- c(median (x),
            w[t+1],
            w[length(w)-t-1])
    ci
}