Skip to content

[R-pkg-devel] Errors in R package - Updated

2 messages · Steven Spiriti, Duncan Murdoch

#
To Whom It May Concern:

    I have created a package called "freeknotsplines" in R, and I need to
make a few updates to it.  I resubmitted the package, and received the
following error message:

checking examples ... ERROR
Running examples in ?freeknotsplines-Ex.R? failed
The error most likely occurred in:
Error: $ operator not defined for this S4 class
Execution halted

     This error occurs in one of the examples in the documentation.   It is
platform dependent, and does not occur on the machine I am using to create
the package.  freekt is a new class, which is apparently an S4 class.

Here is the R code for my package:

#' @export coef freekt
#' @export fitted freekt
#' @export residuals freekt
#' @export plot freekt
#' @export summary freekt
#' @export AIC freekt
#' @export AICc freekt
#' @export BIC freekt
#' @export adjGCV freekt
#' @export adjAIC freekt

library(splines)

setClass("freekt", slots = c(x = "numeric", y = "numeric", degree =
"integer",
    seed = "integer", stream = "integer", lambda = "numeric",
    optknot = "numeric", tracehat = "numeric", GCV = "numeric",
    GSJS = "numeric", call = "call"))

freelsgen <- function(x, y, degree, numknot, seed = 5, stream = 0)
{
  n <- length(x)
  ord <- degree + 1
  optknot <- rep(0, times = numknot)
  tracehat <- 0
  GCV <- 0
  GSJS <- 0
  result <- .C("freelsgen", as.integer(n), as.double(x), as.double(y),
               as.integer(ord), as.integer(numknot), as.integer(seed),
               as.integer(stream), as.double(optknot),
               as.double(tracehat), as.double(GCV), as.double(GSJS))
  answer <- new("freekt", x = x, y = y, degree = as.integer(degree),
                seed = as.integer(seed), stream = as.integer(stream),
                lambda = 0, optknot = result[[8]],
                tracehat = result[[9]], GCV = result[[10]],
                GSJS = result[[11]], call = match.call())
  return(answer)
}

freelsgold <- function(x, y, degree, numknot, seed = 5, stream = 0)
{
  n <- length(x)
  ord <- degree + 1
  optknot <- rep(0, times = numknot)
  tracehat <- 0
  GCV <- 0
  GSJS <- 0
  result <- .C("freelsgold", as.integer(n), as.double(x), as.double(y),
               as.integer(ord), as.integer(numknot), as.integer(seed),
               as.integer(stream), as.double(optknot),
               as.double(tracehat), as.double(GCV), as.double(GSJS))
  answer <- new("freekt", x = x, y = y, degree = as.integer(degree),
                seed = as.integer(seed), stream = as.integer(stream),
                lambda = 0, optknot = result[[8]],
                tracehat = result[[9]], GCV = result[[10]],
                GSJS = result[[11]], call = match.call())
  return(answer)
}

freepsgen <- function(x, y, degree, numknot, seed = 5, stream = 0)
{
  n <- length(x)
  ord <- degree + 1
  optknot <- rep(0, times = numknot)
  lambda <- 0
  tracehat <- 0
  GCV <- 0
  GSJS <- 0
  result <- .C("freepsgen", as.integer(n), as.double(x), as.double(y),
               as.integer(ord), as.integer(numknot), as.integer(seed),
               as.integer(stream), as.double(lambda),
               as.double(optknot), as.double(tracehat),
               as.double(GCV), as.double(GSJS))
  answer <- new("freekt", x = x, y = y, degree = as.integer(degree),
               seed = as.integer(seed), stream = as.integer(stream),
               lambda = result[[8]], optknot = result[[9]],
               tracehat = result[[10]], GCV = result[[11]],
               GSJS = result[[12]], call = match.call())
  return(answer)
}

