ldframe <-
  function(...)
{
  ## $Id: ldframe.S,v 1.1 1997/06/24 09:24:11 smithdm1 Exp $
  obj <- data.frame(...)
  as.ldframe(obj)
}

as.ldframe <- function(dfobj, ...)
{
  UseMethod('as.ldframe')
}

as.ldframe.default <- function(dfobj, ...)
{
  obj <- as.data.frame(obj)
  as.ldframe.data.frame(obj, ...)
}

as.ldframe.data.frame <-
  function(dfobj, Subject, Time, sort.time = T)
{
  ## $Id: as.ldframe.data.frame.S,v 1.3 1997/11/07 16:41:58 smithdm1 Exp $
  ## make the data frame dfobj into an ldframe object
  ## the Subject is taken from dfobj, or can be specified with the Subject
  ##     argument (and evaluated in the frame of dfobj).  If not given, 
  ##     all rows are given Subject 1
  ## the Time column is taken from dfobj, or can be specified with the Time
  ##     arguement.  If not given, unit-spaced timepoints are given within
  ##     each subject
  ## The object is sorted by Subject, and by Time within Subject
  ## Make sure it has a Subject component
#  ## PJ
#  attach(dfobj)
  if(is.null(dfobj$Subject)) {
    if(missing(Subject)) {
      dfobj$Subject <- as.factor(rep(1, nrow(dfobj)))
    }
    else {
      new.s <- substitute(Subject)
      ## PJ
      # dfobj$Subject <- eval(new.s, local = dfobj)
      dfobj$Subject <- eval(new.s, envir = dfobj)
      if(is.name(new.s))
        dfobj[[as.character(new.s)]] <- NULL
    }
  }
  ## Don't allow NA's in Subject
  s <- dfobj$Subject
  if(any(is.na(s))) stop("Missing values not allowed in Subject column.")	
  ## Ensure Subjects are clustered
  ## PJ
  #rle.s <- rle(s)
  rle.s <- rle(as.vector(s))
  uni.s <- unique(s)
  if(length(rle.s$values) > length(uni.s)) {
    k <- length(rle.s$values)
    code <- match(rle.s$values, uni.s)
    ol <- sort.list(rep(code, rle.s$lengths))
    dfobj <- dfobj[ol,  ]
  }
  ## Make sure it hase a Time component
  if(is.null(dfobj$Time)) {
    if(missing(Time)) {
      s <- dfobj$Subject
      ## PJ
      # ns <- rle(s)$lengths
      ns <- rle(as.vector(s))$lengths
      dfobj$Time <- unlist(lapply(ns, seq))
    }
    else {
      new.t <- substitute(Time)
      ## PJ
      ## dfobj$Time <- eval(new.t, local = dfobj)
      dfobj$Time <- eval(new.t, envir = dfobj)
      if(is.name(new.t))
        dfobj[[as.character(new.t)]] <- NULL
    }
  }
  if(sort.time) {
    ## sort by time within subject
    ## PJ
    ## rle.s <- rle(dfobj$Subject)
    rle.s <- rle(as.vector(dfobj$Subject))
    m <- length(rle.s$values)
    ol <- order(rep(1:m, rle.s$lengths), dfobj$Time)
    dfobj <- dfobj[ol,  ]
  }
  ## Make sure Subject is a factor
  if(!is.factor(dfobj$Subject)) dfobj$Subject <- as.factor(dfobj$Subject)
  ## give it ldframe class
  if(!inherits(dfobj, "ldaframe"))
    class(dfobj) <- c("ldframe", "data.frame")
  dfobj
}


is.ldframe <- function(obj)
{
  ## $Id: is.ldframe.S,v 1.3 1997/09/18 16:20:32 smithdm1 Exp $
  ## only T if obj is of the right class, has a block-sorted Subject coln
  ## PJ
  ## inherits(obj, "ldframe") && !is.null(obj[["Subject"]]) && (length(rle(obj$Subject)$values) == length(unique(obj$Subject)))
  inherits(obj, "ldframe") && !is.null(obj[["Subject"]]) && (length(rle(as.vector(obj$Subject))$values) == length(unique(obj$Subject)))
}


