function.txt

# Mode function courtesy of John Verzani
	mode <- function(x) UseMethod("mode")
    	mode.default <- function(x) {
      		tb <- table(x)
      		val <- names(tb[ which(tb == max(tb))])
      	return(val)
    }
    	mode.numeric <- function(x) {
      	val <- NextMethod(x)
      	as.numeric(val)
                                            #  as.numeric(Mode.default(x))
    }
   
   
# Population Variance function by wbw
   	var.pop <- function(x) {
  	tabs <- table(x)
  	nmiss <- length(x) - sum(tabs)
	sigma2 <- var(x, na.rm = TRUE)*((length(x)-nmiss) - 1)/(length(x)- nmiss)
  	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The population variance of the non-missing cases is",sigma2,"\n")
   }

# Function for the average deviation, AD by wbw
  	AD <- function(x) {
   	tabs <- table(x)
	nmiss <- length(x) - sum(tabs)
	dev <- abs(x - mean(x, na.rm = TRUE))
   	AD <- mean(dev, na.rm = TRUE)
   	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The Average Deviation of the non-missing cases is",AD,"\n")
   }
 

# Skewness and Kurtosis functions for Fisher coefficients by wbw for complete data
skewness <-function(x){
	nn <- length(x)
	m3 <- sum((x-mean(x))^3)*nn
	s3 <- sd(x)^3
	g1 <- m3/(s3*(nn-1)*(nn-2))
	g1
  }

SEsk <- function(x) {
	nn <- length(x)
	SEg1 <- sqrt((6*nn*(nn-1))/((nn-2)*(nn+1)*(nn+3)))
	SEg1
  }

kurtosis <-function(x) {
	nn <- length(x)
	m4<-sum((x-mean(x))^4)*nn*(nn+1)
	s4<-var(x)^2
	g2 <- m4/(s4*(nn-1)*(nn-2)*(nn-3)) - (3*(nn-1)^2)/((nn-2)*(nn-3))
	g2
  }

SEku <- function(x) {
	nn <- length(x)
	SEg2 <- sqrt((24*nn*(nn-1)^2)/((nn-3)*(nn-2)*(nn+3)*(nn+5)))
	SEg2
  }

# Skewness and Kurtosis functions for Fisher coefficients by wbw allowing for missing data
skewness.miss <-function(x){
   	tabs <- table(x)
   	nmiss <- length(x) - sum(tabs)
   	nn <- length(x) - nmiss
   	m3 <- sum((x - mean(x, na.rm = TRUE))^3, na.rm = TRUE)*nn
   	s3 <- sd(x, na.rm = TRUE)^3
   	g1 <- m3/(s3*(nn - 1)*(nn - 2))
   	g1
   	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The skew of the non-missing cases is",g1,"\n")
   }

SEsk.miss <- function(x) {
   	tabs <- table(x)
   	nmiss <- length(x) - sum(tabs)
   	nn <- length(x) - nmiss
   	SEg1 <- sqrt((6*nn*(nn - 1))/((nn - 2)*(nn + 1)*(nn + 3)))
   	SEg1
	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The SE of skew of the non-missing cases is",SEg1,"\n")

   }
       
kurtosis.miss <- function(x) {
   	tabs <- table(x)
   	nmiss <- length(x) - sum(tabs)
   	nn <- length(x) - nmiss
   	m4<-sum((x - mean(x, na.rm = TRUE))^4, na.rm = TRUE)*nn*(nn + 1)
   	s4<-var(x, na.rm = TRUE)^2
   	g2 <- m4/(s4*(nn - 1)*(nn - 2)*(nn - 3)) - (3*(nn - 1)^2)/((nn - 2)*(nn - 3))
   	g2
	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The kurtosis of the non-missing cases is",g2,"\n")
   }
       
