stand<-function(x) {(x-mean(x))/sd(x)}

Q1<-function(x)
   {quantile(x,0.25)}

Q2<-function(x)
   {quantile(x,0.50)}

Q3<-function(x)
   {quantile(x,0.75)}

b1<-function(x) 
   {n<-length(x)
    z<-stand(x)
    b1<-mean(z**3)
    p<-pchisq(1,n*b1*b1/6)
    list(sqrt_b1=b1,p=p)}


b2<-function(x) 
   {n<-length(x)
    z<-stand(x)
    b2<-mean(z**4)
    p<-pchisq(1,n*(b2-3)**2/24)
    list(b2=b2, p=p)}
 
hist.with.normal<-function(x)
#    Dalgaard(2003)
    {m<-mean(x)
     s<-sd(x)
     hist(x, freq=F)
     curve(dnorm(x,m,s), add=T)}

boxplot.with.normal<-function(x)
     {m<-mean(x)
     s<-sd(x)
     n<-length(x)  
     x.new<-c(x,rnorm(n,m,s))
     w<-c(rep(1,n),rep(2,n))
     boxplot(x.new ~ w, names=c("real data","normal data"))}


    

table2<-function(x,y)
     {T<-table(x,y)
      TT<-cbind(T,apply(T,1,sum))
      TTT<-rbind(TT,apply(TT,2,sum))
      TTTT=round(100*TTT/TTT[,ncol(TTT)],2)
      list(freq=TTT,pcent=TTTT)}


plot.binom<-function(n,p)
     {x<-0:n
     plot(x,dbinom(x,size=n,prob=p),type="h")}

plot.pois<-function(lambda)
     {x<-0:(3*lambda)
     plot(x,dpois(x,lambda),type="h")}

plot.normal<-function(mu,sigma)
     {x<-sigma*seq(-4,4,.1)+mu
      plot(x,dnorm(x,mu,sigma),type="l")}


plot.exponential<-function(mu)
     {x<-seq(0,5,.1)*mu
      plot(x,dexp(x,1/mu),type="l")}

sample.of.means<-function(nn,n)
#     nn = number of samples
#     n = sample size
     {data<-matrix(rnorm(nn*n),nrow=nn)
      apply(data,1,mean)} 

sample.of.medians<-function(nn,n)
#     nn = number of samples
#     n = sample size
     {data<-matrix(rnorm(nn*n),nrow=nn)
      apply(data,1,median)} 


p.test<-function(x,p=0.5)
#      one proportion: 
#      x is a random sample from Ber(p)  
      {n<-length(x)
       est<-mean(x) 
       se.est<-sqrt(est*(1-est)/n)
       CI.lower<-est-1.96*se.est
       CI.upper<-est+1.96*se.est
       z.stat<-(est-p)/se.est
       p.value<-2*min(pnorm(z.stat),1-pnorm(z.stat))
       round(cbind(est,se.est,CI.lower,CI.upper,z.stat,p.value),3)}

p.test2<-function(x,p=0.5)
#      one proportion
#      x is a random sample from Ber(p) 
#      the test and CI obtained utilizing  logit transformation
      {n<-length(x)
       est<-mean(x) 
       logit.est<-log(est/(1-est))
       logit.p<-log(p/(1-p))
       SE<-1/sqrt(n*est*(1-est))
       z.stat<-(logit.est-logit.p)/SE
       p.value<-2*min(pnorm(z.stat),1-pnorm(z.stat)) 
       lower<-logit.est-1.96*SE
       upper<-logit.est+1.96*SE
       CI.lower<- 1/(1+exp(-lower))
       CI.upper<- 1/(1+exp(-upper))
       round(cbind(est, CI.lower,CI.upper,z.stat,p.value),3)}


RD<-function(twobytwo)
       {n11<-twobytwo[1,1]
        n12<-twobytwo[1,2] 
        n21<-twobytwo[2,1] 
        n22<-twobytwo[2,2]
        n1<-n11+n12
        n2<-n21+n22
        p1<-n12/n1
        p2<-n22/n2
        var1<-p1*(1-p1)/n1
        var2<-p2*(1-p2)/n2
        RD<-p2-p1
        SE<-sqrt(var1+var2)
        lower<-RD-1.96*SE
        upper<-RD+1.96*SE 
        z<-RD/SE
        p.value<-2*min(pnorm(z),1-pnorm(z))
        list(p1.EST=p1, p2.EST=p2, RD.EST=RD, RD.CI=c(lower,upper), p.value=p.value)}
 
RD1<-function(x,y)
        {RD(table(x,y))}

         