plot.ldframe <- function(obj, y.expr, x.expr, line, color, pchar, selector, n.full, legend, ...)
{
### $Id: plot.ldframe.S,v 1.11 1998/01/20 12:02:36 smithdm1 Exp $
  plotline <- function(x, y, type = "l")
    {
      ## x and y are lists -- x should be sorted within lists
      ## pad out x and y with NA's etc and plot using current graphics options
      addna <- function(x)
        c(x, NA)
      y.out <- unlist(lapply(y, addna))
      t.out <- unlist(lapply(x, addna))
      lines(t.out, y.out, type = type)
    }
  missingline <- function(x, y, type = "l")
    {
      ## add lines where there are internal NAs
      ## returns list with elements "dropin" and "dropout"
      K <- length(x)
      x1 <- x2 <- y1 <- y2 <- NULL
      dropin.x <- dropin.y <- NULL
      dropout.x <- dropout.y <- NULL
      for(i in 1:K) {
        yi <- y[[i]]
        xi <- x[[i]]
        n <- length(yi)
        miss <- is.na(yi)
        if(sum(miss) > 0) {
          preidx <- (1:n)[!miss & c(miss[-1], F)]	# before an NA
          postidx <- (1:n)[!miss & c(F, miss[ - n])]	# after an NA
          npre <- length(preidx)
          if(length(postidx) && ((npre == 0) || postidx[1] <= preidx[1])) {
            ## drop-in
            dropin.x <- c(dropin.x, xi[postidx[1]])
            dropin.y <- c(dropin.y, yi[postidx[1]])
            postidx <- postidx[-1]
          }
          npost <- length(postidx)	## (npost==0) || (length(preidx)>length(postidx))
          if(npre > npost) {
            ## drop-out
            dropout.x <- c(dropout.x, xi[preidx[npre]])
            dropout.y <- c(dropout.y, yi[preidx[npre]])
            preidx <- preidx[ - npre]
          }
          x1 <- c(x1, xi[preidx])
          x2 <- c(x2, xi[postidx])
          y1 <- c(y1, yi[preidx])
          y2 <- c(y2, yi[postidx])
        }
      }
      if(((type == "l") || (type == "b") || (type == "l")) && length(x1)) {
        tso.miss <- tsopts("miss.style")[[1]]
        if(!is.null(tso.miss)) {
          ch.miss <- par(tso.miss)
          segments(x1, y1, x2, y2)
          if(!is.null(ch.miss))
            par(ch.miss)
        }
      }
      ## do drop-in and drop-out
      drops <- list(dropin = list(x = dropin.x, y = dropin.y), dropout = list(x = dropout.x, y = dropout.y))
      tso.in <- tsopts("dropin.style")[[1]]
      if(!is.null(tso.in)) {
        ch.in <- par(tso.in)
        points(drops$dropin)
        if(!is.null(ch.in))
          par(ch.in)
      }
      tso.out <- tsopts("dropout.style")[[1]]
      if(!is.null(tso.out)) {
        ch.out <- par(tso.out)
        points(drops$dropout)
        if(!is.null(ch.out))
          par(ch.out)
      }
    }
  parstrip <- function(input, defaults)
    {
      ## delete any occurrences like in defaults from input
      ## return a list of 2 elements:
      ## [[1]] is input with all defaults removed
      ## [[2]] is defaults updated from inputs.  named NULLs are deleted
      for(i in names(defaults)) {
        if(!is.null(names(input[i]))) {
          defaults[[i]] <- input[[i]]
          input[[i]] <- NULL
        }
        if(is.null(defaults[[i]]))
          defaults[[i]] <- NULL
      }
      list(input, defaults)
    }
  ## ensure ldframe status, and sort by Time within Subject
  obj <- as.ldframe(obj)
  Subj <- obj$Subject
  Time <- obj$Time
  tso.change <- tsopts(...)
  on.exit(tsopts(tso.change))
  if(missing(y.expr)) {
    yname <- names(obj)[1]
    y <- obj[[yname]]
  }
  else {
    y.expr <- substitute(y.expr)
    yname <- deparse(y.expr)
    ## PJ
    ##    y <- eval(y.expr, local = obj)
    y <- eval(y.expr, envir = obj)
  }
  if(missing(x.expr)) {
    xname <- "Time"
    x <- obj[[xname]]
  }
  else {
    x.expr <- substitute(x.expr)
    xname <- deparse(x.expr)
    ## PJ
    ## x <- eval(x.expr, local = obj)
    x <- eval(x.expr, local = obj)
  }
  if(is.null(y) || is.null(x) || (length(y) != length(x))) stop(
                                          "No data to plot, or lengths of x and y do not match.")	
  ## split sorts by subject! ack!
  ysplit <- rle.split(y, Subj)
  tsplit <- rle.split(x, Subj)
  m <- length(ysplit)
  yrange <- range(y, na.rm = T)
  trange <- range(x, na.rm = T)
  tso.general <- tsopts("general")[[1]]	## The following options, if included in the "general" tsopts
  ## option, will *not* be passed to par, but instead included in the
  ## top-level plot() call with the given default as shown (or not at
  ## all if NULL)
  defaults <- list(type = "l", xlab = xname, ylab = yname, xlim = NULL, ylim = NULL, axes = NULL, add = 
                   F, main = NULL, sub = NULL, log = NULL)
  general.strip <- parstrip(tso.general, defaults)
  defaults <- general.strip[[2]]
  if(exists("plot.ldframe.type", frame = 1)) {
    ## called from points.ldframe or lines.ldframe
    if(is.null(tso.general$type)) defaults$type <- get("plot.ldframe.type", frame = 1)
    defaults$add <- T
  }
  tso.general <- general.strip[[1]]
  ch.general <- par(tso.general)
  on.exit(if(!is.null(ch.general)) par(ch.general), T)
  type.general <- defaults$type
  defaults$type <- "n"
  defaults$log <- NULL	# too hard to handle!
  if(!(defaults$add)) do.call("plot", c(list(trange, yrange), defaults))	## eliminate shadows
  if(missing(n.full))
    n.full <- m
  else if(n.full < 1)
    n.full <- round(n.full * m)
  if(n.full < m) {
    if(n.full == 0) {
      shadow.x <- tsplit
      shadow.y <- ysplit
    }
    else {
      scores <- if(missing(selector)) runif(m) else sapply(ysplit, selector)
      perm <- sample(1:m)
      scores <- scores[perm]	## non.nas <- sort.list(scores, na.last=NA)
      ## avoid the warning in the above:
      non.nas <- sort.list(scores)
      non.nas <- non.nas[!is.na(scores[non.nas])]
      if(length(non.nas) < n.full) {
        n.full <- length(non.nas)
        if(n.full < 1)
          stop("selector function returned all missing values")
        warning(paste("selector function returned too many NAs: plotting only", n.full,
                      "subjects in full"))
      }
      selected <- non.nas[round(quantile(1:length(non.nas), seq(1, 0, length = n.full)))]	
      ## return subjects to original order
      selected <- perm[selected]
      shadow.x <- tsplit[ - selected]
      shadow.y <- ysplit[ - selected]
    }
    if(length(shadow.x) > 0) {
      type.shadow <- type.general
      tso.shadow <- tsopts("shadows")[[1]]
      if(!is.null(tso.shadow)) {
        ## tsopts(shadows=NULL) prevents shadow plotting
        if(!is.null(tso.shadow$type)) {
          type.shadow <- tso.shadow$type
          tso.shadow$type <- NULL
        }
        ch.shad <- par(tso.shadow)
        plotline(shadow.x, shadow.y, type.shadow)
        missingline(shadow.x, shadow.y, type.shadow)
        if(!is.null(ch.shad))
          par(ch.shad)
      }
    }
  }
  else {
    ## all subhects selected for plotting
    selected <- 1:m
  }
  if(n.full > 0) {
    ## plot selected lines
    obj.ssv <- ssv(obj)
    ## PJ
    ##    line.f <- as.factor(if(missing(line)) rep(1, m) else eval(substitute(line), local = obj.ssv))
    ##   color.f <- as.factor(if(missing(color)) rep(1, m) else eval(substitute(color), local = obj.ssv))
    ##     points.f <- as.factor(if(missing(pchar)) rep(1, m) else eval(substitute(pchar), local = obj.ssv))
    line.f <- as.factor(if(missing(line)) rep(1, m) else eval(substitute(line), envir = obj.ssv))
    color.f <- as.factor(if(missing(color)) rep(1, m) else eval(substitute(color), envir = obj.ssv))
    points.f <- as.factor(if(missing(pchar)) rep(1, m) else eval(substitute(pchar), envir = obj.ssv))
    #line.fc <- codes(line.f)[selected]
    #color.fc <- codes(color.f)[selected]
    #points.fc <- codes(points.f)[selected]
    #valeska
    line.fc <- unclass(line.f)[selected]
    color.fc <- unclass(color.f)[selected]
    points.fc <- unclass(points.f)[selected]


    sel.idx <- selected
    inter.f <- interaction(line.fc, color.fc, points.fc, drop = T)
    chunks <- levels(inter.f)
    ncolor <- length(levels(color.f))
    colors <- rep(tsopts("cols")[[1]], length = ncolor)
    nline <- length(levels(line.f))
    lines <- rep(tsopts("ltys")[[1]], length = nline)
    npchar <- length(levels(points.f))
    pchars <- rep(tsopts("pchs")[[1]], length = npchar)
    tso.line <- tsopts("full")[[1]]
    ch.line <- par(tso.line)
    ch.lcp <- par(c("col", "pch", "lty"))
    for(chunk in chunks) {
      chunkidx <- inter.f == chunk
      #coli <- colors[codes(color.f)[(sel.idx[chunkidx])[1]]]
      #ltyi <- lines[codes(line.f)[(sel.idx[chunkidx])[1]]]
      #pchi <- pchars[codes(points.f)[(sel.idx[chunkidx])[1]]]
      coli <- colors[unclass(color.f)[(sel.idx[chunkidx])[1]]]
      ltyi <- lines[unclass(line.f)[(sel.idx[chunkidx])[1]]]
      pchi <- pchars[unclass(points.f)[(sel.idx[chunkidx])[1]]]

      xi <- tsplit[sel.idx[chunkidx]]
      yi <- ysplit[sel.idx[chunkidx]]
      if(!missing(color))
        par(col = coli)
      if(!missing(line))
        par(lty = ltyi)
      if(!missing(pchar))
        par(pch = pchi)
      plotline(xi, yi, type.general)
      missingline(xi, yi, type.general)
    }
    par(ch.lcp)
    if(!is.null(ch.line))
      par(ch.line)
  }
  ## Now, draw a legend if required
  if(!missing(legend)) {
    lty <- NULL
    pch <- NULL
    col <- NULL
    type <- NULL
    labels <- NULL
    color.l <- levels(as.factor(color.f))
    color.len <- length(color.l)
    if(color.len > 1) {
      lty <- c(lty, rep(1, color.len))
      pch <- c(pch, rep("", color.len))
      col <- c(col, colors[1:color.len])
      type <- c(type, rep("l", color.len))
      labels <- c(labels, color.l)
    }
    pch.l <- levels(as.factor(points.f))
    pch.len <- length(pch.l)
    if(pch.len > 1) {
      lty <- c(lty, rep(1, pch.len))
      pch <- c(pch, pchars[1:pch.len])
      col <- c(col, rep(1, pch.len))
      type <- c(type, rep("p", pch.len))
      labels <- c(labels, pch.l)
    }
    lty.l <- levels(as.factor(line.f))
    lty.len <- length(lty.l)
    if(lty.len > 1) {
      lty <- c(lty, lines[1:lty.len])
      pch <- c(pch, rep("", lty.len))
      col <- c(col, rep(1, lty.len))
      type <- c(type, rep("l", lty.len))
      labels <- c(labels, lty.l)
    }
    if(length(lty)) {
      do.call("key", c(list(line = list(lty = lty, pch = pch, col = col), type = type), 
                       legend, 	## allow text elt in legend arg to supersede default
                       if(is.null(legend$text)) list(text = list(labels))))
    }
  }
  invisible()
}