SEku.miss <- function(x) {
   	tabs <- table(x)
   	nmiss <- length(x) - sum(tabs)
   	nn <- length(x) - nmiss
   	SEg2 <- sqrt((24*nn*(nn - 1)^2)/((nn - 3)*(nn - 2)*(nn + 3)*(nn + 5)))
   	SEg2
	cat("There are",nmiss,"cases with a missing value.","\n")
	cat("The SE of kurtosis of the non-missing cases is",SEg2,"\n")
   }

# Function for index of dispersion written by jf
# First establish a vector of counts of the categories
# Such as vals <- c(  ,  ,  ,  ,  )
indexdisp <- function(x)  {
	c <- length(x)
   	N <- sum(x)
   	nj2 <- sum(x^2)
   	IoD <- (c*((N^2)-nj2))/(N^2*(c-1))
   	IoD
   }

# Syntax for completing the g-squared test for normality by wbw 
g.square <- function(x)  {
	nn <- length(x)
	z1 <- skewness(x)/SEsk(x)
	z2 <- kurtosis(x)/SEku(x)
	g2 <- z1^2 + z2^2
	print ("The value of g-squared")
	print (g2)
	print ("The sample size")
	print (nn)
	g2crit.10 <- 1.950626 + .1556472*log(nn) + 12.4138806*(1/log(nn)) - .4386174*(1/exp(nn/200)) - 9.79946114*(1/sqrt(nn))
 	g2crit.05 <- 14.093477 - .506116*log(nn) - 23.948084*(1/log(nn)) - 56.880456*(1/nn) + 80.242742*(1/sqrt(nn)) - 35.834511*(1/nn^(1/3)) + 187.756861*(1/nn^2)
	g2crit.01 <- -3.741696 + 1.266011*log(nn) - 23.808164*(1/log(nn)) - 92.874186*(1/nn) - 15.247724*(1/sqrt(nn)) + 94.06006*(1/nn^(1/3)) 
	if(g2 >= g2crit.01) print ("Significant at the .01 level") else
	if(g2 >= g2crit.05) print ("Significant at the .05 level") else
	if(g2 >= g2crit.10) print ("Significant at the .10 level") else
	print ("p > .10")
   }  

# Defining the function for the Line Test, lt
line.test <- function(x) {
	n <- length(x)
	rank.x <- rank(x, ties.method = "average")
	perc <- (rank.x - .5)/n
	norm <- qnorm(perc, 0, 1)
	lt <- cor(x, norm)
	cat("The correlation between the actual cdf and normal cdf is", lt,"\n")
	crit.10 <- .996357 -.036317*(1/n)+.003680*(1/exp(n/200))-
	  .432204*(1/sqrt(n))+.160390*(1/n^(1/3))
	crit.05 <- .995539 -.070244*(1/n)+.004649*(1/exp(n/200))-
	  .522094*(1/sqrt(n))+.195069*(1/n^(1/3))
	crit.01 <- .993098 -.102328*(1/n)+.007911*(1/exp(n/200))-
	  .786898*(1/sqrt(n))+.298123*(1/n^(1/3))
	cat("The critical value at the .10 level", crit.10,"\b.\n")
	cat("The critical value at the .05 level", crit.05,"\b.\n")
	cat("The critical value at the .01 level", crit.01,"\b.\n")
	if (lt > crit.10) print("Not significant at the .10 level.")
	if (lt <= crit.10 & lt > crit.05) print("Significant at the .10 level.")
      if (lt <= crit.05 & lt > crit.01) print("Significant at the .05 level.")
      if(lt <= crit.01) print("Significant at the .01 level.")
	
   }

