.packageName <- "betareg"
"betareg" <- function (formula, link = "logit", data)
{
  call <- match.call()
  m <- match.call(expand = FALSE)
  link = c(link)
  m$link <- NULL
  m$x <- m$y
  m[[1]] <- as.name("model.frame")
  m <- eval(m, sys.frame(sys.parent()))
  Terms <- attr(m, "terms")
  Y <- model.extract(m, response)
  X <- model.matrix(Terms, m, contrasts)
  if (min(Y) <= 0 || max(Y) >= 1)
    stop("OUT OF RANGE (0,1)!")
  offset <- model.offset(m)
  linktemp <- substitute(link)
  if (!is.character(linktemp)) {
    linktemp <- deparse(linktemp)
    if (linktemp == "link")
      linktemp <- eval(link)
  }
  if (any(linktemp == c("logit", "probit", "cloglog")))
    stats <- make.link(linktemp)
  else stop(paste(linktemp, "link not available, available links are \"logit\", ",
                  "\"probit\" and \"cloglog\""))
  link1 <- structure(list(link = linktemp, linkfun = stats$linkfun,
                          linkinv = stats$linkinv, mu.eta = stats$mu.eta, diflink = function(t) 1/(stats$mu.eta(stats$linkfun(t))),
                          T = function(etahat) diag(c(stats$mu.eta(etahat)))))
  fit <- c()
  fit1 <- br.fit(X, Y, link1)
  fit$call <- call
  fit$funlink <- link1
  fit$linkinv <- link1$linkinv
  fit$coefficients <- c(fit1$coeff, structure(fit1$phiest,
                                              .Names = c("phi")))
  fit$stder <- c(fit1$stder, fit1$phist)
  fit$zstats <- fit1$zstats
  fit$pvalues <- fit1$pvalues
  fit$fitted.values <- fit1$fitted
  fit$residuals <- fit1$res
  fit$k <- fit1$k
  fit$h <- fit1$h
  fit$GL <- fit1$GL
  fit$terms <- Terms
  fit$x <- X
  fit$y <- Y
  fit$resd <- fit1$resd
  fit$resstd <- fit1$resstd
  fit$Pseudo.R2 <- cbind(fit1$pseudor2)
     fit$etahat <- fit1$etahat
     fit$sigma2 <- fit1$sigma2
     fit$phi <- fit1$phi
  colnames(fit$Pseudo.R2) <- c("")
  rownames(fit$Pseudo.R2) <- c("")
  attr(fit, "na.message") <- attr(m, "na.message")
  class(fit) <- c("betareg", "lm")
  fit
}
"br.fit" <- function (x, y, link)
{
    x <- as.matrix(x)
    y <- as.matrix(y)
    linkfun <- link$linkfun
    linkinv <- link$linkinv
    mu.eta <- link$mu.eta
    diflink <- link$diflink
    T <- link$T
    ynew = linkfun(y)
    ystar = log(y/(1-y))
    ajuste = lm.fit(x, ynew)
    beta = c(ajuste$coef)
    k = length(beta)
    n = length(y)
    mean = fitted(ajuste)
    mean = linkinv(mean)
    dlink = diflink(mean)
    er = residuals(ajuste)
    sigma2 = sum(er^2)/((n - k) * (dlink)^2)
    phi = 1/n * sum(mean * (1 - mean)/sigma2 - 1)
    reg = c(beta, phi)
    loglik <- function(z) {
        z1 = z[1:k]
        z2 = z[k + 1]
        mu = linkinv(x %*% z1)
        sum(lgamma(z2) - lgamma(mu * z2) - lgamma((1 - mu) *
            z2) + (mu * z2 - 1) * log(y) + ((1 - mu) * z2 - 1) *
            log(1 - y))
    }
    loglikt <- function(z) {
        d = length(z) - 1
        z1 = z[1:d]
        z2 = z[d + 1]
        mu = z1
        lgamma(z2) - lgamma(mu * z2) - lgamma((1 - mu) * z2) +
            (mu * z2 - 1) * log(y) + ((1 - mu) * z2 - 1) * log(1 -
            y)
    }
    escore <- function(z) {
        z1 = z[1:k]
        z2 = z[k + 1]
        mu = linkinv(x %*% z1)
        munew = digamma(mu * z2) - digamma((1 - mu) * z2)
        T = diag(c(mu.eta(x %*% z1) ) )
        c(z2 * t(x) %*% T %*% (ystar - munew), sum(digamma(z2) -
            mu * digamma(mu * z2) - (1 - mu) * digamma((1 - mu) *
            z2) + mu * log(y) + (1 - mu) * log(1 - y)))
    }
    opt <- optim(reg, loglik, escore, method = "BFGS", control = list(fnscale = -1,
        maxit = 2000))
    if (opt$conv != 0)
        warning("FUNCTION DID NOT CONVERGE!")
    z <- c()
    coef <- (opt$par)[1:k]
    z$coeff <- coef
    z$beta <- beta
    z$phi <- phi
    etahat <- x %*% coef
    phihat <- opt$par[k + 1]
    muhat = linkinv(etahat)
    z$fitted <- muhat
    z$phiest <- phihat
    psi1 = trigamma(muhat * phihat)
    psi2 = trigamma((1 - muhat) * phihat)
    T1 = T(etahat)
    W = diag(c(phihat * (psi1 + psi2))) %*% T1^2
    vc = phihat * (psi1 * muhat - psi2 * (1 - muhat))
    D = diag(c(psi1 * (muhat^2) + psi2 * (1 - muhat)^2 - trigamma(phihat)))
    tempinv = solve(t(x) %*% W %*% x)
    g = sum(diag(D)) - (1/phihat) * t(vc) %*% t(T1) %*% x %*%
        tempinv %*% t(x) %*% T1 %*% vc
    K1 = tempinv %*% (c(g) * diag(k) + (1/phihat) * t(x) %*%
        T1 %*% vc %*% t(vc) %*% t(T1) %*% x %*% tempinv)
    K2 = -tempinv %*% t(x) %*% T1 %*% vc
    tempmatrix <- (-t(vc) %*% t(T1) %*% x %*% tempinv)
    tempmatrix <- cbind(tempmatrix, phihat)
    fisherinv = (1/(phihat * c(g))) * rbind(cbind(K1, K2), tempmatrix)
    stderrors <- sqrt(diag(fisherinv))[1:k]
    z$stderrors <- stderrors
    phier <- sqrt(diag(fisherinv))[k + 1]
    muhat <- as.vector(muhat)
    H = sqrt(W) %*% x %*% tempinv %*% t(x) %*% sqrt(W)
    h = diag(H)
    z$k = k
    z$h = h
    mustar = digamma(muhat * phihat) - digamma((1 - muhat) *
        phihat)
    Q = (phihat * (trigamma(muhat * phihat) + trigamma((1 - muhat) *
        phihat)) - (ystar - mustar) * etahat/(mu.eta(etahat))) *(mu.eta(etahat))^2
    Q <- as.vector(Q)
    Q <- diag(Q)
    f = vc - (ystar - mustar)
    e = -(y - muhat)/(y * (1 - y))
    XQXinv = solve(t(x) %*% Q %*% x)
    M = 1/(y * (1 - y))
    M = as.vector(M)
    M = diag(M)
    g = sum(diag(D)) - (1/phihat) * t(f) %*% t(T1) %*% x %*% tempinv %*% t(x) %*% T1 %*% f
    GL1 = T1 %*% x %*% XQXinv %*% t(x) %*% T1 %*% M
    GL2 = (1/(c(g) * phihat)) * T1 %*% x %*% XQXinv %*% t(x) %*%
        T1 %*% f %*% (t(f) %*% T1 %*% x %*% XQXinv %*% t(x) %*%
        T1 %*% M - t(e))
    GL = GL1 + GL2
    z$GL = GL
    z$fitted <- muhat
val <- 2 * (loglikt(c(y, phihat)) - loglikt(c(muhat, phihat)))
val[val<0] <- 0
    resd <- sign(y - muhat) * sqrt(val)
resstd <- sqrt(1.0+phihat)*(y-muhat)/sqrt((1.0-h)*muhat*(1-muhat));
z$resstd <- resstd
    z$resd <- resd
    z$phistd <- phier
    z$zstats <- coef/stderrors
    res <- y - muhat
    res <- as.vector(res)
    z$res <- res
    z$pvalues <- 2 * (1 - pnorm(abs(coef/stderrors)))
    pseudor2 <- cor(etahat, ynew)^2
    z$pseudor2 <- pseudor2
       z$etahat <- etahat
       sigma2 = sum(res^2)/((n - k) * (diflink(muhat)^2))
       z$sigma2 <- sigma2
    z
}
envelope.beta <- function(model=fit.model,sim=100,conf=.90, pch="+",font.main=1, cex.main=1.5, type = c("standardized","deviance")) {
   type <- match.arg(type)
  if(!any(class(model)=="betareg"))
  {stop("The model must be from the class betareg")}
  main = switch(type, standardized = "Half-Normal Plot of Standardized Residuals", deviance = "Half-Normal Plot of Deviance Residuals")
  ylab= switch(type, standardized = "Absolute Values of Standardized Residuals", deviance = "Absolute Values of Deviance Residuals")
  xlab="Half-Normal Quantile"
  alfa <-(1-conf)/2
  X <- model$x
  y <-model$y
  n <- nrow(X)
  p <- ncol(X)
  H <- X%*%solve(t(X)%*%X)%*%t(X)
  h <- diag(H)
  m <- model$fitted
  
  
###The line below is to avoid division by 0 when studentize the residuals, but trying to keep the leverage value high.
  h[round(h,15)==1]<-0.999999999999999
  
  si <- model$sigma2
  r <- model$residuals
  res <- switch(type, standardized = model$resstd, deviance = model$resd)
  e <- matrix(0,n,sim)
  e1 <- numeric(n)
  e2 <- numeric(n)
  phi <- model$phi
  
  for(i in 1:sim) {
    resp <- rbeta(n, m*phi, (1-m)*phi)
    fit <- betareg(resp~X-1)
    ti <- fit$residuals/(model$sigma2*sqrt(1-h))
    eo <- switch(type,standardized = sort(abs(fit$resstd)), deviance = sort(abs(fit$resd)))
    e[,i] <- eo
  }
  
  for(i in 1:n) {
    eo <- sort(e[i,])
    e1[i] <- quantile(eo,alfa)
    e2[i] <- quantile(eo,1-alfa)
  }
  
  med <- apply(e,1,median)
  qq <- qnorm((n+1:n+.5)/(2*n+1.125))
  plot(qq, sort(abs(res)), ylim=range(abs(res), e1, e2), pch=pch,
       main=main, xlab=xlab, ylab=ylab, cex.main=cex.main, font.main=font.main)
  
  lines(qq,e1,lty=1)
  lines(qq,e2,lty=1)
  lines(qq,med,lty=2) 
}

																												
"plot.betareg" <-
function (x, which = 1:4, caption = c("Deviance residuals vs indices of obs.", 
    "Standardized residuals vs indices of obs.", "Generalized leverage vs. Predicted values", "Cook's distance plot"), 
    panel = points, sub.caption = deparse(x$call), main = "", 
    ask = prod(par("mfcol")) < length(which) && dev.interactive(), 
    ..., id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75) 
{
    if (!inherits(x, "betareg")) 
        stop("Use only with 'betareg' objects")
    if (!is.numeric(which) || any(which < 1) || any(which > 4)) 
        stop("`which' must be in 1:4")
    rdev <- residuals(x,type="deviance")
    n <- length(rdev)
    h <- x$h
    k <- x$k
    gl <- x$GL
    gl <- diag(gl)
    yh <- predict(x)
    sr <- residuals(x)
    show <- rep(FALSE, 4)
    show[which] <- TRUE
    one.fig <- prod(par("mfcol")) == 1
    if (ask) {
        op <- par(ask = TRUE)
        on.exit(par(op))
    }
    if (show[1]) {
        plot(rdev, xlab = "Indices of obs.", ylab = "Deviance residuals", main = main, ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(caption[1], 3, 0.25)
        abline(h = 0, lty = 3, col = "gray")
    }
    if (show[2]) {
        plot(sr, xlab = "Obs. number", ylab = "Standardized residuals", main = main, ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(caption[2], 3, 0.25)
        abline(h = 0, lty = 3, col = "gray")
    }
    if (show[3]) {
	        plot(yh, gl, xlab = "Fitted values", ylab = "Generalized leverage", main = main,...)
        panel(yh, gl, ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(caption[3], 3, 0.25)
    }
    if (show[4]) {
	plot(h*(sr^2)/(k*(1-h)^2),xlab="Obs. number",ylab = "Cook's distance",type="h",main=main)
	mtext(caption[4], 3, 0.25)
       if (one.fig) 
            title(sub = sub.caption, ...)
    }
    if (!one.fig && par("oma")[3] >= 1) 
        mtext(sub.caption, outer = TRUE, cex = 1.25)
    invisible()
}
"predict.betareg" <-
function(object, terms = object$x[,-1],... ) 
{
   coef <- (object$coeff)[1:object$k]
   x <- terms 
   x <- matrix(x,length(x)/(ncol(object$x)-1),ncol(object$x)-1)
   x <- cbind(1,x)
   if(ncol(x) != ncol(object$x))
       stop("The number of columns must be equal to the number of coefficients.")
   pred = x%*%coef
   pred = object$linkinv(pred)
   pred = as.vector(pred)
   pred
}

"residuals.betareg" <-
function(object, type = c("standardized","usual","deviance"),...)
{
    type <- match.arg(type)
    muhat <- object$fitted
    r = object$y - muhat
    rd <- object$resd
    rd = as.vector(rd)
    phihat <- c(object$coef[object$k+1])
    sqrtvar <- sqrt(muhat*(1-muhat)/(1+phihat))
    r = as.vector(r)
    res <- switch(type, standardized = as.vector(r/sqrtvar),usual=r,deviance = rd)
        if (is.null(object$na.action)) 
        res
    else naresid(object$na.action, as.vector(r/sqrtvar))
}
"summary.betareg" <-
function (object,...) 
{
    z <- object
    ans <<- z[c("call", "terms")]
    class(ans) <<- "summary.betareg"
    ans$coefficients <<- z$coef
    ans$std <<- z$stder
    ans$zstats <<- z$zstats
    ans$pvalues <<- z$pvalues
    k <- z$k
    ans$est <<- cbind(format(ans$coeff),format(ans$std),c(format(ans$zstats[1:k]),""),c(format(ans$pvalues[1:k]),""))
    dimnames(ans$est) <- list(NULL, c("estimates", "std. errors","z-stats", "p-value"))
    rownames(ans$est) <- names(ans$coef)
    ans$Pseudo.R2 <- z$Pseudo.R2
    cat("\nCall:\n", deparse(ans$call), "\n\n", sep = "")
    if (length(coef(ans))) {
        cat("Coefficients:\n")
        print.default(format(ans$est),print.gap = 2, quote = FALSE)
        cat("\n")
        cat("Pseudo R^2:")
        print.default(format(ans$Pseudo.R2), print.gap = 2, quote = FALSE)
    }
    invisible(ans)
}