tsopts <- function(...)
{
### $Id: tsopts.S,v 1.5 1997/09/18 16:23:09 smithdm1 Exp $
  nonmembers <- function(new, matchset)
    {
      n <- length(new)
      mr <- match(new, matchset, nomatch = NA)
      idx <- (1:n)[is.na(mr)]
      paste(new[idx], sep = "", collapse = ", ")
    }
  if(nargs() == 0) return(.ldats.options)
  current <- .ldats.options
  valid <- c("axes", "tslinetype", "tsline", "miss", "drop", "drop.style", "miss.style", "pchars", 
             "ltypes", "colors", "tschosen.type", "tschosen", "tsnotchosen.type", "tsnotchosen", 	
             ## new for plot.ldframe
             "shadows", "dropin.style", "dropout.style", "general", "full", "cols", "pchs", "ltys")
  temp <- list(...)
  if((length(temp) == 1) && is.null(names(temp)) && is.character(temp[[1]])) {
    arg <- temp[[1]]
    bad <- nonmembers(arg, valid)
    if(nchar(bad) > 0)
      stop(paste("Invalid option name(s): ", bad))
    else return(.ldats.options[arg])
  }
  ## a single list may be used as an arg
  if((length(temp) == 1) && is.list(temp[[1]]) && is.null(names(temp)))
    temp <- temp[[1]]
  n <- names(temp)
  if(is.null(n) || any(nchar(n) == 0))
    stop("options must be given by name")
  bad <- nonmembers(n, valid)
  if(nchar(bad) > 0)
    stop(paste("Invalid option(s): ", bad))
  changed <- current[n]
  current[n] <- temp
  .ldats.options <<- current
  invisible(changed)
}

