# my function to plot ecdf
mycdf<-function(x) {
       x<-sort(na.omit(x))
       y<-(1:length(x)-0.5)/length(x)
       list(x=x,y=y)
       }
       
myround<-function(x) {
         y<-round(x)
         y[!(x+0.5)%%1]<-round(x[!(x+0.5)%%1]+0.1)
         y
         }

# The follofing functions are defined for HD-modeling purposes
# Transformation for natural spline regression
# Eq 2.25 in Harrells book, p. 20
# t is a vector of knots, j the component to be calculated
natural.spline.comp<-function(X,t,j) {
                     k<-length(t)
                     pmax(0,X-t[j])^3-pmax(0,X-t[k-1])^3*(t[k]-t[j])/(t[k]-t[k-1])+
                     pmax(0,X-t[k])^3*(t[k-1]-t[j])/(t[k]-t[k-1])
                     }

# This is a handy function for looking for trends in residuals etc.
mywhiskers<-function(x,y,nclass=10,limits=NA,add=FALSE,se=TRUE,main="",xlab="x",ylab="y",ylim=NA,lwd=1) {
            away<-is.na(x+y)
            x<-x[!away]
            y<-y[!away]
            if (is.na(limits[1])) limits<-seq(min(x),max(x)+1e-10,length=nclass+1) else nclass=length(limits)-1
            means<-sapply(1:nclass,function(i) mean(y[x>=limits[i]&x<limits[i+1]]))
            if (se) { 
                    ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]])/
                                                         sqrt(sum(x>=limits[i]&x<limits[i+1])))
                    } else {
                    ses<-sapply(1:nclass,function(i) sd(y[x>=limits[i]&x<limits[i+1]]))
                    }
            lb<-means-1.96*ses
            ub<-means+1.96*ses
            xclass<-1/2*(limits[-1]+limits[-nclass-1])
            if (add) {
               points(xclass,means) 
               } else {
               if(is.na(ylim[1])) ylim<-c(min(lb),max(ub))
               plot(xclass,means,ylim=ylim,main=main,xlab=xlab,ylab=ylab,xlim=range(x))
               }
            sapply(1:nclass,function(i) lines(xclass[c(i,i)],c(lb[i],ub[i]),lwd=lwd))
            }

# Mean squared error 
mse<-function(y) {
     sum(y^2)/length(y)
     }

# Root mean squared error 
rmse<-function(y) {
     sqrt(sum(y^2)/length(y))
     }

# This is for plotting the response against age and connecting the observations of the same subject
# in a longitudinal study.
plot.linesbygroups<-function(x,y,group,xlab="x",ylab="y",main="",cex=0.5,pch=19,col=1,col.lin=1,
                             lw=FALSE,ylim=NULL,xlim=NULL,add=FALSE,lty="solid") {
    col=rep(col,length(x))
    col.lin=rep(col.lin,length(x))
    if (!add) {
    plot(x,y,type="n",xlab=xlab,ylab=ylab,main=main,cex=cex,pch=pch,col=col,ylim=ylim,xlim=xlim)
    }
    apu<-unique(group)
    for (i in 1:length(apu)) {
        lines(x[group==apu[i]],y[group==apu[i]],col=col.lin[group==apu[i]],lty=lty)
        }
    points(x,y,col=col,cex=cex)
    if (lw) {
       apu<-unique(col)
       for (i in 1:length(apu))
           lines(lowess(x[col==apu[i]],y[col==apu[i]],f=0.5,iter=5),col=apu[i],lwd=3)
       }
    }

# searches the root of equation fn(x)==0 between u and b using a simple up-and-down algorithm    
updown<-function(l,u,fn,crit=6) {
        fnu<-fn(u)
        fnl<-fn(l)
        if (fnl*fnu>0) return(list(par=NA,value=NA))
        value<-fn((u+l)/2)
        while (round(value,crit)!=0) {
              if (fnu*value>0) {
                 u<-(u+l)/2 
                 fnu<-value 
                 } else {
                 l<-(u+l)/2
                 fnl<-value
                 }
              value<-fn((u+l)/2)             
#              cat(l,u,value,"\n")
              }
        list(par=(l+u)/2,value=value)
        }

# a handy function for plotting residuals
resplot<-function(x,res,limits=NA,ylim=range(res),xlab="x",ylab="res",main="",points=TRUE) {
         if (points) {
            plot(x,res,col="red",ylim=ylim,xlab=xlab,ylab=ylab,main=main)
            mywhiskers(x,res,se=FALSE,add=TRUE,lwd=1,limits=limits)
            } else {
            mywhiskers(x,res,se=FALSE,add=FALSE,lwd=1,limits=limits,xlab=xlab,ylab=ylab)
            }
         mywhiskers(x,res,add=TRUE,lwd=5,limits=limits)         
         abline(0,0,col="gray")
         }

# this function forms a symmetric (dim by dim) matrix from vector x that includes only the lower triangle of the matrix.
# Useful when the output of sas mixed is used in R.  
 magicmatrix<-function(x,dim) {
              lengths<-1:dim
              start<-c(0,sapply(1:(dim-1),function(i) sum(lengths[1:i])))+1
              end<-sapply(1:dim,function(i) sum(lengths[1:i]))
              mat<-c()
              for (i in 1:dim) {
                  mat<-rbind(mat,c(x[start[i]:end[i]],rep(NA,dim-lengths[i])))
                  }
              for (i in 1:(dim-1)) for (j in (i+1):dim) mat[i,j]<-mat[j,i]
              mat
              }

         

# searches the root of equation fn(x)==0 between u and b using a simple up-and-down algorithm 
# stops when either the window width or value goes to 0, useful for functions with discontinuities 
# or very large slopes
updown2<-function(l,u,fn,crit=10) {
        fnu<-fn(u)
        fnl<-fn(l)
        if (fnl*fnu>0) return(list(par=NA,value=NA))
        value<-fn((u+l)/2)
        while (round(value,crit)!=0&round(l-u,crit)!=0) {
              if (fnu*value>0) {
                 u<-(u+l)/2 
                 fnu<-value 
                 } else {
                 l<-(u+l)/2
                 fnl<-value
                 }
              value<-fn((u+l)/2)             
#              cat(l,u,value,"\n")
              }
        list(par=(l+u)/2,value=value)
        }