freepsgold <- function(x, y, degree, numknot, seed = 5, stream = 0)
{
  n <- length(x)
  ord <- degree + 1
  optknot <- rep(0, times = numknot)
  lambda <- 0
  tracehat <- 0
  GCV <- 0
  GSJS <- 0
  result <- .C("freepsgold", as.integer(n), as.double(x), as.double(y),
               as.integer(ord), as.integer(numknot), as.integer(seed),
               as.integer(stream), as.double(lambda),
               as.double(optknot), as.double(tracehat),
               as.double(GCV), as.double(GSJS))
  answer <- new("freekt", x = x, y = y,  degree = as.integer(degree),
               seed = as.integer(seed), stream = as.integer(stream),
               lambda = result[[8]], optknot = result[[9]],
               tracehat = result[[10]], GCV = result[[11]],
               GSJS = result[[12]], call = match.call())
  return(answer)
}

chgbasismat <- function(knot, ord)
{
  dimmat <- length(knot) - ord
  answer <- matrix(0, nrow = dimmat, ncol = dimmat)
  for (j in 0:(ord-1))
  {
      brow <- splineDesign(knot, knot[1], ord, j)
      brow <- as.vector(brow/factorial(j))
      answer[j + 1, ] <- brow
  }
  nknot <- sort(-1*knot)
  for (j in 1:(dimmat - ord))
  {
      brow <- splineDesign(knot, knot[ord + j], ord, ord - 1)
      brow2 <- splineDesign(nknot, nknot[length(knot) - ord - (j - 1)],
               ord, ord - 1)
      brow2 <- brow2[dimmat:1]
      brow <- brow + (-1)^ord * brow2
      brow <- as.vector(brow/factorial(ord - 1))
      answer[ord + j, ] <- brow
  }
  return(answer)
}

coef.freekt <- function(object, ...)
{
  xdat <- object at x
  ydat <- object at y
  optknot <- object at optknot
  ord <- object at degree + 1
  lambda <- object at lambda
  fulloptknot <- c(rep(min(xdat), ord), optknot, rep(max(xdat), ord))  #
includes endpoints
  Xmat <- splineDesign(fulloptknot, xdat, ord)
  if ((lambda == 0) | (length(optknot) == 0))
    coef <- solve(t(Xmat)%*%Xmat, t(Xmat)%*%ydat)
  else
    {
      numknots <- length(optknot)
      Amat <- chgbasismat(fulloptknot, ord)
      Istar <- diag(c(rep(0, times = ord), rep(1, times = numknots)))
      coef <- solve(t(Xmat)%*%Xmat + lambda*t(Amat)%*%Istar%*%Amat,
                    t(Xmat)%*%ydat)
    }
  return(coef)
}

fitted.freekt <- function(object, xfit = object at x, ...)
{
  xdat <- object at x
  ydat <- object at y
  optknot <- object at optknot
  ord <- object at degree + 1
  fulloptknot <- c(rep(min(xdat), ord), optknot, rep(max(xdat), ord))  #
includes endpoints
  coef <- coef.freekt(object)
  yfit <- splineDesign(fulloptknot, xfit, ord) %*%coef
  return(yfit)
}

residuals.freekt <- function(object, ...)
{
  fit <- fitted.freekt(object)
  return(object at y - fit)
}

plot.freekt <- function(x, xfit = x at x, linecolor="blue", lwd = 1, lty = 1,
...)
{
  xfit <- as.vector(xfit)
  yfit <- fitted.freekt(x, xfit)
  plot(x at x, x at y, ...)
  lines(xfit[order(xfit)], yfit[order(xfit)], col=linecolor,
        lwd = lwd, lty = lty)
}

summary.freekt <- function(object, ...)
{
  answer <- NULL
  if (object at lambda != 0)
  {
     currline <- c(object at lambda,
        rep(NA, times = length(object at optknot)-1))
     answer <- rbind(answer, currline)
  }
  currline <- object at optknot
  answer <- rbind(answer, currline)
  currline <- c(object at GCV,
        rep(NA, times = length(object at optknot)-1))
  answer <- rbind(answer, currline)
  RSS <- sum((residuals(object))^2)
  currline <- c(RSS, rep(NA, times = length(object at optknot)-1))
  answer <- rbind(answer, currline)
  if (object at lambda != 0)
       rownames(answer) <-
           c("Optimal lambda", "Optimal knots", "GCV", "RSS")
  else
       rownames(answer) <- c("Optimal knots", "GCV", "RSS")
  colnames(answer) <- rep("", times = length(object at optknot))
  class(answer) <- "table"
  print(answer)
}