.ldats.options<-
list(axes = list(), tslinetype = "l", tsline = list(lwd = 2, col = 5), miss = T,
	drop = T, drop.style = list(pch = 4, cex = 3), miss.style = list(lty = 
	2), tschosen = list(lwd = 2, cex = 2), tschosen.type = "o", tsnotchosen
	 = list(lwd = 1, lty = 1, col = 7), tsnotchosen.type = "l", pchars = 
	".", ltypes = 1, colors = c(2, 3, 4, 5, 6), shadows = list(col = 16), 
	dropin.style = NULL, dropout.style = list(pch = "x", col = 8), general
	 = list(), full = list(), cols = c(2, 3, 4, 5, 6, 7), pchs = c("a", "b",
	"c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", 
	"q", "r", "s", "t", "u", "v", "w", "x", "y", "z"), ltys = c(1, 2, 4, 3, 
	5, 6, 7, 8))

rle.split <- function(data, group)
{
## $Id: rle.split.S,v 1.1 1997/11/14 12:04:39 smithdm1 Exp $
## An internal function to Oswald, to correct deficiencies in split.default.
## This is a safer version of split, based on rle:
## - the order of the list will be the same as rle(group)
## - the result list will have as many elements as rle(group)
## - unrepresented levels in group are ignored
	rl <- rle(as.vector(group))
	m <- length(rl$lengths)
	splitcode <- rep(1:m, rl$lengths)
	result <- split(data, splitcode)
	names(result) <- NULL
	result
}

