home *** CD-ROM | disk | FTP | other *** search
/ Big Green CD 8 / BGCD_8_Dev.iso / NEXTSTEP / UNIX / Educational / R-0.49-MI / R-0.49-M / library / gee < prev    next >
Encoding:
Text File  |  1997-09-13  |  12.3 KB  |  397 lines
  1. library.dynam("gee.so")
  2.  
  3. "model.extract" <- function (frame, component) 
  4. {
  5.     component <- as.character(substitute(component))
  6.     rval <- switch(component, response = model.response(frame), offset = model.offset(frame
  7.     ), weights = frame$"(weights)", start = frame$"(start)",id=frame$"(id)")
  8.     if (length(rval) == nrow(frame)) 
  9.         names(rval) <- attr(frame, "row.names")
  10.     else if (is.matrix(rval) && nrow(rval) == nrow(frame)) {
  11.         t1 <- dimnames(rval)
  12.         dimnames(rval) <- list(attr(frame, "row.names"), t1[[2]])
  13.     }
  14.     return(rval)
  15. }
  16. "as.family" <-
  17. function (family) 
  18. {
  19.     if (is.character(family)) 
  20.         family <- get(family)
  21.     if (is.function(family)) 
  22.         family <- family()
  23.     if (is.null(family$family)) 
  24.         stop("family not recognised")
  25.         family
  26. }
  27. # gee Splus support @(#) geeformula.q 4.4 96/09/27
  28. print.gee <- function(x, digits = NULL, quote = F, prefix = "")
  29. {
  30. # gee Splus support  @(#) geeformula.q 4.4 96/09/27
  31.         if(is.null(digits)) digits <- options()$digits else options(digits =
  32.                         digits)
  33.         cat("\n", x$title)
  34.         cat("\n", x$version, "\n")
  35.         cat("\nModel:\n")
  36.         cat(" Link:                     ", x$model$link, "\n")
  37.         cat(" Variance to Mean Relation:", x$model$varfun, "\n")
  38.         if(!is.null(x$model$M))
  39.                 cat(" Correlation Structure:    ", x$model$corstr, ", M =", x$
  40.                         model$M, "\n")
  41.         else cat(" Correlation Structure:    ", x$model$corstr, "\n")
  42.         cat("\nCall:\n")
  43.         dput(x$call)    #       cat("\nTerms:\n")
  44. ###        ys <- matrix(rep(as.matrix(x$id, ncol = 1), 5), ncol = 5)
  45.         ys <- matrix(rep(matrix(x$id, ncol = 1), 5), ncol = 5)
  46.         ys[, 2] <- x$y
  47.         ys[, 3] <- x$linear.predictors
  48.         ys[, 4] <- x$fitted.values
  49.         ys[, 5] <- x$residuals
  50.         dimnames(ys) <- list(1:length(x$y), c("ID", "Y", "LP", "fitted",
  51.                 "Residual"))    #       cat("\nFitted Values:\n")
  52.         cat("\nNumber of observations : ", x$nobs, "\n")
  53.         cat("\nMaximum cluster size   : ", x$max.id, "\n")
  54.         nas <- x$nas
  55.         if(any(nas))
  56.                 cat("\n\nCoefficients: (", sum(nas),
  57.                         " not defined because of singularities)\n", sep = "")
  58.         else cat("\n\nCoefficients:\n")
  59.         print(x$coefficients, digits = digits)
  60.         cat("\nEstimated Scale Parameter: ", format(round(x$scale, digits)))
  61.         cat("\nNumber of Iterations: ", x$iterations)
  62. ###<TSL>
  63.     four<-min(4,dim(x$working.correlation)[[1]])
  64.         cat("\n\nWorking Correlation[1:",four,",1:",four,"]\n")
  65.         print(x$working.correlation[1:four, 1:four], digits = digits)
  66. ###</TSL>
  67.         cat("\n\nReturned Error Value:\n")
  68.         print(x$error)
  69.         invisible(x)
  70. }
  71.  
  72. print.summary.gee <- function(x, digits = NULL, quote = F, prefix = "" )
  73. {
  74. # gee Splus support @(#) geeformula.q 4.4 96/09/27
  75.     if(is.null(digits))
  76.         digits <- options()$digits
  77.     else options(digits = digits)
  78.     cat("\n",x$title)
  79.     cat("\n",x$version,"\n")
  80.     cat("\nModel:\n")
  81.     cat(" Link:                     ",x$model$link,"\n")
  82.     cat(" Variance to Mean Relation:",x$model$varfun,"\n")
  83.     if(!is.null(x$model$M))
  84.             cat(" Correlation Structure:    ",x$model$corstr,", M =",x$model$M,"\n")
  85.     else     cat(" Correlation Structure:    ",x$model$corstr,"\n")
  86.     cat("\nCall:\n")
  87.     dput(x$call)
  88.     cat("\nSummary of Residuals:\n")
  89.     print(x$residual.summary, digits = digits)
  90.     nas <- x$nas
  91. ###    if(any(nas))
  92.     if(!is.null(nas) && any(nas))
  93.         cat("\n\nCoefficients: (", sum(nas), 
  94.             " not defined because of singularities)\n", sep = "")
  95.     else cat("\n\nCoefficients:\n")
  96.     print(x$coefficients, digits = digits)
  97.     cat("\nEstimated Scale Parameter: ", format(round(x$scale,digits)))
  98.     cat("\nNumber of Iterations: ", x$iterations)
  99.     cat("\n\nWorking Correlation\n")
  100.     print(x$working.correlation,digits=digits)
  101.     if(!is.null(x$naive.correlation)) {
  102.         cat("\n\nNaive Correlation of Estimates:\n")
  103.         print(x$naive.correlation)
  104.     }
  105.     if(!is.null(x$robust.correlation)) {
  106.         cat("\n\nRobust Correlation of Estimates:\n")
  107.         print(x$robust.correlation)
  108.     }
  109.     invisible(x)
  110. }
  111. summary.gee <- function(object, correlation = T)
  112. {
  113. # gee Splus support @(#) geeformula.q 4.4 96/09/27
  114.     coef <- object$coefficients
  115.     resid <- object$residuals
  116.     n <- length(resid)
  117.     p <- object$rank
  118.     if(is.null(p))
  119.         p <- sum(!is.na(coef))
  120.     if(all(is.na(coef))) {
  121.         warning("This model has zero rank --- no summary is provided")
  122.         return(object)
  123.     }
  124.     nas <- is.na(coef)
  125.     cnames <- names(coef[!nas])
  126.     coef <- matrix(rep(coef[!nas], 5), ncol = 5)
  127.     dimnames(coef) <- list(cnames, c("Estimate", 
  128.             "Naive S.E.",  "Naive z", 
  129.                 "Robust S.E.", "Robust z"))
  130.     rse <- sqrt(diag(object$robust.variance))
  131.     nse <- sqrt(diag(object$naive.variance))
  132.     coef[,2] <- nse
  133.     coef[,3] <- coef[,1]/coef[,2]
  134.     coef[,4] <- rse
  135.     coef[,5] <- coef[,1]/coef[,4]
  136.     summary <- list()
  137.     summary$call <- object$call
  138.     summary$version <- object$version
  139.     summary$nobs <- object$nobs
  140.     summary$residual.summary <- quantile(as.vector(object$residuals))
  141.     names(summary$residual.summary) <- c("Min", "1Q", "Median", "3Q", "Max")
  142.     summary$model<- object$model
  143.     summary$title <- object$title
  144.     summary$coefficients <- coef
  145.     summary$working.correlation <- object$working.correlation
  146.     summary$scale <- object$scale
  147.     summary$error <- paste("Error code was", object$error)
  148.     summary$working.correlation <- object$working.correlation
  149.     summary$iterations <- object$iterations
  150.     if ( correlation ) {
  151.     #    rob.var <- object$robust.variance
  152.     #    nai.var <- object$naive.variance
  153.     #    summary$robust.correlation <- rob.var /
  154.     #    outer(sqrt(diag(rob.var)),sqrt(diag(rob.var)))
  155.     #    dimnames(summary$robust.correlation) <- list(object$xnames,object$xnames)
  156.     #    summary$naive.correlation <- nai.var /
  157.     #    outer(sqrt(diag(nai.var)),sqrt(diag(nai.var)))
  158.     #    dimnames(summary$naive.correlation) <- list(object$xnames,object$xnames)
  159.     }
  160.     attr(summary,"class") <- "summary.gee"
  161.     summary
  162. }
  163. gee <- function(formula = formula(data), id = id, data = sys.parent(), 
  164.     subset, na.action, offset = NA, R = NA, b = NA, tol = 
  165.     0.001, maxiter = as.integer(25), family = gaussian, corstr = 
  166.     "independence", Mv = 1, silent = T, contrasts = NULL, scale.fix = F, 
  167.     scale.value = 1)
  168. {
  169. # gee Splus support @(#) geeformula.q 4.4 96/09/27
  170.     print("Beginning Cgee S-function, @(#) geeformula.q 4.4 96/09/27"
  171.         )
  172.     call <- match.call()
  173.     m <- match.call(expand = F)
  174.     m$R <- m$b <- m$tol <- m$maxiter <- m$link <- m$varfun <-
  175.                m$corstr <- m$Mv <- m$silent <- m$contrasts <-
  176.                m$family <- m$scale.fix <- m$scale.value <- NULL
  177.     if(is.null(m$id)) {
  178.         m$id <- as.name("id")
  179.     }
  180.     if(is.null(m$na.action)) {
  181.     }
  182.     else {
  183.         if(m$na.action != "na.omit") {
  184.             print("Warning: Only na.omit is implemented for gee")
  185.             print("         continuing with na.action=na.omit")
  186.             m$na.action <- as.name("na.omit")
  187.         }
  188.     }
  189.     m[[1]] <- as.name("model.frame")
  190. ##<tsl>
  191.     m$use.data<-T
  192. ##</tsl>
  193.     m <- eval(m, sys.parent())
  194.     Terms <- attr(m, "terms")
  195.     y <- as.matrix(model.extract(m, response))
  196. ## <tsl>    x <- model.matrix(Terms, m, contrasts)
  197. ### Rchange
  198.     if (is.null(contrasts))
  199.         x <- model.matrix(Terms, m)
  200.     else stop("R model.matrix() can't handle contrasts argument")
  201. ### </tsl>
  202.     N <- rep(1, length(y))
  203.     if (dim(y)[2]==2) {
  204.          N <- as.vector(y %*% c(1,1))
  205.          y <- y[,1]
  206.     }
  207.     else {
  208.              if (dim(y)[2]>2)
  209.          stop("Only binomial response matrices (2 columns)")
  210.     }
  211.     offset <- model.extract(m, offset)
  212.     id <- model.extract(m, id)
  213.     if(is.null(id)) {
  214.         stop("Id variable not found")
  215.     }
  216.     nobs <- nrow(x)
  217.     nc <- ncol(x)
  218.     xnames <- dimnames(x)[[2]]
  219.     if(is.null(xnames)) {
  220.         xnames <- paste("x", 1:nc, sep = "")
  221.         dimnames(x) <- list(NULL, xnames)
  222.     }
  223.     family <- as.family(family)
  224.     if(deparse(substitute(b)) != "NA") {
  225.         b <- matrix(b, ncol = 1)
  226.         Rb <- nrow(b)
  227.         if(Rb != nc) {
  228.             stop("Dim beta != ncol(x)")
  229.         }
  230.         start.beta <- as.integer(1)
  231.         beta <- matrix(as.double(b), ncol = 1)
  232.     }
  233.     else {  print("running glm to get initial regression estimate")
  234. ## <tsl>    beta <- as.numeric(glm(m, family = family)$coef)
  235.         mm <- match.call(expand = F)
  236.         mm$R <- mm$b <- mm$tol <- mm$maxiter <- mm$link <- mm$varfun <-mm$corstr <- mm$Mv <- mm$silent <- mm$contrasts <-mm$scale.fix <- mm$scale.value <- mm$id<-NULL
  237.         mm[[1]]<-as.name("glm")
  238.         beta <- as.numeric(eval(mm,sys.frame(sys.parent()))$coef)
  239. ###    </tsl>
  240.         print(beta)
  241.     }
  242.     if(length(id) != length(y)) {
  243.         stop("Id and y not same length")
  244.     }
  245.     p <- ncol(x)
  246.     maxclsz <- as.integer(max(unlist(lapply(split(id, id), 
  247.         "length"))))
  248.     maxiter <- as.integer(maxiter)
  249.     silent <- as.integer(silent)
  250.     if(is.null(offset)) {
  251.         offset <- rep(0, length(y))
  252.     }
  253.     else if(length(offset) == 1 && offset == 0) {
  254.         offset <- rep(0, length(y))
  255.     }
  256.     else if(length(offset) != length(y)) {
  257.         print("Warning: offset and y are not of same length")
  258.         print("         offset being replicated to length of y")
  259.         offset <- rep(offset, length(y))
  260.     }
  261.     offset <- as.double(offset)
  262.     if(!is.na(R[1])) {
  263.         Rr <- nrow(R)
  264.         Rc <- ncol(R)
  265.         R <- matrix(R, ncol = Rc)
  266.         if(Rr != Rc) {
  267.             stop("R is not square!")
  268.         }
  269.         if(Rr < maxclsz) {
  270.             stop("R not big enough to accommodate some clusters.")
  271.         }
  272.     }
  273.     else {
  274.         R <- matrix(as.double(rep(0, maxclsz * maxclsz)), nrow = 
  275.             maxclsz)
  276.     }
  277.     links <- c("identity", "logarithm", "logit", "reciprocal", "probit",
  278. "cloglog")
  279.     varfuns <- c("gaussian", "poisson", "binomial", "gamma")
  280.     corstrs <- c("independence", "fixed", "stat_M_dep", "non_stat_M_dep", 
  281.         "exchangeable", "AR-M", "unstructured")
  282. #    defam <- function(x)
  283. #    {
  284. # function to get the essential components from family arg
  285. #        xfam <- unclass(family(x))$family
  286. #        c(xfam["link"], xfam["variance"])
  287. #    }
  288. #    glmlinks <- unlist(glm.links["name",  ])
  289. #    geelinks <- glmlinks[c(1, 5, 2, 6, 4, 3)]    # match ordering known to cgee
  290. #    glmvars <- unlist(glm.variances["name",  ])
  291. #    geevars <- glmvars[c(1, 3, 2, 4)]
  292.     thisfam <- family$family
  293.     link <- family$link
  294.     varfun <- family$family
  295.     linkv <- as.integer(match(c(link), links, -1))
  296.     varfunv <- as.integer(match(c(varfun), varfuns, -1))
  297.     corstrv <- as.integer(match(c(corstr), corstrs, -1))
  298.     if(linkv < 1)
  299.         stop("unknown link.")
  300.     if(varfunv < 1)
  301.         stop("unknown varfun.")
  302.     if(corstrv < 1)
  303.         stop("unknown corstr.")
  304.     naivvar <- matrix(rep(0, p * p), nrow = p)
  305.     robvar <- matrix(rep(0, p * p), nrow = p)
  306.     phi <- as.double(scale.value)
  307.     scale.fix <- as.integer(scale.fix)
  308.     errstate <- as.integer(1)
  309.     tol <- as.double(tol)
  310.     Mv <- as.integer(Mv)
  311.     maxcl <- as.integer(0)
  312.     if(!(is.double(x)))
  313.         x <- as.double(x)
  314.     if(!(is.double(y)))
  315.         y <- as.double(y)
  316.     if(!(is.double(id)))
  317.         id <- as.double(id)
  318.     if(!(is.double(N)))
  319.         N <- as.double(N)
  320.     modvec <- as.integer(c(linkv, varfunv, corstrv))
  321.     z <- .C("Cgee",
  322.         x,
  323.         y,
  324.         id,
  325.         N,
  326.         offset,
  327.         nobs,
  328.         p,
  329.         modvec,
  330.         Mv,
  331.         estb = beta,
  332.         nv = naivvar,
  333.         rv = robvar,
  334.         sc = phi,
  335.         wcor = R,
  336.         tol,
  337.         mc = maxcl,
  338.         iter = maxiter,
  339.         silent,
  340.         err = errstate,
  341.         scale.fix)
  342.     if(z$err != 0)
  343.         warning(paste("Note: Cgee had an error (code=", z$err, 
  344.             ").  Results suspect."))
  345.     if(min(eigen(z$wcor)$values) < 0) {
  346.         warning("Working correlation estimate not positive definite")
  347.         z$err <- z$err + 1000
  348.     }
  349.     fit <- list()
  350. ## <tsl>    attr(fit, "class") <- c("gee", "glm")
  351.     fit$title <- "GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA"
  352.     fit$version <- "gee S-function, version 4.4 modified 96/09/27 (1996)"
  353.     links <- c("Identity", "Logarithm", "Logit", "Reciprocal", "Probit",
  354. "Cloglog")
  355.     varfuns <- c("Gaussian", "Poisson", "Binomial", "Gamma")
  356.     corstrs <- c("Independent", "Fixed", "Stationary M-dependent", 
  357.         "Non-Stationary M-dependent", "Exchangeable", "AR-M", 
  358.         "Unstructured")
  359.     fit$model <- list()
  360.     fit$model$link <- links[linkv]
  361.     fit$model$varfun <- varfuns[varfunv]
  362.     fit$model$corstr <- corstrs[corstrv]
  363.     if(!is.na(match(c(corstrv), c(3, 4, 6))))
  364.         fit$model$M <- Mv
  365.     fit$call <- call
  366.     fit$terms <- Terms
  367.     fit$formula <- as.vector(attr(Terms, "formula"))
  368.     fit$contrasts <- attr(x, "contrasts")
  369.     fit$nobs <- nobs
  370.     fit$iterations <- z$iter
  371.     fit$coefficients <- as.vector(z$estb)
  372.     fit$nas <- is.na(fit$coefficients)
  373.     names(fit$coefficients) <- xnames
  374.     eta <- as.vector(x %*% fit$coefficients)
  375.     fit$linear.predictors <- eta
  376. ##Rchange
  377.     mu <- as.vector(family$linkinv(eta))
  378. ##
  379.     fit$fitted.values <- mu
  380.     fit$residuals <- y - mu
  381.     fit$family <- family
  382.     fit$y <- as.vector(y)
  383.     fit$id <- as.vector(id)
  384.     fit$max.id <- z$mc
  385.     z$wcor <- matrix(z$wcor, ncol = maxclsz)
  386.     fit$working.correlation <- z$wcor
  387.     fit$scale <- z$sc
  388.     fit$robust.variance <- z$rv
  389.     fit$naive.variance <- z$nv
  390.     fit$xnames <- xnames
  391.     fit$error <- z$err
  392.     dimnames(fit$robust.variance) <- list(xnames, xnames)
  393.     dimnames(fit$naive.variance) <- list(xnames, xnames)
  394.     class(fit)<-c("gee", "glm")
  395.     fit
  396. }
  397.