R/asp.R

########## R-function: asp ##########

# For fitting adaptive semiparametric regression models.

# Last changed: 30 JUL 2007

"asp" <-   function (form, adap = TRUE, random = NULL, group = NULL, family = "gaussian", 
            spar.method = "REML", omit.missing = NULL, niter = 20, niter.var = 50,tol=1e-6, returnFit=FALSE, 
            weights = NULL,correlation=NULL,control=NULL) 
{
  require("SemiPar")

#read the formula

  spm.info <- aspmFormRead(form, omit.missing)
  design.info <- spmDesign(spm.info)

#get the design matrices

  Xb <- design.info$X
  Zb <- design.info$Z
  Wb <- cbind(Xb, Zb)
  asp.info <- NULL
  model.matrices <- list(Xb = Xb, Zb = Zb, Wb = Wb)

#define weights for grouped binomial data

  if ((family == "binomial") & (is.null(weights))) 
    weights = rep(1, nrow(Xb))
  if (!is.null(dim(spm.info$y))) {
    weights <- apply(spm.info$y, 1, sum)
    spm.info$y <- (spm.info$y)[, 1]/weights
  }

#get non-adaptive fit 

  aspm.obj <- aspm(spm.info, random, group, family, spar.method, omit.missing,correlation=correlation,weights=weights,control=control)

  if (family != "gaussian") 
    aspm.obj$fit$fitted <- NULL
  fit <- aspm.obj$fit$fitted

#redefine correlation matrix in case of correlated errors

  if(!is.null(correlation))
    y.cov <- aspm.obj$fit$sigma^2*corMatrix(aspm.obj$fit$modelStruct$corStruct)
  else
    if(!is.null(aspm.obj$fit$sigma))
      y.cov <- aspm.obj$fit$sigma^2*diag(rep(1,nrow(Xb)))
  else
    y.cov <- NULL

#get the initial estimates for the variance of random effects

  b.hat <- as.vector(unlist(aspm.obj$fit$coef$random))
  kc <- kb <- theta <- sigma.theta <- knots.theta <- adap.ind <- theta.ind <- NULL
  Wc <- list()
  stb <- 1
  l <- 0
  if (!is.null(spm.info$pen)) 
    adap.ind <- spm.info$pen$adap
  if (!is.null(spm.info$krige)) 
    adap.ind <- c(adap.ind, spm.info$krige$adap)
  if (sum(adap.ind) < 1)
    adap <- FALSE
  if ((!is.null(spm.info$pen) | !is.null(spm.info$krige)) & 
      (adap == TRUE)) {
    if (!is.null(spm.info$pen)) 
      for (l in 1:length(spm.info$pen$name)) {
        knot <- spm.info$pen$knots[[l]]
        var.knots <- spm.info$pen$var.knot[[l]]
        theta.frame <- list(knot = knot, var.knots = var.knots)
        attach(theta.frame)
        if (!is.null(var.knots)) 
          design.info.theta <- spmDesign(aspmFormRead(as.formula(paste("1~f(knot,knots=var.knots)")), omit.missing))
        else design.info.theta <- spmDesign(aspmFormRead(as.formula(paste("1~f(knot)")), omit.missing))
        detach(theta.frame)
        Xc <- design.info.theta$X
        Zc <- design.info.theta$Z
        Wc[[l]] <- cbind(Xc, Zc)
        kc <- c(kc, ncol(design.info.theta$X), ncol(design.info.theta$Z))
        kb <- c(kb, nrow(design.info.theta$Z))
        all <- rep(1, nrow(Zc))
        assign("Xc", Xc)
        assign("Zc", Zc)
        bb <- b.hat[stb:(kb[l] + stb - 1)]
        b.hat.var <- log((bb - mean(bb))^2/var(bb))
        b.fit <- lme(b.hat.var ~ Xc - 1, random = list(all = pdIdent(~Zc - 1)))
        theta <- c(theta, c(as.vector(b.fit$coef$fixed), as.vector(unlist(b.fit$coef$random))))
        theta.ind <- c(theta.ind, rep(adap.ind[l], ncol(Wc[[l]])))
        sigma.theta <- c(sigma.theta,  as.vector(b.fit$sigma^2 * exp(2 * unlist(b.fit$modelStruct$reStruct))))
        knots.theta <- c(knots.theta, list(design.info.theta$spm.info$pen$knots[[1]]))
        stb <- stb + kb[l]
      }
    if (!is.null(spm.info$krige)) {
      knot <- spm.info$krige$knots
      var.knots <- spm.info$krige$var.knot
      knots1 <- knot[, 1]
      knots2 <- knot[, 2]
      theta.frame <- list(knots1 = knots1, knots2 = knots2, 
                          var.knots = var.knots)
      attach(theta.frame)
      if (!is.null(var.knots)) 
        design.info.theta <- spmDesign(aspmFormRead(as.formula(paste("1~f(knots1,knots2,knots=var.knots)")), omit.missing))
      else design.info.theta <- spmDesign(aspmFormRead(as.formula(paste("1~f(knots1,knots2)")), omit.missing))
      Xc <- design.info.theta$X
      Zc <- design.info.theta$Z
      Wc[[l + 1]] <- cbind(Xc, Zc)
      kc <- c(kc, ncol(design.info.theta$X), ncol(design.info.theta$Z))
      kb <- c(kb, nrow(design.info.theta$Z))
      all <- rep(1, nrow(Zc))
      assign("Xc", Xc)
      assign("Zc", Zc)
      bb <- b.hat[stb:(kb[l + 1] + stb - 1)]
      b.hat.var <- log((bb - mean(bb))^2/var(bb))
      b.fit <- lme(b.hat.var ~ Xc - 1, random = list(all = pdIdent(~Zc - 1)))
      theta <- c(theta, c(as.vector(b.fit$coef$fixed), as.vector(unlist(b.fit$coef$random))))
      theta.ind <- c(theta.ind, rep(adap.ind[l + 1], length(theta)))
      sigma.theta <- c(sigma.theta, as.vector(b.fit$sigma^2 * exp(2 * unlist(b.fit$modelStruct$reStruct))))
      knots.theta <- c(knots.theta, list(design.info.theta$spm.info$krige$knots[[1]]))
      detach(theta.frame)
    }


#store the model matrices
    
    Wc <- bdiag(Wc)
  if (family != "gaussian")
    {
      if (family == "poisson") 
      V <- c(fit)
    if (family == "binomial") 
      V <- c(fit * (1 - fit) * weights)
    V12 <- sqrt(V)
    Zb <- t(t(Zb) * V12)
    }
    ZVZ <- t(Zb)%*%Zb
    XVX <- t(Xb)%*%Xb
    ZVX <- t(Zb)%*%Xb
    
    model.matrices <- list(Xb = Xb, Zb = Zb, Wc = Wc, ZVZ=ZVZ, XVX=XVX, ZVX=ZVX, kb = kb, kc = kc)

#iterate for the next estimate of the variance of random effects

    for (i in 1:niter) {
      theta.info <- list(theta = theta, sigma.theta = sigma.theta, asp.info = aspm.obj, model.matrices = model.matrices)

	#get the new estimate for the spline coefficients

      theta.obj <- aspGetTheta(theta.info, niter.var, family, weights,spar.method)

	#obtain variance of random effects

      sigma.theta <- theta.obj$sigma.theta
      theta <- theta.obj$theta
      theta1 <- c(theta)
      theta1[!theta.ind] <- 0
      ran.var <- exp(c(Wc %*% theta1))
      
	#obtain remaining estimate with the estimated variance of random effects

      aspm.obj <- aspm(spm.info, random, group, family, spar.method, omit.missing, Si.b = ran.var, weights = weights,correlation=correlation,control=control)
      if (family != "gaussian") 
        aspm.obj$fit$fitted <- NULL
      fit1 <- aspm.obj$fit$fitted
      epsilon.fit <- sum((fit - fit1)^2)/sum(fit^2)
      fit <- fit1

	#exit if epsilon is small

      if (epsilon.fit <= tol) 
        break
      if (i == niter)
        
          if (returnFit==TRUE) {cat("Warning: Iteration limit reached without convergence");break}
          else
            stop ("Iteration limit reached without convergence")
    }

#extract objects for the output related to the adaptive fit

    theta <- theta1
    subknots.pen <- knots.theta
    subknots.krige <- design.info.theta$info$krige$knots
    random.var.fix <- as.vector(aspm.obj$fit$sigma^2 * exp(2 * unlist(aspm.obj$fit$modelStruct$reStruct)))
    coeff.random.var <- c(theta) + c(ginv(Wc) %*% rep(log(random.var.fix), kb))
    model.matrices <- list(Xb = Xb, Zb = Zb, Wb = cbind(Xb, Zb), Wc = Wc)
    if(!is.null(correlation))
      y.cov <- aspm.obj$fit$sigma^2*corMatrix(aspm.obj$fit$modelStruct$corStruct)
    else
      y.cov <- aspm.obj$fit$sigma^2*diag(rep(1,nrow(Xb)))
    asp.info <- list(subknots = list(subknots.pen = subknots.pen, subknots.krige = subknots.krige), coef.random.var = coeff.random.var,  var.random.var = sigma.theta)
  }

#extract remaining objects for the output

  fitted <- aspm.obj$fit$fitted
  coeff.mean.fixed <- aspm.obj$fit$coef$fixed
  coeff.mean.random <- as.vector(unlist(aspm.obj$fit$coef$random))
  lin.x <- aspm.obj$info$lin$x
  pen.x <- aspm.obj$info$pen$x
  krige.x <- aspm.obj$info$krige$x
  pen.knots <- aspm.obj$info$pen$knots
  krige.knots <- aspm.obj$info$krige$knots
  random.var <- diag(aspm.obj$aux$random.var)

#construct ouput

  asp.info <- c(list(fitted = fitted, coef.mean = c(coeff.mean.fixed, coeff.mean.random), design.matrices = model.matrices, 
                     x = list(lin.x = lin.x, pen.x = pen.x, krige.x = krige.x), 
                     knots = list(pen.knots = pen.knots, krige.knots = krige.knots), 
                     y.cov=y.cov, random.var = random.var),asp.info)

  asp.fit.obj <- c(asp.info, aspm.obj)
  class(asp.fit.obj) <- "spm"
  return(asp.fit.obj)
}

Try the AdaptFit package in your browser

Any scripts or data that you put into this service are public.

AdaptFit documentation built on May 1, 2019, 9:45 p.m.