ssv <- function(obj, name, warn)
{
  ## $Id: ssv.S,v 1.6 1998/01/22 14:38:06 smithdm1 Exp $
  ## return time-independent variables
  ## if name is not given, return all time-independent variables
  ## if name is given, return only those specified
  ## if warn is T, give warnings for nonexistent or non-time-ind variables
  ## warn is T by default unless name is missing
  if(!is.ldframe(obj)) stop("First argument must be ldframe object")
  if(missing(warn))
    warn <- !missing(name)
  if(missing(name))
    name <- names(obj)
  result <- list()
  resnames <- NULL
  s <- obj[["Subject"]]	## using as.numeric(s) below prevents problems if there are excess levels
  ## PJ
  ## rle.full.s <- rle(s)
  rle.full.s <- rle(as.vector(s))
  m <- length(rle.full.s$values)
  N <- length(s)
  for(n in name) {
    v <- obj[[n]]
    if(!is.null(v)) {
      ## first check that v *is* time-independent
      naidx <- is.na(v)
      s.nona <- s[!naidx]
      v.nona <- v[!naidx]
      rle.s <- if(any(naidx)) rle(as.vector(s.nona)) else rle.full.s
      nsubj <- length(rle.s$values)
      if(nsubj > 0) {
        if(nsubj == 1)
          firstidx <- 1
        else firstidx <- c(1, 1 + cumsum(rle.s$lengths[ - nsubj]))
        ideal <- rep(v.nona[firstidx], rle.s$lengths)
        if(all(ideal == v.nona)) {
          result.index <- rep(N + 1, m)	# N+1 index will give an NA
          nona.index <- ((1:N)[!naidx])[firstidx]
          non.na.subject.index <- match(rle.s$values, rle.full.s$values)
          result.index[non.na.subject.index] <- nona.index
          result <- c(result, list(v[result.index]))
          resnames <- c(resnames, n)
        }
        else if(warn)
          warning(paste("Variable", n, "exists, but is not subject-specific"))
      }
      else {
        ## nsubj==0 -- v is all NAs
        ## result <- c(result, list(rep(NA,m)))
        ## resnames <- c(resnames,n)
        if(warn) warning(paste("Variable", n, "exists, but is all NAs"))
      }
    }
    else if(warn)
      warning(paste("Variable", n, "does not exist"))
  }
  names(result) <- resnames
  result <- as.data.frame(result)
  if(nrow(result) == 0)
    result <- NULL
  result
}