# Effect sizes for oneway ANOVA
oneway.effect.size <- function(dv, group) {
	tab <- table(group)
	f.gp <- factor(group)
	k <- length(tab)
	result <- aov(dv ~ f.gp)
	sum <- anova(result)
	ssag <- sum[1,2]
	sswg <- sum[2,2]
	dfag <- sum[1,1]
	mswg <- sum[2,3]
	sstot <- ssag + sswg
	eta.square <- ssag/sstot	
	omega.square <- (ssag - dfag*mswg)/(mswg + sstot)
	cat("Eta Squared = ",eta.square,"\n")
	cat("Omega-hat Squared = ",omega.square,"\n")
}	
# This assumes that you have read in data.frame with dv and group
# Call the function
#oneway.effect.size(recall, group)


# Function for completing contrasts, planned or Scheffe
# To use the function, read in the data.frame with data <- read.table(...)
# attach(data)
# Create factor level of the grouping variable
# data$f.gp <- factor(group)  #Name must be f.gp
# names(data)
# Supply the coefficients for the contrast.  For example,
# weights <- c(1, -1, 0)
# Call the function where dv is the name of the dependent variable in the data.frame
# mcp(dv, f.gp, weights)
mcp <- function(dv, f.gp, weights)  {
	samsiz <- table(f.gp)
	k <- length(samsiz)
	means <- tapply(dv, f.gp, mean)
	model <- aov(dv ~ f.gp)
	anova.table <- anova(model)
	dfwg <- anova.table[2,1]
	mswg <- anova.table[2,3]
	# Test for length
	if(length(weights) != k) print("Number of coefficients not equal to k")
	# Test for contrast
	if(sum(weights) != 0) print("Warning: Sum of coefficients not equal to zero")
	d <- sum(means*weights)
	denom <- mswg*sum((weights*weights/samsiz))
	f.stat <- d^2/denom
	p.value <- pf(f.stat, 1, dfwg, lower.tail = FALSE)
	cat("The F-statistic for the contrast is", f.stat, "\n")
	cat("The p-value for the contrast is", p.value, "\n")
}


# Effect sizes for oneway ANOVA
oneway.effect.size <- function(dv, group) {
	tab <- table(group)
	f.gp <- factor(group)
	k <- length(tab)
	result <- aov(dv ~ f.gp)
	sum <- anova(result)
	ssag <- sum[1,2]
	sswg <- sum[2,2]
	dfag <- sum[1,1]
	mswg <- sum[2,3]
	sstot <- ssag + sswg
	eta.square <- ssag/sstot	
	omega.square <- (ssag - dfag*mswg)/(mswg + sstot)
	cat("Eta Squared = ",eta.square,"\n")
	cat("Omega-hat Squared = ",omega.square,"\n")
}	
# This assumes that you have read in data.frame with dv and group
# Call the function
#oneway.effect.size(recall, group)

# Function for Univariate Statistics
univ.describe <- function(x) {
	source("c:/r/functions.txt")
		cat(" N = ", length(x), "\n",
		"Mean", mean(x), "\n",
		"Median", median(x), "\n",
		"Mode(s):", mode(x), "\n",
		"Minimum", min(x), "\n",
		"Maximum", max(x), "\n",
		"Range", max(x) - min(x), "\n",
		"Interquartile Range", IQR(x), "\n",
		"Variance", var(x), "\n",
		"Standard Deviation", sd(x), "\n",
		"Coefficient of Variation", 100*sd(x)/mean(x), "\n",
		"Skewness", round(skewness(x), 3), "  St. Err.",
			round(SEsk(x), 3), "\n",
		"Kurtosis", round(kurtosis(x)< 3), "  St. Err.", 
			round(SEku(x), 3), "\n", "\n", "\n",
		"Quantiles", "\n", "\n",
		"q99", quantile(x, .99), "\n",
		"q95", quantile(x, .95), "\n",
		"q90", quantile(x, .90), "\n",
		"q75", quantile(x, .75), "\n",
		"q50", quantile(x, .50), "\n",
		"q25", quantile(x, .25), "\n",
		"q10", quantile(x, .10), "\n",
		"q05", quantile(x, .05), "\n",
		"q01", quantile(x, .01), "\n")
		stem(x)
		boxplot(x)
}