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"))}


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")}

plot.gamma<-function(mu,nu)
     {x<-seq(0,5,.1)*mu
      plot(x,dgamma(x,shape=nu, scale=mu/nu),type="l")}



lm.estimates<-function(formula)
       {results<-lm(formula)
        round(summary(results)$coef,3)} 

lm.compare<-function(formula1,formula2)
     {results1<-lm(formula1)
      results2<-lm(formula2)
      df1<-results1$df.residual
      df2<-results2$df.residual
      df<- df1-df2
      res1<-results1$residuals   
      res2<-results2$residuals 
      n<-length(res1)
      var1<-mean(res1**2)
      var2<-mean(res2**2)
      F<- (var1-var2)*df2/var2/df
      chisq<-n*log(var1/var2)
      df<-df1-df2
      p.value<-1-pf(F,df,df2)
      p.value.appr<-1-pchisq(chisq,df) 
      list(var1=var1,var2=var2,F=F, chisq=chisq,s=df,p.value=p.value,p.value.appr=p.value.appr)
      }   

dev.compare<-function(fit1,fit2)
     {
      df1<-fit1$df.residual
      df2<-fit2$df.residual
      dev1<-fit1$deviance   
      dev2<-fit2$deviance 
      chisq<-dev1-dev2
      df<-df1-df2
      p.value<-1-pchisq(chisq,df) 
      list(chisq=chisq,df=df,p.value=p.value)
      }   



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)} 

OR.CI<-function(formula)
       {results<-glm(formula,family=binomial(link=logit))
        A<-summary(results)$coef
        OR.EST<- exp(A[,1])
        OR.LOWER<- exp(A[,1]-1.96*A[,2])
        OR.UPPER<- exp(A[,1]+1.96*A[,2])
        round(cbind(OR.EST,OR.LOWER,OR.UPPER),4)}
 
OR.estimates<-function(formula, weight=n)
    {results<-glm(formula,family=binomial(link=logit), weight=n)
        A<-summary(results)$coef
        OR.EST<- exp(A[,1])
        OR.LOWER<- exp(A[,1]-1.96*A[,2])
        OR.UPPER<- exp(A[,1]+1.96*A[,2])
        round(cbind(OR.EST,OR.LOWER,OR.UPPER),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)}

RR.estimates<-function(fit)
        {A<-summary(fit)$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)}