identify.ldframe <- function(obj, y.expr, x.expr, labels, ...)
{
  ## $Id: identify.ldframe.S,v 1.2 1998/03/10 14:09:13 smithdm1 Exp $
  if(missing(labels)) labels <- paste(row.names(obj), "[", obj$Subject, "]", sep = "")
  if(missing(y.expr)) {
    yname <- names(obj)[1]
    y <- obj[[yname]]
  }
  else {
    y.expr <- substitute(y.expr)
    yname <- deparse(y.expr)
    y <- eval(y.expr, envir = obj)
  }
  if(missing(x.expr)) {
    xname <- "Time"
    x <- obj[[xname]]
  }
  else {
    x.expr <- substitute(x.expr)
    xname <- deparse(x.expr)
    x <- eval(x.expr, envir = obj)
  }
  identify(x, y, labels, ...)
}

#---------------------------------------
# Funçoes das bibliotecas Gee e Geepack
# Modificadas por Valeska
#---------------------------------------


summary.gee <- function (object, correlation = TRUE, ...) 
{
    coef <- object$coefficients
    resid <- object$residuals
    n <- length(resid)
    p <- object$rank
    if (is.null(p)) 
        p <- sum(!is.na(coef))
    if (!p) {
        warning("This model has zero rank --- no summary is provided")
        return(object)
    }
    nas <- is.na(coef)
    cnames <- names(coef[!nas])
    coef <- matrix(rep(coef[!nas], 5), ncol = 5)
    dimnames(coef) <- list(cnames, c("Estimate", "Naive S.E.", 
         "Robust S.E.", "Robust z", "P-valor"))
    rse <- sqrt(diag(object$robust.variance))
    nse <- sqrt(diag(object$naive.variance))
    coef[, 2] <- nse
    coef[, 3] <- rse
    coef[, 4] <- coef[, 1]/coef[, 3]
    coef[, 5] <- 2*(1-pnorm(abs(coef[, 4])))
    summary <- list()
    summary$call <- object$call
    summary$version <- object$version
    summary$nobs <- object$nobs
    summary$residual.summary <- quantile(as.vector(object$residuals))
    names(summary$residual.summary) <- c("Min", "1Q", "Median", 
        "3Q", "Max")
    summary$model <- object$model
    summary$title <- object$title
    summary$coefficients <- coef
    summary$working.correlation <- object$working.correlation
    summary$scale <- object$scale
    summary$error <- paste("Error code was", object$error)
    summary$working.correlation <- object$working.correlation
    summary$iterations <- object$iterations
    if (correlation) {
    }
    attr(summary, "class") <- "summary.gee"
    summary
}


