R/logit.survey.R

In IQSS/Zelig4: Everyone's Statistical Software

#' @export
zelig2logit.survey <- function(
                               formula,
                               weights=NULL, 
                               ids=NULL,
                               probs=NULL,
                               strata = NULL,  
                               fpc=NULL,
                               nest = FALSE,
                               check.strata = !nest,
                               repweights = NULL,
                               type,
                               combined.weights=FALSE,
                               rho = NULL,
                               bootstrap.average=NULL, 
                               scale=NULL,
                               rscales=NULL,
                               fpctype="fraction",
                               return.replicates=FALSE,
                               na.action="na.omit",
                               start=NULL,
                               etastart=NULL, 
                               mustart=NULL,
                               offset=NULL, 	      		
                               model1=TRUE,
                               method="glm.fit",
                               x=FALSE,
                               y=TRUE,
                               contrasts=NULL,
                               design=NULL,
                               data
                               ) {

  loadDependencies("survey")

  if (is.null(ids))
    ids <- ~1

  # the following lines designate the design
  # NOTE: nothing truly special goes on here;
  #       the below just makes sure the design is created correctly
  #       for whether or not the replication weights are set
  design <- if (is.null(repweights))
    svydesign(
              data=data,
              ids=ids,
              probs=probs,
              strata=strata,
              fpc=fpc,
              nest=nest,
              check.strata=check.strata,
              weights=weights
              )

  else {
    .survey.prob.weights <- weights
    
    svrepdesign(
                data=data,
                repweights=repweights, 	
                type=type,
                weights=weights,
                combined.weights=combined.weights, 
                rho=rho,
                bootstrap.average=bootstrap.average,
                scale=scale,
                rscales=rscales,
                fpctype=fpctype,
                fpc=fpc
                )
  }

  # we cannot plug in family=Gamma yet because of weird issues
  # with glm. Uncomment the below lines for an explanation:

  ## fails:
  # test <- Gauss
  # svyglm(formula=formula, design=design, family=test)

  ## works:
  # svyglm(formula=formula, design=design, family=Gauss)

  # this is because of how glm is written (it evaluates the
  # family variable as a function in the parent.frame)

  z(.function = svyglm,
    formula = formula,
    design  = design,
    family  = quasibinomial(link="logit")
    )
}

#' @S3method param logit.survey
param.logit.survey <- function(obj, num=1000, ...) {
  list(
       simulations = mvrnorm(num, coef(obj), vcov(obj)),
       alpha = NULL,
       fam   = binomial(link="logit")
       )
}

#' @S3method qi logit.survey
qi.logit.survey <- function(obj, x, x1=NULL, y=NULL, num=1000, param=NULL) {

  model <- GetObject(obj)

  coef <- coef(param)
  alpha <- alpha(param)

  eta <- coef %*% t(x)

  link.inverse <- linkinv(param)

  theta <- matrix(link.inverse(eta), nrow=nrow(coef))

  pr <- ev <- matrix(NA, nrow=nrow(theta), ncol(theta))

  dimnames(pr) <- dimnames(ev) <- dimnames(theta)

  ev <- theta

  for (k in 1:ncol(theta)) {
    pr[,k] <- rbinom(length(ev[,k]), 1, ev[,k])
    pr[,k] <- as.character(pr[,k])
  }

  levels(pr) <- c("0", "1")
  
  if (!is.null(y) && NCOL(y))
    y <- y[,1]


  # invisiblify 
  pr1 <- ev1 <- fd <- rr <- NA

  
  if (!is.null(x1)) {
    ev1 <- theta1 <- matrix(link.inverse(coef %*% t(x1)),
                            nrow = nrow(coef)
                            )


    pr1 <- matrix(NA, nrow=nrow(theta), ncol(theta))

    for (k in 1:ncol(theta)) {
      pr1[,k] <- rbinom(length(ev1[,k]), 1, ev1[,k])
      pr1[,k] <- as.character(pr1[,k])
    }

    levels(pr1) <- c("0", "1")
    
    fd <- ev1-ev
    rr <- ev1/ev
  }


  att.ev <- att.pr <- NA

  if (!is.null(y)) {

    yvar <- matrix(rep(y, nrow(coef)),
                   nrow = nrow(coef)
                   )

    tmp.ev <- yvar - ev
    tmp.pr <- yvar - as.integer(pr)

    att.ev <- matrix(apply(tmp.ev, 1, mean), nrow=nrow(coef))
    att.pr <- matrix(apply(tmp.pr, 1, mean), nrow=nrow(coef))
  }

  list(
       "Expected Values: E(Y|X)" = ev,
       "Expected Values: E(Y|X1)" = ev1,
       "Predicted Values: Y|X" = pr,
       "Predicted Values: Y|X1" = pr1,
       "First Differences: E(Y|X1) - E(Y|X)" = fd,
       "Risk Ratios: P(Y=1|X1)/P(Y=0|X)" = rr,
       "Average Treatment Effect: Y - EV" = att.ev,
       "Average Treatment Effect: Y - PR" = att.pr
       )
}

#' @S3method describe logit.survey
describe.logit.survey <- function(...) {
  list(
       authors = "Nicholas Carnes",
       year = 2008,
       description = "Survey-Weighted Logitistic Regression for Continuous, Positive Dependent Variables"
       )
}