.packageName <- "SMIR"
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/NPL.bands.R"
NPL.bands <- function(x,conf.level=c(0.95,0.99)){
# Function to compute nonparametric likelihood confidence bands
# as described by Owen, A. (1997) JASA 90:516-521.
# Requires a vector of numeric values 
# and a confidence level (either 0.95 or 0.99).
# The function stores the limits in a list of three variables
# lower and upper along with the distinct values  of x.
  if (!is.numeric(x)) stop("argument must be numeric")
  yn <- table(x) # table of frequencies of distinct value (sorted)
  yi <- as.numeric(names(yn)) # distinct values
  cn <- as.numeric(cumsum(yn)) # cumulated counts
  nn <- rep(sum(yn)+1,length(cn))
  p <- as.numeric(cn/nn)
  if (conf.level==0.95)
    {
      lambda <- ifelse(nn<=100,(3.0123+0.4835*log(nn)-0.00957*(log(nn))^2
                         -0.001488*(log(nn))^3),
                       (3.0806+0.4894*log(nn)-0.02086*(log(nn))^2))}
  else {if (conf.level==0.99)
          {lambda <- ifelse(nn<=100,(-4.626-0.541*log(nn)+0.0242*(log(nn))^2),
                            (-4.71-0.512*log(nn)+0.0219*(log(nn))^2))}
  else stop("Must be either 0.95 or 0.99")}
  lambda <- sqrt(2*lambda)
  phi <- pbeta(p,1/3,1/3)
  se <- 1/(5.3*sqrt(nn*(p*(1-p))^(1/3)))
  phiu <- phi + lambda*se
  phiu <- ifelse(phiu>1,1,phiu)
  phil <- phi - lambda*se
  phil <- ifelse(phil<0,0,phil)
  pu <- qbeta(phiu,1/3,1/3)
  pl <- qbeta(phil,1/3,1/3)
  list(x=yi,lower=pl,upper=pu)
}

  
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/R2.R"
R2 <- function(model) {
    if (!any(class(model)=="lm")) stop("model must be class ``lm'' or ``glm''")
    yhat <- fitted(model)
    ehat <- residuals(model)
    y <- yhat + ehat
    cor(yhat,y)^2
}  
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/R2CV.R"
R2CV <- function(model){
  if (!any(class(model)=="lm")) stop("model must be class ``lm'' or ``glm''")
  yhat <- fitted(model)
  ehat <- residuals(model)
  y <- yhat + ehat
  hi <- influence(model)$hat
  yhatcve <- (yhat-hi*y)/(1-hi)
  cor(y,yhatcve)^2
}
  
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/Rsq.R"
Rsq <- function(model) cor(model$y,fitted(model))^2

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/coxph.disparity.R"
## but need to account for smallest observations being censored!!!
coxph.disparity<-function(fit){
  if (class(fit)!="coxph") stop("Argument needs to be class ``coxph''")
  require(survival)
bhz<-basehaz(fit,centered=TRUE)
bhz$haz<-diff(c(0,bhz$hazard))
bhz$tint<-diff(c(0,bhz$time))
ti <- apply(outer(fit$y[,1],bhz$time,">="),1,sum)
  -2*(sum(fit$y[,2]*(ti!=0)*(log(bhz$haz[ti+(ti==0)]/bhz$tint[ti+(ti==0)])+fit$linear)-(ti!=0)*bhz$hazard[ti+(ti==0)]*exp(fit$linear)))
}
#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/disparity.R"
disparity <- function(model)
  UseMethod("disparity")
disparity.lm <- function(model){
  e <- resid(model)
  n <- length(e)
  c <- n*(1+log(2*pi)-log(n))
  RSS <- sum(e^2)
  -2*(c+n*log(RSS))
}
disparity.glm <- function(model){
  family <- model$family$family
  logL <- switch(family,
                 "gaussian"=function(model){
                   sigma <- sqrt(summary(model)$dispersion)
                   e <- resid(model)
                   sum(dnorm(e,0,sigma,log=TRUE))},
                 "poisson"=function(model){
                   lambda <- fitted(model)
                   sum(dnorm(model$y,lambda,log=TRUE))},
                 "binomial"=function(model){
                   mu <- fitted(model)
                   size <- model$prior.weights
                   x <- model$y*size
                   sum(dbinom(x,size,mu,log=TRUE))},
                 "Gamma"=function(model){
                   library(MASS)
                   shape <- gamma.shape(model)$alpha
                   sum(dgamma(model$y,shape=shape,
                          scale=fitted(model,
                            type="response"),log=TRUE))}
                 )
  -2*logL(model)}