RR<-function(twobytwo)
       {n11<-twobytwo[1,1]
        n12<-twobytwo[1,2] 
        n21<-twobytwo[2,1] 
        n22<-twobytwo[2,2]
        n1<-n11+n12
        n2<-n21+n22
        p1<-n12/n1
        p2<-n22/n2
        logp1<-log(p1)
        logp2<-log(p2)
        var1<-(1-p1)/p1/n1
        var2<-(1-p2)/p2/n2
        logRR<-logp2-logp1
        SE<-sqrt(var1+var2)
        lower<-logRR-1.96*SE
        upper<-logRR+1.96*SE 
        z<-logRR/SE
        p.value<-2*min(pnorm(z),1-pnorm(z))
        list(p1.EST=p1, p2.EST=p2, RR.EST=exp(logRR), RR.CI=c(exp(lower),exp(upper)), p.value=p.value)}
 
RR1<-function(x,y)
        {RR(table(x,y))}         

OR<-function(twobytwo)
       {n11<-twobytwo[1,1]
        n12<-twobytwo[1,2] 
        n21<-twobytwo[2,1] 
        n22<-twobytwo[2,2]
        n1<-n11+n12
        n2<-n21+n22
        p1<-n12/n1
        p2<-n22/n2
        logitp1<-log(p1/(1-p1))
        logitp2<-log(p2/(1-p2))
        var1<-1/(1-p1)/p1/n1
        var2<-1/(1-p2)/p2/n2
        logOR<-logitp2-logitp1
        SE<-sqrt(var1+var2)
        lower<-logOR-1.96*SE
        upper<-logOR+1.96*SE 
        z<-logOR/SE
        p.value<-2* min(pnorm(z),1-pnorm(z))
        list(p1.EST=p1, p2.EST=p2, OR.EST=exp(logOR),
         OR.CI=c(exp(lower),exp(upper)), p.value=p.value)}
          
OR1<-function(x,y)
        {OR(table(x,y))}


logistic<-function(formula)
       {results<-glm(formula,family=binomial(link=logit))
        round(summary(results)$coef,4)} 


 

RR.CI<-function(formula)
       {results<-glm(formula,family=binomial(link=log))
        A<-summary(results)$coef
        RR.EST<- exp(A[,1])
        RR.LOWER<- exp(A[,1]-1.96*A[,2])
        RR.UPPER<- exp(A[,1]+1.96*A[,2])
        round(cbind(RR.EST,RR.LOWER,RR.UPPER),4)}

normal.model<-function(formula)
       {A<-summary(lm(formula))$coefficients
        est<- A[,1]
        se.est<- A[,2]
        CI.lower<- A[,1]-1.96*A[,2]
        CI.upper<- A[,1]+1.96*A[,2]
        t.stat<-est/se.est
        p.value<- A[,4]
        round(cbind(est,se.est,CI.lower,CI.upper,t.stat,p.value),3)}

normal.model.residuals<-function(formula)
       {res<-summary(lm(formula))$residuals
        res}

hist.residuals<-function(formula)
        {res<- normal.model.residuals(formula)
         hist.with.normal(res)}
        

median.test<-function(x,mu=0)
      {n<-length(x)
       s<-sum(x>mu)
       sign.test<-s
       p.value<-2*min(pbinom(s,n,0.5),1-pbinom(s-1,n,0.5),0.5)
       i<-qbinom(0.025,n,0.5)
       x<-sort(x)
       conf.level<-1-2*pbinom(i-1,n,0.5) 
       median<-median(x)
       CI.lower<-x[i]
       CI.upper<-x[n+1-i]
       round(cbind(median,CI.lower,CI.upper,conf.level,sign.test,p.value),3)}
       

logistic.model<-function(formula)
       {results<-glm(formula,family=binomial(link=logit))
        A<-summary(results)$coef
        OR.est<- exp(A[,1])
        CI.lower<- exp(A[,1]-1.96*A[,2])
        CI.upper<- exp(A[,1]+1.96*A[,2])
        z.stat<-A[,3]
        p.value<-A[,4] 
        round(cbind(OR.est,CI.lower,CI.upper,z.stat,p.value),3)}


OR.model<-function(formula)
       {results<-glm(formula,family=binomial(link=logit))
        A<-summary(results)$coef
        OR.est<- exp(A[,1])
        CI.lower<- exp(A[,1]-1.96*A[,2])
        CI.upper<- exp(A[,1]+1.96*A[,2])
        z.stat<-A[,3]
        p.value<-A[,4] 
        round(cbind(OR.est,CI.lower,CI.upper,z.stat,p.value),3)}


library(survival)

cox.model<-function(formula)
      {results<-coxph(formula)
        A<-summary(results)$coef
        HR.est<- exp(A[,1])
        HR.lower<- exp(A[,1]-1.96*A[,3])
        HR.upper<- exp(A[,1]+1.96*A[,3])
        z.stat<-A[,4]
        p.value<-A[,5] 
        round(cbind(HR.est,HR.lower,HR.upper,z.stat,p.value),3)}

survival.curves<-function(formula)        
       {results <- survfit(formula)
        plot(results)}