summary.geese<-function(object,...){
	estimate<-object$beta
	naive.se<-sqrt(diag(object$vbeta.naiv))
	#z.naive<-estimate/naive.se
	robust.se<-sqrt(diag(object$vbeta))
	z.robust<-estimate/robust.se
	p.value<-2*(1-pnorm(abs(z.robust)))

	mean.sum<-data.frame(estimate, naive.se, robust.se, z.robust, p.value)
	rownames(mean.sum) <- object$xnames
	
	corr.sum <- data.frame(estimate = object$alpha, robust.se = sqrt(diag(object$valpha)))
    	corr.sum$wald.z <- corr.sum$estimate/corr.sum$robust.se
    	corr.sum$p <- 1 - pnorm(corr.sum$wald.z)
    	if (object$model$corstr != "independence"){
    		rownames(corr.sum) <- object$zcor.names
    	}
    	
    	scale.sum <- data.frame(estimate = object$gamma, robust.se = sqrt(diag(object$vgamma)))
    	scale.sum$wald.z <- scale.sum$estimate/scale.sum$robust.se
    	scale.sum$p <- 1 - pnorm(scale.sum$wald.z)
    	
    	ans <- list(mean = mean.sum, correlation = corr.sum, scale = scale.sum, 
        call = object$call, model = object$model, control = object$control, 
        error = object$err, clusz = object$clusz)
    	class(ans) <- "summary.geese"
    	ans
}

geese<-
function (formula = formula(data), sformula = ~1, id, waves = NULL, 
    data = parent.frame(), subset = NULL, na.action = na.omit, 
    contrasts = NULL, weights = NULL, zcor = NULL, corp = NULL, 
    control = geese.control(...), b = NA, alpha = NA, gm = NA, 
    family = gaussian(), mean.link = NULL, variance = NULL, cor.link = "identity", 
    sca.link = "identity", link.same = TRUE, scale.fix = FALSE, 
    scale.value = 1, corstr = "independence", ...) 
{
    scall <- match.call()
    mnames <- c("", "formula", "data", "offset", "weights", "subset", 
        "na.action", "id", "waves")
    cnames <- names(scall)
    cnames <- cnames[match(mnames, cnames, 0)]
    mcall <- scall[cnames]
    if (is.null(mcall$id)) 
        mcall$id <- as.name("id")
    mcall[[1]] <- as.name("model.frame")
    m <- eval(mcall, parent.frame())
    y <- model.extract(m, response)
    if (is.null(dim(y))) 
        N <- length(y)
    else N <- dim(y)[1]
    mterms <- attr(m, "terms")
    x <- model.matrix(mterms, m, contrasts)
    offset <- model.extract(m, offset)
    if (is.null(offset)) 
        offset <- rep(0, N)
    w <- model.extract(m, weights)
    if (is.null(w)) 
        w <- rep(1, N)
    id <- model.extract(m, id)
    waves <- model.extract(m, waves)
    if (is.null(id)) 
        stop("id variable not found.")
    mcall$formula[3] <- switch(match(length(sformula), c(0, 2, 
        3)), 1, sformula[2], sformula[3])
    m <- eval(mcall, parent.frame())
    terms <- attr(m, "terms")
    zsca <- model.matrix(terms, m, contrasts)
    soffset <- model.extract(m, offset)
    if (is.null(soffset)) 
        soffset <- rep(0, N)
    if (is.character(family)) 
        family <- get(family)
    if (is.function(family)) 
        family <- family()
    ans <- geese.fit(x, y, id, offset, soffset, w, waves, zsca, 
        zcor, corp, control, b, alpha, gm, family, mean.link, 
        variance, cor.link, sca.link, link.same, scale.fix, scale.value, 
        corstr, ...)
    #fit <- x %*% ans$beta
    ans <- c(ans, list(call = scall, formula = formula))
    class(ans) <- "geese"
    ans
}

resumo.grupo <- function (x, y, na.rm=F)
{
	min <- tapply(x, y, min, na.rm=na.rm)
	max <- tapply(x, y, max, na.rm=na.rm)
	mediana <- tapply(x, y, median, na.rm=na.rm)
	media <- tapply(x, y, mean, na.rm=na.rm)
	sd <- tapply(x, y, sd, na.rm=na.rm)
	cbind(Mínimo=min, Máximo=max, Mediana=mediana, Média=media, DP=sd)
}


