#Mike's Package
library(boot)
library(psych)


#BOOTD
#Boot d requires as input a simple dataset in the form of a continuous DV and Grouping variable

bootd <- function(data, i){
    datadiff <- data
    split1=data.frame(split(datadiff [i,], datadiff [,2])[1])
    split2=data.frame(split(datadiff [i,], datadiff [,2])[2])
    g1=split1[,1]
    g2=split2[,1]
    n1=length(g1)
    n2=length(g2)
    df=n1+n2-2
  meandiff= mean(g1)-mean(g2)
  poolsd <- sqrt((sd(g1)^2 * (n1 - 1) + sd(g2)^2 * (n2 - 1))/(df))
  Cohend = meandiff/poolsd
 }


#BOOTCOR
#Bootstrap correlation.  Will work on first two variables of a dataset.
bootcor <- function(data, i) cor(data[i, 1], data[i, 2], use="complete.obs")


#BOOTREGCOEFF. With the appropriate libraries loaded, one can change lm to lmRob (robust library), lmrob (robustbase), 
#or rlm (MASS)
bootregcoeff <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data=d)
  return(coef(fit))
}

#BOOTRSQ
bootRsq <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data=d)
  return(summary(fit)$r.square)
}

#VECTDIFF
#Vector Difference Norm.  Only works with a complete data matrix for two samples.  Each column is a different sample.
vectdiff= function(data,indices){
x=data[indices,1]
y=data[indices,2]
out=t(y)%*%y-2*t(y)%*%x+t(x)%*%x
return(out)
}

#Harmonic Mean with NAs for row and column.  When doing rows, creates a column of the means. When columns, the last row are the harmonic means.
hmeanrow=function(vars){
harmonic.mean=1/rowMeans(1/vars,na.rm=T)
data.frame(vars,harmonic.mean)
}
hmeancol=function(vars){
harmonic.mean=1/colMeans(1/vars,na.rm=T)
data.frame(rbind(vars,harmonic.mean))
}

#INF.CONF.INT
#Tryon's inferential confidence interval with equivalence test.
inf.conf.int=function(g1,g2,conf.level=.95,paired=F){
      if(paired==FALSE){g1<-na.omit(g1); g2<-na.omit(g2)}
      if(paired==TRUE){dataomit=na.omit(data.frame(g1,g2));g1=dataomit[,1] ; g2=dataomit[,2]}
      m1=mean(g1)
      m2=mean(g2)
      s1=sd(g1)
      s2=sd(g2)
      n1=length(g1)
      n2=length(g2)
      se1=s1/sqrt(n1)
      se2=s2/sqrt(n2)
      df1=n1-1
      df2=n2-1
      tcrit1=qt(1-(1-conf.level)/2,df1)
      tcrit2=qt(1-(1-conf.level)/2,df2)
      adjust=  sqrt(se1^2+se2^2)/(se1+se2)
        if(paired==TRUE) {adjust=sqrt(se1^2+se2^2-2*cor(g1,g2,use="complete.obs")*se1*se2)/(se1+se2)}
        if(paired==FALSE){pooledsd= sqrt((df1*var(g1) + df2*var(g2))/(df1+df2))}
      margerr1=tcrit1*adjust*se1
      margerr2=tcrit2*adjust*se2
      LL1=m1-margerr1
      UL1=m1+margerr1
      LL2=m2-margerr2
      UL2=m2+margerr2
      if(paired==FALSE){
            if(m1<=m2){
              deltalowm=m1
              deltamargerrlo=margerr1
              deltamargerrhi=margerr2}
            if(m1>=m2){
              deltalowm=m2
              deltamargerrlo=margerr2
              deltamargerrhi=margerr2}
            deltahim=deltalowm+(pooledsd*.2)
            deltaLL=deltalowm-deltamargerrlo
            deltaUL=deltahim+deltamargerrhi
            deltarange=data.frame(LL=deltaLL,UL=deltaUL,Range=deltaUL-deltaLL)
            }
      if(m1>m2){overlap=LL1-UL2}
      if(m1<m2){overlap=LL2-UL1}
if(m1>m2){equirange=data.frame(LL2,UL1,Range=UL1-LL2)}
if(m1<m2){equirange=data.frame(LL1,UL2,Range=UL2-LL1)}
if(overlap>0){note="The intervals do not overlap. The samples are statistically different at the specified alpha level."}
if(overlap<0){note="The intervals overlap. The samples are not statistically different at the specified alpha level."}
if(paired==F){
if(equirange[,3]<=deltarange[,3]){equinote="The samples are statistically equivalent. The observed range is less than one based on a small effect of Cohen's d =.20."}
if(equirange[,3]>deltarange[,3]){equinote="Equivalence not established.  The observed range is greater than one based on a small effect of Cohen's d =.20."}
  return(list("Lower Limit Group 1"=LL1,"Upper Limit Group 1"=UL1,"Lower Limit Group 2"=LL2,"Upper Limit Group 2"=UL2, "Difference Outcome"=note,"Equivalence Range"=equirange, "Delta Range"=deltarange,"Equivalence Outcome"=equinote))
}
if(paired==T){
  return(list("Lower Limit Group 1"=LL1,"Upper Limit Group 1"=UL1,"Lower Limit Group 2"=LL2,"Upper Limit Group 2"=UL2, "Outcome"=note,"Equivalence Range"=equirange))}
}

#WINVAR
# Copied from Wilcox package. tr is the amount of Winsorization which defaults to .2.

winvar<-function(x,tr=.2,na.rm=F){
   if(na.rm)x<-x[!is.na(x)]
   y<-sort(x)
   n<-length(x)
   ibot<-floor(tr*n)+1
   itop<-length(x)-ibot+1
   xbot<-y[ibot]
   xtop<-y[itop]
   y<-ifelse(y<=xbot,xbot,y)
   y<-ifelse(y>=xtop,xtop,y)
   winvar<-var(y)
   winvar
}

#ROBUSTD.  For independent samples. See Algina, Keselman, and Penfield 2005.
robustd= function(g1,g2,trim=.2){
  g1<-g1[!is.na(g1)]  # Remove missing values
  g2<-g2[!is.na(g2)] 
  m1=mean(g1,tr=trim)
  m2=mean(g2,tr=trim)
  n1=length(g1)
  n2=length(g2)
  df1=n1-1
  df2=n2-1
  var1=winvar(g1)
  var2=winvar(g2)
    pooledwinsd= (sqrt((df1*var1 + df2*var2)/(df1+df2)))
    robustd= (m1-m2)/pooledwinsd
  return(list("Robust Cohen's d"=robustd))
  }