AIC.freekt <- function(object, ..., k = 2)
{
  answer <- 0
  RSS <- sum((residuals(object))^2)
  n <- length(object at x)
  npar <- object at tracehat
  answer <- n * log(RSS/n) + k * npar
  return(answer)
}

AICc.freekt <- function(object)
{
  answer <- 0
  RSS <- sum((residuals(object))^2)
  n <- length(object at x)
  npar <- object at tracehat
  answer <- n * log(RSS/n) + 2 * npar +
            2 * npar * (npar + 1) /(n - npar - 1)
  return(answer)
}

BIC.freekt <- function(object, ...)
{
  answer <- 0
  RSS <- sum((residuals(object))^2)
  n <- length(object at x)
  npar <- object at tracehat
  answer <- n * log(RSS/n) + log(n) * npar
  return(answer)
}

adjGCV.freekt <- function(object, d = 3)
{
  RSS <- sum((residuals(object))^2)
  n <- length(object at x)
  adjtrace <- object at tracehat + d * length(object at optknot)
  answer <- (RSS/n) / (1 - adjtrace / n)^2
  return(answer)
}

adjAIC.freekt <- function(object)
{
  answer <- 0
  RSS <- sum((residuals(object))^2)
  n <- length(object at x)
  npar <- object at tracehat
  effdim <- 2 * npar - object at degree - 1
  answer <- n * log(RSS/n) +  2 * effdim
  return(answer)
}

fit.search.numknots <- function(x, y, degree, minknot = 1, maxknot = 5,
                                alg = "LS", search = "genetic",
                                knotnumcrit = "adjGCV", k = 2,
                                d = 3, seed = 5, stream = 0)
{
  bestcrit <- Inf
  funcname <- ""
  answer <- NULL
  if ((alg == "LS") && (search == "genetic"))
    funcname <- "freelsgen"
  if ((alg == "LS") && (search == "golden"))
    funcname <- "freelsgold"
  if ((alg == "PS") && (search == "genetic"))
    funcname <- "freepsgen"
  if ((alg == "PS") && (search == "golden"))
    funcname <- "freepsgold"
  for (numknot in seq(from = minknot, to = maxknot))
  {
      currcall <- call(funcname, x, y, degree, numknot, seed, stream)
      currfit <- eval(currcall)
      currcrit <- switch(knotnumcrit, GCV = currfit at GCV,
          AIC = AIC(currfit, k = k), adjAIC = adjAIC.freekt(currfit),
          AICc = AICc.freekt(currfit), BIC = BIC(currfit),
          adjGCV = adjGCV.freekt(currfit, d))
      print(paste("Number of knots = ", numknot, ", ", knotnumcrit,
          " = ", currcrit, sep = ""), quote = FALSE)
      if (currcrit < bestcrit)
      {
          bestcrit <- currcrit
          answer <- currfit
      }
  }
  return(answer)
}

     I am still not clear on what is causing this error.

Sincerely,

Steven Spiriti
#
On 25/05/2018 3:25 PM, Steven Spiriti wrote:
The problem is that you are mixing up S3 and S4 stuff.  You have an S4 
class "freekt", but you are defining functions like coef.freekt and 
expecting them to act as methods for that class.  That's not how to 
define methods in S4.  You use the methods::setMethod() function to do 
that.  Using the naming convention generic.class() is the way to do it 
in S3.

Since your S4 class has no "coef" methods, you'll get the old S3 generic 
handling your coef(xy.freekt) call, and it doesn't know what to do with 
it.  It will pass the object to stats:::coef.default, which tries to 
extract the coefficients using object$coefficients, and that fails when 
object is xy.freekt.

Duncan Murdoch