#line 1 "/private/var/folders/JA/JAolSeAcFueOWF1Uy33XRk+++TI/-Tmp-/Rtmpsps4Rg/R.INSTALL10d63af1/SMIR/R/treg.R"
treg <- function(lm.object, r, verbose=TRUE){
    if (class(lm.object)!="lm") stop("model must be class ``lm''")
    X <- model.matrix(lm.object)
  y <- lm.object$model$y
  nu <- length(y)
  w <- rep(1,nu)
  d <- 0
  convergence.criteria <- 0.0001
  converged <- FALSE
  fit <- lm.object
  while (!converged) {
    rss <- sum(w*resid(fit)^2)
    phi <- (r + 1) * rss/nu
    sigma <- sqrt(phi/r)
    t <- resid(fit)/sqrt(phi)
    w <- 1/(1+t^2)
    e <- -2*(nu*( lgamma((r+1)/2) -0.5*log(pi) - lgamma(r/2) -
                0.5*log(phi)) + (r + 1)/2*sum(log(w)))
    if (verbose)cat("-2 log Lmax =",e," with sigma = ",sigma," and scale parameter ",r,"\n")
    converged <- abs(d-e) < convergence.criteria
    d <- e
    fit <- lm.wfit(X,y,w)
#    cat(r,"\n")
  }
      fit <- c(lm.object,  list(weights = w, disparity=e, tcoef = coef(fit), r=r, sigma = sigma))
  class(fit) <- c("lm","treg")
  fit
  }
#
summary.treg <- function(object, ...)
#  summary.lm
#function (object, correlation = FALSE, symbolic.cor = FALSE,
#    ...)
{
    z <- object
    p <- z$rank
    Qr <- object$qr
#    if (is.null(z$terms) || is.null(Qr))
#        stop("invalid 'lm' object:  no 'terms' nor 'qr' component")
    n <- NROW(Qr$qr)
    rdf <- n - p
    if (is.na(z$df.residual) || rdf != z$df.residual)
        warning("residual degrees of freedom in object suggest this is not an \"lm\" fit")
    p1 <- 1:p
    r <- z$residuals
    f <- z$fitted
    w <- z$weights
        mss <- if (attr(z$terms, "intercept")) {
            m <- sum(w * f/sum(w))
            sum(w * (f - m)^2)
        }
        else sum(w * f^2)
        rss <- sum(w * r^2)
        r <- sqrt(w) * r
#
    resvar <- rss/rdf
    R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    se <- sqrt(diag(R) * resvar)
    est <- z$coefficients[Qr$pivot[p1]]
    tval <- est/se
    ans <- z[c("call", "terms")]
    ans$residuals <- r
    ans$coefficients <- cbind(est, se, tval, 2 * pt(abs(tval),
        rdf, lower.tail = FALSE))
    dimnames(ans$coefficients) <- list(names(z$coefficients)[Qr$pivot[p1]],
        c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
    ans$aliased <- is.na(coef(object))
    ans$sigma <- sqrt(resvar)
    ans$df <- c(p, rdf, NCOL(Qr$qr))
    if (p != attr(z$terms, "intercept")) {
        df.int <- if (attr(z$terms, "intercept"))
            1
        else 0
        ans$r.squared <- mss/(mss + rss)
        ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n -
            df.int)/rdf)
        ans$fstatistic <- c(value = (mss/(p - df.int))/resvar,
            numdf = p - df.int, dendf = rdf)
    }
    else ans$r.squared <- ans$adj.r.squared <- 0
    ans$cov.unscaled <- R
    dimnames(ans$cov.unscaled) <- dimnames(ans$coefficients)[c(1,
        1)]
#    if (correlation) {
#        ans$correlation <- (R * resvar)/outer(se, se)
#        dimnames(ans$correlation) <- dimnames(ans$cov.unscaled)
#        ans$symbolic.cor <- symbolic.cor
#    }
    print(ans$call)
    print(ans$coefficients)
    cat("Disparity = ",signif(z$disparity,5),"\n")
    cat("r value   = ", z$r,"\n")
    if (!is.null(z$na.action))
        ans$na.action <- z$na.action
    class(ans) <- c("summary.lm","summary.treg")
    ans
}