or.glm <- function(x, ic=0.95)
{
	if (!all(class(x)==c("glm", "lm"))){stop("Essa função só pode ser usada em objetos da classe 'glm'")}
	x<-summary.glm(x)
	tabela <- cbind(Coeficientes=x$coefficients[,1], "Erro Padrão"=x$coefficients[,2], OR=exp(x$coefficients[,1]), "Inf"=rep(0, length(x$coefficients[,1])), Sup=rep(0, length(x$coefficients[,1])), "p-valor"=x$coefficients[,4])
	tabela[,4]<-exp(tabela[,1]-(tabela[,2]*qnorm(1-((1-ic)/2))))
	tabela[,5]<-exp(tabela[,1]+(tabela[,2]*qnorm(1-((1-ic)/2))))
	tabela[1,3:5]<-rep(NaN,3)
	round(tabela,4)
}


or.gee <- function(x, ic=0.95)
{
	if (!all(class(x)==c("gee", "glm"))){stop("Essa função só pode ser usada em objetos da classe 'gee'")}
	x<-summary.gee(x)
	tabela <- cbind(Coeficientes=x$coefficients[,1], "Erro Padrão"=x$coefficients[,3], OR=exp(x$coefficients[,1]), "Inf"=rep(0, length(x$coefficients[,1])), Sup=rep(0, length(x$coefficients[,1])), "p-valor"=x$coefficients[,5])
	tabela[,4]<-exp(tabela[,1]-(tabela[,2]*qnorm(1-((1-ic)/2))))
	tabela[,5]<-exp(tabela[,1]+(tabela[,2]*qnorm(1-((1-ic)/2))))
	tabela[1,3:5]<-rep(NaN,3)
	round(tabela,4)
}

#-------------------
perfil.std<-function(id,y,group=NULL, flag=TRUE){
	#argumentos: id = identificação do indivíduo, y= variável resposta
	# no formato matrix, e group = grupos
	# retorna grafico se flag=TRUE e retorna os dados padronizados
	#os argumentos vem de um dado np formato wide
	if (is.null(group)){
		media<-	matrix(apply(y,2,mean, na.rm=T),ncol=ncol(y),nrow=nrow(y),byrow=T)
		dp<-matrix(apply(y,2,sd, na.rm=T),ncol=ncol(y),nrow=nrow(y),byrow=T)
		ystd<-(y-media)/dp
		ystd<-data.frame(id,ystd)
		ylong<-reshape(ystd,idvar="id",varying=list(names(ystd)[2:ncol(ystd)]),
		       timevar="time",direction='long',v.names="y")
		ylong<-ylong[order(ylong$id),]
		ylong<-as.ldframe(ylong,id,time)
		if (flag){
			plot(ylong,y)
			title('Perfil Individual Padronizado')
		}
	}
	else{
		namegroup<- unique(group)
		ycomp<-data.frame(NULL)
		for (i in 1:length(namegroup)){
			ygroup<-y[group==namegroup[i],]
			idgroup<-id[group==namegroup[i]]
			media<-	matrix(apply(ygroup,2,mean, na.rm=T),ncol=ncol(ygroup),nrow=nrow(ygroup),byrow=T)
			dp<-matrix(apply(ygroup,2,sd, na.rm=T),ncol=ncol(ygroup),nrow=nrow(ygroup),byrow=T)
			ystd<-(ygroup-media)/dp
			ystd<-data.frame(idgroup,ystd)
			ycomp<-rbind(ycomp,ystd)
			ylong<-reshape(ystd,idvar="idgroup",varying=list(names(ystd)[2:ncol(ystd)]),
		       		timevar="time",direction='long',v.names="y")
			ylong<-ylong[order(ylong$id),]
			ylong<-as.ldframe(ylong,idgroup,time)
			if (flag){
				win.graph()
				plot(ylong,y,cols=i)
				title(paste('Grupo',namegroup[i],sep=''))
			}
		}
		ystd<-ycomp
	}
	ystd
}