Nothing
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
cao.fit <-
function(x, y, w = rep_len(1, length(x[, 1])),
etastart = NULL, mustart = NULL, coefstart = NULL,
offset = 0, family,
control = cao.control(...), criterion = "coefficients",
qr.arg = FALSE, constraints = NULL, extra = NULL,
Terms = Terms, function.name = "cao", ...) {
maxitl <- NULL
fv <- NULL
eff.n <- nrow(x) # + sum(abs(w[1:nrow(x)]))
specialCM <- NULL
post <- list()
check.rank <- TRUE
nonparametric <- TRUE
optim.maxit <- control$optim.maxit
save.weights <- control$save.weights
trace <- control$trace
minimize.criterion <- control$min.criterion
n <- dim(x)[1]
copy.X.vlm <- FALSE # May be overwritten in @initialize
X.vlm.save <- NULL
intercept.only <- ncol(x) == 1 && dimnames(x)[[2]] == "(Intercept)"
y.names <- predictors.names <- NULL # May be overwritten in @initialize
n.save <- n
Rank <- control$Rank
rrcontrol <- control #
if (length(family@initialize))
eval(family@initialize) # Initialize mu and M (and optionally w)
n <- n.save
modelno <- switch(family@vfamily[1], "poissonff" = 2,
"binomialff" = 1, "quasipoissonff" = 0,
"quasibinomialff" = 0, "negbinomial" = 3,
"gamma2" = 5, "gaussianff" = 8,
0) # stop("cannot fit this model using fast algorithm")
if (!modelno)
stop("the family function does not work with cao()")
if (modelno == 1)
modelno <- get("modelno", envir = VGAMenv)
eval(rrr.init.expression)
if (length(etastart)) {
eta <- etastart
mu <- if (length(mustart)) mustart else
if (length(body(slot(family, "linkinv"))))
slot(family, "linkinv")(eta, extra) else
warning("argument 'etastart' assigned a value ",
"but there is no 'linkinv' slot to use it")
}
if (length(mustart)) {
mu <- mustart
if (length(body(slot(family, "linkfun")))) {
eta <- slot(family, "linkfun")(mu, extra)
} else {
warning("argument 'mustart' assigned a value ",
"but there is no 'link' slot to use it")
}
}
M <- if (is.matrix(eta)) ncol(eta) else 1
if (length(family@constraints))
eval(family@constraints)
special.matrix <- matrix(-34956.125, M, M) # An unlikely used matrix
just.testing <- cm.VGAM(special.matrix, x, rrcontrol$noRRR, constraints)
findex <- trivial.constraints(just.testing, special.matrix)
if (all(findex == 1))
stop("No covariates to form latent variables from.")
colx1.index <- names.colx1.index <- NULL
dx2 <- dimnames(x)[[2]]
if (sum(findex)) {
asx <- attr(x, "assign")
for (ii in names(findex))
if (findex[ii]) {
names.colx1.index <- c(names.colx1.index, dx2[asx[[ii]]])
colx1.index <- c(colx1.index, asx[[ii]])
}
names(colx1.index) <- names.colx1.index
}
rrcontrol$colx1.index <- control$colx1.index <- colx1.index
colx2.index <- 1:ncol(x)
names(colx2.index) <- dx2
colx2.index <- colx2.index[-colx1.index]
p1 <- length(colx1.index)
p2 <- length(colx2.index)
rrcontrol$colx2.index <- control$colx2.index <- colx2.index
Cmat <- if (length(rrcontrol$Cinit)) {
matrix(rrcontrol$Cinit, p2, Rank)
} else {
if (!rrcontrol$Use.Init.Poisson.QO) {
matrix(rnorm(p2 * Rank, sd = rrcontrol$sd.Cinit), p2, Rank)
} else {
.Init.Poisson.QO(ymat = as.matrix(y),
X1 = x[, colx1.index, drop = FALSE],
X2 = x[, colx2.index, drop = FALSE],
Rank = rrcontrol$Rank, trace = rrcontrol$trace,
max.ncol.etamat = rrcontrol$Etamat.colmax,
Crow1positive = rrcontrol$Crow1positive,
constwt = any(family@vfamily[1] ==
c("negbinomial", "gamma2", "gaussianff")),
takelog = any(family@vfamily[1] != c("gaussianff")))
}
}
rrcontrol$Cinit <- control$Cinit <- Cmat # Good for valt0()
Hlist <- process.constraints(constraints, x, M, specialCM = specialCM)
nice31 <- checkCMCO(Hlist, control = control, modelno = modelno)
if (nice31 != 1)
stop("not nice")
ncolHlist <- unlist(lapply(Hlist, ncol))
latvar.mat <- x[, colx2.index, drop = FALSE] %*% Cmat
rmfromVGAMenv(c("etamat", "beta"), prefix = ".VGAM.CAO.")
Nice21 <- length(names.colx1.index) == 1 &&
names.colx1.index == "(Intercept)"
if (!Nice21)
stop("'noRRR = ~ 1' is supported only, without constraints")
NOS <- ifelse(modelno %in% c(3, 5), M/2, M)
p1star. <- if (Nice21) ifelse(modelno %in% c(3, 5), 2, 1) else M
p2star. <- if (Nice21) Rank else stop("not Nice21")
pstar. <- p1star. + p2star.
nstar <- if (Nice21) ifelse(modelno %in% c(3, 5), n * 2, n) else n * M
lenbeta <- pstar. * ifelse(Nice21, NOS, 1)
othint <-
c(Rank, control$eq.tol, pstar. ,
dim2wz = 1, inited = 0, # w(, dimw) cols
modelno, maxitl = control$maxitl,
actnits = 0, twice = 0, p1star. ,
p2star. , Nice21, lenbeta,
controlI.tolerances = 0, control$trace,
p1, p2 = p2, imethod = control$imethod, bchat = 0)
othdbl <- c(small = control$SmallNo, fseps = control$epsilon,
.Machine$double.eps,
iKvector = rep_len(control$iKvector, NOS),
iShape = rep_len(control$iShape, NOS),
resss = 0, bfeps = control$bf.epsilon, hstep = 0.1)
for (iter in 1:optim.maxit) {
if (control$trace) {
cat("\nIteration", iter, "\n")
flush.console()
}
conjgrad <- optim(par = c(Cmat), fn = callcaoc,
gr = if (control$GradientFunction) calldcaoc else NULL,
method = "BFGS",
control = list(fnscale = 1,
trace = as.integer(control$trace),
maxit = control$Maxit.optim,
REPORT = 10),
etamat = eta, xmat = x, ymat = y, # as.matrix(y),
wvec = w, modelno = modelno,
Control = control,
Nice21 = Nice21,
p1star. = p1star. , p2star. = p2star. ,
n = n, M = M,
othint = othint, othdbl = othdbl,
alldump = FALSE)
Cmat <- matrix(conjgrad$par, p2, Rank) # old becoz of scale(cmatrix)
if (converged <- (conjgrad$convergence == 0))
break
}
if (!converged) {
if (control$maxitl > 1) {
warning("convergence not obtained in ", control$maxitl,
" iterations.")
} else {
warning("convergence not obtained")
}
} else {
}
Cmat <- crow1C(Cmat, control$Crow1positive) # Make sure signs are right
flush.console()
temp9 <- callcaoc(cmatrix = Cmat,
etamat = eta, xmat = x, ymat = y,
wvec = w, modelno = modelno,
Control = control,
Nice21 = Nice21,
p1star. = p1star. , p2star. = p2star. ,
n = n, M = M,
othint = othint, othdbl = othdbl,
alldump = TRUE)
if (!is.list(extra))
extra <- list()
extra$Cmat <- temp9$Cmat
ynames <- dimnames(y)[[2]]
extra$df1.nl <- temp9$df1.nl
extra$lambda1 <- temp9$lambda1
extra$spar1 <- temp9$spar1
names(extra$df1.nl) <-
names(extra$lambda1) <-
names(extra$spar1) <- ynames
if (Rank == 2) {
extra$spar2 <- temp9$spar2
extra$lambda2 <- temp9$lambda2
extra$df2.nl <- temp9$df2.nl
names(extra$df2.nl) <-
names(extra$lambda2) <-
names(extra$spar2) <- ynames
}
extra$alldeviance <- temp9$alldeviance
names(extra$alldeviance) <- ynames
mu <- matrix(temp9$fitted, n, NOS, byrow = TRUE)
dn <- labels(x)
yn <- dn[[1]]
if (is.matrix(mu)) {
if (length(dimnames(y)[[2]])) {
y.names <- dimnames(y)[[2]]
}
if (length(dimnames(mu)[[2]])) {
y.names <- dimnames(mu)[[2]]
}
dimnames(mu) <- list(yn, y.names)
} else {
names(mu) <- names(fv)
}
fit <- list(
fitted.values = mu,
Cmatrix = Cmat,
terms = Terms) # terms: This used to be done in vglm()
misc <- list(
criterion = criterion,
intercept.only = intercept.only,
predictors.names = predictors.names,
M = M,
n = n,
nonparametric = nonparametric,
p = ncol(x),
ynames = ynames)
crit.list <- list()
crit.list$deviance <- temp9$deviance
if (w[1] != 1 || any(w != w[1]))
fit$prior.weights <- w
if (length(family@last))
eval(family@last)
structure(c(fit,
temp9,
list(
contrasts = attr(x, "contrasts"),
control = control,
crit.list = crit.list,
extra = extra,
family = family,
iter = iter,
misc = misc,
post = post,
x = x,
y = y)),
vclass = family@vfamily)
}
cao.control <- function(Rank = 1,
all.knots = FALSE,
criterion = "deviance",
Cinit = NULL,
Crow1positive = TRUE,
epsilon = 1.0e-05,
Etamat.colmax = 10,
GradientFunction = FALSE, # For now 20041224
iKvector = 0.1,
iShape = 0.1,
noRRR = ~ 1,
Norrr = NA,
SmallNo = 5.0e-13,
Use.Init.Poisson.QO = TRUE,
Bestof = if (length(Cinit)) 1 else 10,
maxitl = 10, # was 40 prior to 20100420
imethod = 1,
bf.epsilon = 1.0e-7,
bf.maxit = 10, # was 40 prior to 20100420
Maxit.optim = 250,
optim.maxit = 20,
sd.sitescores = 1.0,
sd.Cinit = 0.02,
suppress.warnings = TRUE,
trace = TRUE,
df1.nl = 2.5, # About 1.5--2.5 gives flexibility of a quadratic
df2.nl = 2.5, # About 1.5--2.5 gives flexibility of a quadratic
spar1 = 0, # 0 means df1.nl is used
spar2 = 0, # 0 means df2.nl is used
...) {
if (length(Norrr) != 1 || !is.na(Norrr)) {
warning("argument 'Norrr' has been replaced by 'noRRR'. ",
"Assigning the latter but using 'Norrr' will become ",
"an error in the next VGAM version soon.")
noRRR <- Norrr
}
if (!is.Numeric(iShape, positive = TRUE))
stop("bad input for argument 'iShape'")
if (!is.Numeric(iKvector, positive = TRUE))
stop("bad input for argument 'iKvector'")
if (!is.Numeric(imethod, positive = TRUE, length.arg = 1,
integer.valued = TRUE))
stop("bad input for argument 'imethod'")
if (criterion != "deviance")
stop("'criterion' must be 'deviance'")
if (GradientFunction)
stop("20050114; GradientFunction = TRUE not working yet")
se.fit <- as.logical(FALSE)
if (se.fit)
stop("se.fit = FALSE handled only")
if (length(Cinit) && !is.Numeric(Cinit))
stop("Bad input for argument 'Cinit'")
if (!is.Numeric(Bestof, length.arg = 1,
integer.valued = TRUE, positive = TRUE))
stop("Bad input for argument 'Bestof'")
if (!is.Numeric(maxitl, length.arg = 1,
integer.valued = TRUE, positive = TRUE))
stop("Bad input for argument 'maxitl'")
if (!is.Numeric(bf.epsilon, length.arg = 1,
positive = TRUE))
stop("Bad input for argument 'bf.epsilon'")
if (!is.Numeric(bf.maxit, integer.valued = TRUE,
positive = TRUE, length.arg = 1))
stop("Bad input for argument 'bf.maxit'")
if (!is.Numeric(Etamat.colmax, positive = TRUE,
length.arg = 1) ||
Etamat.colmax < Rank)
stop("bad input for argument 'Etamat.colmax'")
if (!is.Numeric(Maxit.optim, integer.valued = TRUE,
positive = TRUE, length.arg = 1))
stop("Bad input for argument 'Maxit.optim'")
if (!is.Numeric(optim.maxit, length.arg = 1,
integer.valued = TRUE, positive = TRUE))
stop("Bad input for argument 'optim.maxit'")
if (!is.Numeric(sd.sitescores, length.arg = 1,
positive = TRUE))
stop("Bad input for argument 'sd.sitescores'")
if (!is.Numeric(sd.Cinit, length.arg = 1,
positive = TRUE))
stop("Bad input for argument 'sd.Cinit'")
if (!is.Numeric(df1.nl) || any(df1.nl < 0))
stop("Bad input for argument 'df1.nl'")
if (any(df1.nl >= 0 & df1.nl < 0.05)) {
warning("'df1.nl' values between 0 and 0.05 converted to 0.05")
df1.nl[df1.nl < 0.05] <- 0.05
}
if (any(df1.nl > 3.5)) {
warning("'df1.nl' values > 3.5 are excessive")
}
if (!is.Numeric(df2.nl) || any(df2.nl < 0))
stop("Bad input for argument 'df2.nl'")
if (any(df2.nl >= 0 & df2.nl < 0.05)) {
warning("'df2.nl' values between 0 and 0.05 converted to 0.05")
df2.nl[df2.nl < 0.05] <- 0.05
}
if (!is.Numeric(spar1) || any(spar1 < 0))
stop("Bad input for argument 'spar1'")
if (!is.Numeric(spar2) || any(spar2 < 0))
stop("Bad input for argument 'spar2'")
if (!is.Numeric(epsilon, positive = TRUE, length.arg = 1))
stop("Bad input for argument 'epsilon'")
if (!is.Numeric(SmallNo, positive = TRUE, length.arg = 1))
stop("Bad input for argument 'SmallNo'")
if ((SmallNo < .Machine$double.eps) ||
(SmallNo > .0001))
stop("'SmallNo' is out of range")
ans <- list(
Corner = FALSE, # A constant, not a control parameter; unneeded?
eq.tolerances = FALSE, # A constant, not a control parameter; needed
I.tolerances = FALSE, # A constant, not a control parameter; unneeded?
Quadratic = FALSE, # A constant, not a control parameter; unneeded?
all.knots = as.logical(all.knots)[1],
Bestof = Bestof,
Cinit = Cinit,
ConstrainedO = TRUE, # A constant, not a control parameter
criterion = criterion,
Crow1positive = as.logical(rep_len(Crow1positive, Rank)),
epsilon = epsilon,
Etamat.colmax = Etamat.colmax,
FastAlgorithm = TRUE, # A constant, not a control parameter
GradientFunction = as.logical(GradientFunction),
maxitl = maxitl,
bf.epsilon = bf.epsilon,
bf.maxit = bf.maxit,
imethod = imethod,
Maxit.optim = Maxit.optim,
optim.maxit = optim.maxit,
noRRR = noRRR,
Rank = Rank,
sd.sitescores = sd.sitescores,
sd.Cinit = sd.Cinit,
se.fit = se.fit, # If TRUE, then would need storage for S QR fits
SmallNo = SmallNo,
suppress.warnings = as.logical(suppress.warnings),
trace = as.integer(trace),
Use.Init.Poisson.QO = Use.Init.Poisson.QO,
iKvector = as.numeric(iKvector),
iShape = as.numeric(iShape),
DF1 = 2.5, # Used as Default value if df1.nl has no default
DF2 = 2.5, # Used as Default value if df2.nl has no default
SPAR1 = 0, # Used as Default value if spar1 has no default
SPAR2 = 0, # Used as Default value if spar2 has no default
df1.nl = df1.nl,
df2.nl = df2.nl,
spar1 = spar1,
spar2 = spar2)
ans
}
create.cms <- function(Rank = 1, M, MSratio = 1, which, p1 = 1) {
if (!is.Numeric(p1, length.arg = 1,
integer.valued = TRUE, positive = TRUE))
stop("bad input for argument 'p1'")
Hlist. <- vector("list", p1 + Rank)
for (rr in 1:(p1+Rank))
Hlist.[[rr]] <- diag(M)
names(Hlist.) <- if (p1 == 1) c("(Intercept)", names(which)) else stop()
if (MSratio == 2) {
for (r in 1:Rank)
Hlist.[[p1+r]] <- eijfun(1, M)
}
Hlist.
}
callcaoc <- function(cmatrix,
etamat, xmat, ymat, wvec, modelno,
Control, Nice21 = TRUE,
p1star. = if (modelno %in% c(3, 5)) 2 else 1,
p2star. = Rank,
n, M,
othint, othdbl,
alldump = FALSE) {
flush.console()
control <- Control
Rank <- control$Rank
p1 <- length(control$colx1.index)
p2 <- length(control$colx2.index)
yn <- dimnames(ymat)[[2]]
if (length(yn) != ncol(ymat))
stop("the column names of argument 'ymat' must be given")
queue <- qbig <- Rank # 20051019; number of smooths per species
NOS <- if (modelno %in% c(3, 5)) M/2 else M
df1.nl <- procVec(control$df1.nl, yn = yn , Default = control$DF1)
spar1 <- procVec(control$spar1, yn = yn , Default = control$SPAR1)
df2.nl <- procVec(control$df2.nl, yn = yn , Default = control$DF2)
spar2 <- procVec(control$spar2, yn = yn , Default = control$SPAR2)
if (any(c(length(spar1), length(spar2), length(df1.nl),
length(df2.nl)) != NOS))
stop("wrong length in at least one of arguments ",
"'df1.nl', 'df2.nl', 'spar1', 'spar2'")
cmatrix <- matrix(cmatrix, p2, Rank) # crow1C() needs a matrix as input
cmatrix <- crow1C(cmatrix, crow1positive = control$Crow1positive)
numat <- xmat[, control$colx2.index, drop = FALSE] %*% cmatrix
evnu <- eigen(var(numat), symmetric = TRUE)
temp7 <- if (Rank > 1) evnu$vector %*% diag(evnu$value^(-0.5)) else
evnu$vector %*% evnu$value^(-0.5)
cmatrix <- cmatrix %*% temp7
cmatrix <- crow1C(cmatrix, crow1positive = control$Crow1positive)
numat <- xmat[, control$colx2.index, drop = FALSE] %*% cmatrix
dim(numat) <- c(n, Rank)
mynames5 <- param.names("latvar", Rank, skip1 = TRUE)
nu1mat <- cbind("(Intercept)" = 1, latvar = numat)
dimnames(nu1mat) <- list(dimnames(xmat)[[1]],
c("(Intercept)", mynames5))
temp.smooth.frame <- vector("list", p1+Rank) # Temporary makeshift frame
names(temp.smooth.frame) <- c(names(control$colx1.index), mynames5)
for (uu in 1:(p1+Rank)) {
temp.smooth.frame[[uu]] <- nu1mat[, uu]
}
temp.smooth.frame <- data.frame(temp.smooth.frame)
for (uu in 1:Rank) {
attr(temp.smooth.frame[,uu+p1], "spar") <- 0 # this value unused
attr(temp.smooth.frame[,uu+p1], "df") <- 4 # this value unused
}
pstar. <- p1star. + p2star. # = Mdot + Rank
nstar <- if (Nice21) ifelse(modelno %in% c(3, 5), n * 2, n) else n * M
lenbeta <- pstar. * ifelse(Nice21, NOS, 1) # Holds the linear coeffs
inited <- if (exists(".VGAM.CAO.etamat", envir = VGAMenv)) 1 else 0
usethiseta <-
if (inited == 1)
getfromVGAMenv("etamat", prefix = ".VGAM.CAO.") else
t(etamat)
if (anyNA(usethiseta)) {
usethiseta <- t(etamat) # So that dim(usethiseta) == c(M,n)
rmfromVGAMenv("etamat", prefix = ".VGAM.CAO.")
}
usethisbeta <- if (inited == 2)
getfromVGAMenv("beta", prefix = ".VGAM.CAO.") else
double(lenbeta)
othint[5] <- inited # Refine initialization within C
pstar <- NOS * pstar.
bnumat <- if (Nice21) matrix(0, nstar, pstar.) else
stop("code not written here")
M. <- MSratio <- M / NOS # 1 or 2 usually
which <- p1 + (1:Rank) # These columns are smoothed
nwhich <- names(which) <- mynames5
origHlist <-
Hlist. <- create.cms(Rank = Rank, M = M., MSratio = MSratio,
which = which, p1 = p1) # For 1 species only
ncolHlist. <- unlist(lapply(Hlist. , ncol))
smooth.frame <- s.vam(x = nu1mat, zedd = NULL,
wz = NULL, smomat = NULL,
which = which,
smooth.frame = temp.smooth.frame,
bf.maxit = control$bf.maxit,
bf.epsilon = control$bf.epsilon,
trace = FALSE, se.fit = control$se.fit,
X.vlm.save = bnumat, Hlist = Hlist. ,
ncolHlist = ncolHlist. ,
M = M. ,
qbig = NULL, Umat = NULL, # NULL ==> unneeded
all.knots = control$all.knots, nk = NULL,
sf.only = TRUE)
ldk <- 3 * max(ncolHlist.[nwhich]) + 1 # 20020711
dimw. <- M. # Smoothing one spp. at a time
dim1U. <- M.
wz. <- matrix(0, n, dimw. )
if (names(Hlist.)[1] != "(Intercept)")
stop("something wrong here")
Hlist.[[1]] <- NULL
trivc <- rep_len(2 - M. , queue)
ncbvec <- ncolHlist.[nwhich]
ncolb <- max(ncbvec)
qbig. <- NOS * qbig # == NOS * Rank; holds all the smooths
if (!all(as.vector(ncbvec) == rep_len(1, queue)))
stop("'ncbvec' not right---should be a queue-vector of ones")
pbig <- pstar. #
contr.sp <- list(low = -1.5, ## low = 0. was default till R 1.3.x
high = 1.5,
tol = 1e-4, ## tol = 0.001 was default till R 1.3.x
eps = 2e-8, ## eps = 0.00244 was default till R 1.3.x
maxit = 500 )
npetc <-
c(n = nrow(nu1mat), p. = ncol(nu1mat), q = length(which),
se.fit = control$se.fit, 0,
control$bf.maxit, qrank = 0, M = M. , nbig = nstar, pbig = pbig,
qbig = qbig, dim2wz = dimw. , dim1U = dim1U. ,
ierror = 0, ldk = ldk,
contr.sp$maxit, iinfo = 0)
if (Rank == 2) {
smopar <- (c(spar1, spar2))[interleave.VGAM(4 * NOS, M1 = 2)]
dofvec <- (1.0 + c(df1.nl, df2.nl))[interleave.VGAM(4 * NOS, M1 = 2)]
lamvec <- 0 * dofvec
stop("20100414; havent got Rank = 2 going yet")
} else {
smopar <- c(spar1, spar2)
dofvec <- c(df1.nl, df2.nl) + 1.0
lamvec <- 0 * dofvec
}
ans1 <- .C("vcao6",
numat = as.double(numat), ymat = as.double(ymat), wvec = as.double(wvec),
etamat = as.double(usethiseta), fv = double(NOS*n), zedd = double(n*M),
wz = double(n*M), U = double(M*n), # bnumat = as.double(bnumat),
qr = double(nstar*pstar.), qraux = double(pstar.),
qpivot = integer(pstar.),
n = as.integer(n), M = as.integer(M), NOS = as.integer(NOS),
nstar = as.integer(nstar), dim1U = as.integer( M ), # for U, not U.
errcode = integer(1), othint = as.integer(othint),
deviance = double(1 + NOS), # NOS more elts added 20100413
beta = as.double(usethisbeta),
othdbl = as.double(othdbl),
npetc = as.integer(npetc), M. = as.integer( M. ),
dofvec = as.double(dofvec),
lamvec = as.double(lamvec),
smopar = as.double(smopar),
match = as.integer(smooth.frame$matcho),
as.integer(smooth.frame$nef),
which = as.integer(which),
smomat = as.double(matrix(0, n, qbig. )),
nu1mat = as.double(nu1mat),
Hlist = as.double(unlist( Hlist. )),
as.integer(ncbvec),
smap = as.integer(1:(Rank+1)), #
trivc = as.integer(trivc),
levmat = double(NOS * sum(smooth.frame$neffec * ncbvec)),
bcoefficients = double(NOS * sum(smooth.frame$nknots*ncbvec)),
xknots = as.double(unlist(smooth.frame$knots)),
bindex = as.integer(smooth.frame$bindex),
lindex = as.integer(smooth.frame$lindex),
nknots = as.integer(smooth.frame$nknots),
kindex = as.integer(smooth.frame$kindex))
flush.console()
if (ans1$errcode == 0) {
assign2VGAMenv(c("etamat", "beta"), ans1, prefix = ".VGAM.CAO.")
assign(".VGAM.CAO.cmatrix", matrix(cmatrix, p2, Rank),
envir = VGAMenv)
} else {
if (!control$suppress.warnings) {
cat("warning in callcaoc: error code = ", ans1$errcode, "\n")
cat("warning in callcaoc: npetc[14] = ", ans1$npetc[14], "\n")
flush.console()
}
rmfromVGAMenv(c("etamat", "beta"), prefix = ".VGAM.CAO.")
}
returnans <- if (alldump) {
bindex <- ans1$bindex
ncolHlist <- ncbvec
Bspline2 <- vector("list", NOS)
names(Bspline2) <- dimnames(ymat)[[2]]
Bspline <- vector("list", length(nwhich))
names(Bspline) <- nwhich
ind9 <- 0 # moving index
for (sppno in 1:NOS) {
for (ii in seq_along(nwhich)) {
ind7 <- (smooth.frame$bindex[ii]):(smooth.frame$bindex[ii+1]-1)
ans <- ans1$bcoeff[ind9+ind7]
ans <- matrix(ans, ncol = ncolHlist[nwhich[ii]])
Bspline[[ii]] <-
new(Class = "vsmooth.spline.fit",
"Bcoefficients" = ans,
"xmax" = smooth.frame$xmax[ii],
"xmin" = smooth.frame$xmin[ii],
"knots" = as.vector(smooth.frame$knots[[ii]]))
}
ind9 <- ind9 + smooth.frame$bindex[length(nwhich)+1] - 1
Bspline2[[sppno]] <- Bspline
}
qrank <- npetc[7] # Assume all species have the same qrank value
dim(ans1$etamat) <- c(M, n) # was c(n, M) prior to 20060822
df1.nl <- ans1$dofvec[1:NOS] - 1.0
lambda1 <- ans1$lamvec[1:NOS]
spar1 <- ans1$smopar[1:NOS]
if (Rank == 2) {
stop("20100414; this is not working yet")
df2.nl <- ans1$dofvec[NOS + (1:NOS)] - 1.0
lambda2 <- ans1$lamvec[NOS + (1:NOS)]
spar2 <- ans1$smopar[NOS + (1:NOS)]
}
list(deviance = ans1$deviance[1],
alldeviance = ans1$deviance[-1],
bcoefficients = ans1$bcoefficients,
bindex = ans1$bindex,
Bspline = Bspline2,
Cmat = matrix(cmatrix, p2, Rank, dimnames = list(
names(control$colx2.index), mynames5)),
coefficients = ans1$beta,
df1.nl = df1.nl,
df2.nl = if (Rank == 2) df2.nl else NULL,
df.residual = n*M - qrank - sum(ans1$df - 1),
fitted = ans1$fv, # NOS x n
kindex = ans1$kindex,
lambda1 = lambda1,
lambda2 = if (Rank == 2) lambda2 else NULL,
predictors = matrix(ans1$etamat, n, M, byrow = TRUE),
wresiduals = ans1$zedd - t(ans1$etamat), # n x M
spar1 = spar1,
spar2 = if (Rank == 2) spar2 else NULL)
} else {
ans1$deviance[1]
}
flush.console()
returnans
}
calldcaoc <- function(cmatrix,
etamat, xmat, ymat, wvec, modelno,
Control, Nice21 = TRUE,
p1star. = if (modelno %in% c(3, 5)) 2 else 1,
p2star. = Rank,
n, M,
othint, othdbl,
alldump = FALSE) {
if (alldump)
stop("really used?")
flush.console()
U <- NULL
if (!Nice21)
stop("'Nice21' must be TRUE")
control <- Control
Rank <- control$Rank
p2 <- length(control$colx2.index)
yn <- dimnames(ymat)[[2]]
if (!length( yn ))
yn <- param.names("Y", ncol(ymat))
cmatrix <- scale(cmatrix)
xmat2 <- xmat[, control$colx2.index, drop = FALSE] #ccc
numat <- xmat2 %*% matrix(cmatrix, p2, Rank)
dim(numat) <- c(nrow(xmat), Rank)
temp.smooth.frame <- vector("list", 1+Rank) # Temporary makeshift frame
mynames5 <- param.names("latvar", Rank, skip1 = TRUE)
names(temp.smooth.frame) <- c("(Intercept)", mynames5)
temp.smooth.frame[[1]] <- rep_len(1, n)
for (uu in 1:Rank) {
temp.smooth.frame[[uu+1]] <- numat[, uu]
}
temp.smooth.frame <- data.frame(temp.smooth.frame)
for (uu in 1:Rank) {
attr(temp.smooth.frame[,uu+1], "spar") <- 0 # any old value
attr(temp.smooth.frame[,uu+1], "df") <- 4 # any old value
}
pstar. <- p1star. + p2star.
nstar <- if (Nice21) ifelse(modelno %in% c(3, 5), n * 2, n) else n * M
NOS <- ifelse(modelno %in% c(3, 5), M / 2, M)
lenbeta <- pstar. * ifelse(Nice21, NOS, 1)
if (TRUE) {
inited <- if (exists(".VGAM.CAO.etamat", envir = VGAMenv))
1 else 0
usethiseta <- if (inited == 1)
get(".VGAM.CAO.etamat", envir = VGAMenv) else
t(etamat)
}
usethisbeta <- if (inited == 2)
get(".VGAM.CAO.beta", envir = VGAMenv) else
double(lenbeta)
pstar <- NOS * pstar.
bnumat <- if (Nice21)
matrix(0, nstar, pstar) else stop("need 'Nice21'")
M. <- MSratio <- M / NOS # 1 or 2 usually
p1 <- 1
which <- p1 + (1:Rank) # The first 1 is the intercept term
nwhich <- names(which) <- mynames5
origHlist <- Hlist. <-
create.cms(Rank = Rank, M = M., MSratio = MSratio,
which = which, p1 = p1) # For 1 species
ncolHlist. <- unlist(lapply(Hlist. , ncol))
nu1mat <- cbind("(Intercept)" = 1, latvar = numat)
dimnames(nu1mat) <- list(dimnames(xmat)[[1]],
c("(Intercept)", "latvar"))
smooth.frame <- s.vam(x = nu1mat, zedd = NULL, wz = NULL,
smomat = NULL, which = which,
smooth.frame = temp.smooth.frame,
bf.maxit = control$bf.maxit,
bf.epsilon = control$bf.epsilon,
trace = FALSE, se.fit = control$se.fit,
X.vlm.save = bnumat, Hlist = Hlist.,
ncolHlist = ncolHlist. ,
M = M. , qbig = NULL,
Umat = U, # NULL value ==> not needed
all.knots = control$all.knots, nk = NULL,
sf.only = TRUE)
ldk <- 4 * max(ncolHlist.[nwhich]) # was M; # Prior to 20020711
ldk <- 3 * max(ncolHlist.[nwhich]) + 1 # 20020711
wz. <- matrix(0, n, M. ) # not sure
dimw. <- if (is.matrix( wz. )) ncol( wz. ) else 1
dim1U. <- M. # 20100410
queue <- qbig <- Rank # 20051019; number of smooths per species
Hlist.[[1]] <- NULL
trivc <- rep_len(2 - M. , queue)
ncbvec <- ncolHlist.[nwhich]
ncolb <- max(ncbvec)
qbig. <- NOS * qbig # == NOS * Rank
pbig <- pstar. # Not sure
if (FALSE) {
df1.nl <- rep_len(control$df1.nl, NOS) # This is used
df2.nl <- rep_len(control$df2.nl, NOS) # This is used
spar1 <- rep_len(control$spar1, NOS) # This is used
spar2 <- rep_len(control$spar2, NOS) # This is used
} else {
df1.nl <- procVec(control$df1.nl, yn = yn , Default = control$DF1)
df2.nl <- df1.nl # 20100417; stopgap
spar1 <- procVec(control$spar1, yn = yn , Default = control$SPAR1)
spar2 <- spar1 # 20100417; stopgap
dofvec <- c(df1.nl, df2.nl)
lamvec <- 0 * dofvec
smopar <- c(spar1, spar2)
}
contr.sp <- list(low = -1.5, ## low = 0. was default till R 1.3.x
high = 1.5,
tol = 1e-4, ## tol = 0.001 was default till R 1.3.x
eps = 2e-8, ## eps = 0.00244 was default till R 1.3.x
maxit = 500 )
warning("20100405; this is old:")
npetc <-
c(n = n, p = 1+Rank, length(which), se.fit = control$se.fit, 0,
maxitl = control$maxitl, qrank = 0, M = M. , n.M = n* M. ,
pbig = sum( ncolHlist.),
qbig = qbig, dimw = dimw. , dim1U = dim1U. ,
ierror = 0, ldk = ldk)
warning("20100405; this is new:")
npetc <- c(n = nrow(nu1mat), p. = ncol(nu1mat),
q = length(which),
se.fit = control$se.fit, 0,
control$bf.maxit, qrank = 0, M = M. , nbig = nstar, pbig = pbig,
qbig = qbig, dim2wz = dimw. , dim1U = dim1U. , ierror = 0, ldk = ldk,
contr.sp$maxit, iinfo = 0)
flush.console()
if (!Nice21)
stop("need 'Nice21'")
ans1 <- .C("vdcao6",
numat = as.double(numat), as.double(ymat), as.double(wvec),
etamat = as.double(usethiseta), fv = double(NOS*n),
zedd = double(n*M),
wz = double(n*M), U = double(M*n), # bnumat = as.double(bnumat),
qr = double(nstar*pstar.), qraux = double(pstar.),
qpivot = integer(pstar.),
as.integer(n), as.integer(M), NOS = as.integer(NOS),
as.integer(nstar), dim1U = as.integer(M),
errcode = integer(1), othint = as.integer(othint),
deviance = double(1 + NOS), beta = as.double(usethisbeta),
othdbl = as.double(othdbl),
as.double(xmat2),
cmat = as.double(cmatrix),
p2 = as.integer(p2), deriv = double(p2 * Rank),
betasave = double(lenbeta),
npetc = as.integer(npetc), M. = as.integer( M. ),
dofvec = as.double(dofvec + 1.0),
lamvec = as.double(0 * dofvec),
smopar = as.double(smopar),
match = as.integer(smooth.frame$matcho),
as.integer(smooth.frame$nef),
as.integer(which),
smomat = as.double(matrix(0, n, qbig. )),
nu1mat = as.double(nu1mat),
as.double(unlist( Hlist. )),
as.integer(ncbvec), smap = as.integer(1:(Rank+1)),
trivc = as.integer(trivc),
levmat = double(NOS * sum(smooth.frame$neffec * ncbvec)),
bcoefficients = double(NOS * sum(smooth.frame$nknots * ncbvec)),
xknots = as.double(unlist(smooth.frame$knots)),
bindex = as.integer(smooth.frame$bindex),
lindex = as.integer(smooth.frame$lindex),
nknots = as.integer(smooth.frame$nknots),
kindex = as.integer(smooth.frame$kindex))
flush.console()
assign(".VGAM.CAO.etamat", ans1$etamat, envir = VGAMenv)
assign(".VGAM.CAO.z", ans1$zedd, envir = VGAMenv)
assign(".VGAM.CAO.U", ans1$U, envir = VGAMenv) # U
if (ans1$errcode == 0) {
} else {
cat("warning in calldcaoc: error code = ", ans1$errcode, "\n")
flush.console()
}
returnans <- if (alldump) {
bindex <- ans1$bindex
ncolHlist <- ncbvec
Bspline2 <- vector("list", NOS)
names(Bspline2) <- dimnames(ymat)[[2]]
Bspline <- vector("list", length(nwhich))
names(Bspline) <- nwhich
ind9 <- 0 # moving index
for (jay in 1:NOS) {
for (ii in seq_along(nwhich)) {
ind9 <- ind9[length(ind9)] + (bindex[ii]):(bindex[ii+1]-1)
ans <- ans1$bcoeff[ind9]
ans <- matrix(ans, ncol = ncolHlist[nwhich[ii]])
Bspline[[ii]] <-
new(Class = "vsmooth.spline.fit",
"Bcoefficients" = ans,
"xmax" = smooth.frame$xmax[ii],
"xmin" = smooth.frame$xmin[ii],
"knots" = as.vector(smooth.frame$knots[[ii]]))
}
Bspline2[[jay]] <- Bspline
}
qrank <- npetc[7] # Assume all species have the same qrank value
dim(ans1$etamat) <- c(M,n) # bug: was c(n,M) prior to 20060822
list(deviance = ans1$deviance[1],
alldeviance = ans1$deviance[-1],
bcoefficients = ans1$bcoefficients,
bindex = ans1$bindex,
Bspline = Bspline2,
Cmat = matrix(cmatrix, p2, Rank, dimnames = list(
names(control$colx2.index), mynames5)),
coefficients = ans1$beta,
df1.nl = ans1$dofvec[1:NOS] - 1,
df2.nl = if (Rank == 2)
ans1$dofvec[2 * (1:NOS) - 1] - 1 else NULL,
lambda1 = ans1$lambda[1:NOS],
lambda2 = if (Rank == 2)
ans1$lambda[2 * (1:NOS) - 1] else NULL,
df.residual = n * M - qrank - sum(ans1$df - 1),
fitted = ans1$fv,
kindex = ans1$kindex,
predictors=matrix(ans1$etamat, n, M, byrow = TRUE),
wresiduals = ans1$zedd - t(ans1$etamat), # n x M
spar1 = ans1$smopar[1:NOS],
spar2 = if (Rank == 2) ans1$smopar[2 * (1:NOS) - 1] else NULL)
} else {
ans1$deriv
}
flush.console()
returnans
}
setClass(Class = "Coef.rrvgam", representation(
"Bspline" = "list",
"C" = "matrix",
"Constrained" = "logical",
"df1.nl" = "numeric",
"df2.nl" = "numeric",
"dispersion" = "numeric",
"eta2" = "matrix",
"latvar" = "matrix",
"latvar.order" = "matrix",
"M" = "numeric",
"Maximum" = "numeric",
"NOS" = "numeric",
"Optimum" = "matrix",
"Optimum.order"= "matrix",
"Rank" = "numeric",
"spar1" = "numeric",
"spar2" = "numeric"))
Coef.rrvgam <- function(object,
epsOptimum = 0.00001, # Determines how accurately Optimum is estimated
gridlen = 40, # Number of points on the grid (one level at a time)
maxgriditer = 10, # Maximum number of iters allowed for grid search
smallno = 0.05, ...) {
if (!is.Numeric(epsOptimum, positive = TRUE, length.arg = 1))
stop("bad input for argument 'epsOptimum'")
if (!is.Numeric(gridlen, positive = TRUE, integer.valued = TRUE) ||
gridlen < 5)
stop("bad input for argument 'gridlen'")
if (!is.Numeric(maxgriditer, positive = TRUE,
length.arg = 1, integer.valued = TRUE) ||
maxgriditer < 3)
stop("bad input for argument 'maxgriditer'")
if (!is.logical(ConstrainedO <- object@control$ConstrainedO))
stop("cannot determine whether the model is constrained or not")
if (!is.Numeric(smallno, positive = TRUE, length.arg = 1) ||
smallno > 0.5 || smallno < 0.0001)
stop("bad input for argument 'smallno'")
ocontrol <- object@control
if ((Rank <- ocontrol$Rank) > 2) stop("'Rank' must be 1 or 2")
gridlen <- rep_len(gridlen, Rank)
M <- if (any(slotNames(object) == "predictors") &&
is.matrix(object@predictors))
ncol(object@predictors) else
object@misc$M
NOS <- if (length(object@y)) ncol(object@y) else M
MSratio <- M / NOS
nice21 <- (length(ocontrol$colx1.index) == 1) &&
(names(ocontrol$colx1.index) == "(Intercept)")
if (!nice21)
stop("Can only handle 'noRRR = ~ 1'")
p1 <- length(ocontrol$colx1.index)
p2 <- length(ocontrol$colx2.index)
modelno <- object@control$modelno # 1,2,3,... or 0
ynames <- object@misc$ynames
if (!length(ynames)) ynames <- object@misc$predictors.names
if (!length(ynames)) ynames <- object@misc$ynames
if (!length(ynames)) ynames <- param.names("Y", NOS)
lp.names <- object@misc$predictors.names
if (!length(lp.names)) lp.names <- NULL
latvar.names <- param.names("latvar", Rank, skip1 = TRUE)
Cmat <- object@extra$Cmat # p2 x Rank (provided maxitl > 1)
if (ConstrainedO)
dimnames(Cmat) <- list(names(ocontrol$colx2.index), latvar.names)
latvar.mat <- if (ConstrainedO) {
object@x[, ocontrol$colx2.index, drop = FALSE] %*% Cmat
} else {
object@latvar
}
optimum <- matrix(NA_real_, Rank, NOS,
dimnames = list(latvar.names, ynames))
extents <- apply(latvar.mat, 2, range) # 2 by R
maximum <- rep_len(NA_real_, NOS)
which.species <- 1:NOS # Do it for all species
if (Rank == 1) {
gridd <- cbind(seq(extents[1, 1], extents[2, 1], len = gridlen))
eta2matrix <- matrix(0, NOS, 1) # Added 20160716
} else {
gridd <-
expand.grid(seq(extents[1, 1], extents[2, 1], len = gridlen[1]),
seq(extents[1, 2], extents[2, 2], len = gridlen[2]))
eta2matrix <- matrix(0, NOS, 1)
}
gridd.orig <- gridd
for (sppno in seq_along(which.species)) {
gridd <- gridd.orig
gridres1 <- gridd[2, 1] - gridd[1, 1]
gridres2 <- if (Rank == 2) gridd[2, 2] - gridd[1, 2] else 0
griditer <- 1
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], ynames) else which.species[sppno]
if (is.na(indexSpecies))
stop("mismatch found in 'which.species'")
while (griditer == 1 ||
((griditer <= maxgriditer) &&
((gridres1 > epsOptimum) ||
(gridres2 > epsOptimum)))) {
temp <- predictrrvgam(object, grid = gridd, sppno = thisSpecies,
Rank = Rank, deriv = 0, MSratio = MSratio)
yvals <- temp$yvals # gridlen-vector
xvals <- temp$xvals # gridlen x Rank; gridd
if (length(temp$eta2))
eta2matrix[sppno, 1] <- temp$eta2
nnn <- length(yvals)
index <- (1:nnn)[yvals == max(yvals)]
if (length(index) != 1)
warning("could not find a single maximum")
if (Rank == 2) {
initvalue <- rep_len(xvals[index,], Rank) # for optim()
if (abs(initvalue[1] - extents[1, 1]) < smallno)
initvalue[1] <- extents[1, 1] + smallno
if (abs(initvalue[1] - extents[2, 1]) < smallno)
initvalue[1] <- extents[2, 1] - smallno
if (abs(initvalue[2] - extents[1, 2]) < smallno)
initvalue[2] <- extents[1, 2] + smallno
if (abs(initvalue[2] - extents[2, 2]) < smallno)
initvalue[2] <- extents[2, 2] - smallno
break
}
if (index == 1 || index == nnn) {
maximum[sppno] <- optimum[1, sppno] <- NA
gridres1 <- epsOptimum + 1 # equivalent to a break
break # just in case
} else {
maximum[sppno] <- yvals[index] # On the eta scale
optimum[1, sppno] <- xvals[index, 1]
gridd[, 1] <- seq(
max(extents[1, 1], optimum[1, sppno] - gridres1),
min(extents[2, 1], optimum[1, sppno] + gridres1),
len = gridlen)
gridres1 <- gridd[2, 1] - gridd[1, 1]
griditer <- griditer + 1
}
} # of while
if (Rank == 2) {
myfun <- function(x, object, sppno, Rank = 1,
deriv = 0, MSratio = 1) {
x <- matrix(x, 1, length(x))
temp <- predictrrvgam(object, grid = x, sppno = sppno,
Rank = Rank, deriv = deriv, MSratio = MSratio)
temp$yval
}
answer <- optim(initvalue, myfun, gr = NULL, method = "L-BFGS-B",
lower = extents[1, ], upper = extents[2, ],
control = list(fnscale = -1), # maximize!
object = object, sppno = sppno, Rank = Rank,
deriv = 0, MSratio = MSratio)
for (rindex in 1:Rank)
if (abs(answer$par[rindex] - extents[1, rindex]) > smallno &&
abs(answer$par[rindex] - extents[2, rindex]) > smallno) {
optimum[rindex,sppno] <- answer$par[rindex]
maximum[sppno] <- answer$value
}
} # end of Rank = 2
} # end of sppno
myetamat <- rbind(maximum)
if (MSratio == 2)
myetamat <- kronecker(myetamat, matrix(1:0, 1, 2))
maximum <- object@family@linkinv(eta = myetamat,
extra = object@extra)
maximum <- c(maximum) # Convert from matrix to vector
names(maximum) <- ynames
ans <- new(Class = "Coef.rrvgam",
Bspline = object@Bspline,
Constrained = ConstrainedO,
df1.nl = object@extra$df1.nl,
latvar = latvar.mat,
latvar.order = latvar.mat,
Maximum = maximum,
M = M,
NOS = NOS,
Optimum = optimum,
Optimum.order = optimum,
Rank = Rank,
spar1 = object@extra$spar1)
if (ConstrainedO) {
ans@C <- Cmat
} else {
Cmat <- NULL
}
if (Rank == 2) {
dimnames(eta2matrix) <-
list(object@misc$predictors.names[c(FALSE, TRUE)], " ")
ans@eta2 <- eta2matrix
ans@df2.nl <- object@extra$df2.nl
ans@spar2 <- object@extra$spar2
}
for (rindex in 1:Rank) {
ans@Optimum.order[rindex, ] <- order(ans@Optimum[rindex, ])
ans@latvar.order[, rindex] <- order(ans@latvar[, rindex])
}
if (length(object@misc$estimated.dispersion) &&
object@misc$estimated.dispersion) {
p <- length(object@coefficients)
n <- object@misc$n
M <- object@misc$M
NOS <- if (length(object@y)) ncol(object@y) else M
pstar <- p + length(Cmat) # Adjustment
adjusted.dispersion <- object@misc$dispersion *
(n * M - p) / (n * M - pstar)
ans@dispersion <- adjusted.dispersion
}
if (MSratio == 2) {
lcoef <- object@coefficients
temp <- lcoef[((1:NOS)-1) * (2+Rank)+2]
names(temp) <- object@misc$predictors.names[2 * (1:NOS)]
ans@dispersion <- temp
}
dimnames(ans@Optimum) <- list(latvar.names, ynames)
ans
}
show.Coef.rrvgam <- function(object,
digits = max(2, options()$digits-2), ...) {
Rank <- object@Rank
NOS <- object@NOS
M <- object@M
Maximum <- if (length(object@Maximum))
cbind(Maximum = object@Maximum) else NULL
optmat <- cbind(t(object@Optimum))
dimnames(optmat) <- list(dimnames(optmat)[[1]],
if (Rank > 1)
paste("Optimum",
dimnames(optmat)[[2]], sep = ".") else
"Optimum")
if ( object@Constrained ) {
cat("\nC matrix (constrained/canonical coefficients)\n")
print(object@C, digits = digits, ...)
}
cat("\nOptimums and maximums\n")
print(cbind(Optimum = optmat,
Maximum), digits = max(1, digits-1))
cat("\nNonlinear degrees of freedom\n")
if (Rank == 1) {
print(cbind(df1.nl = object@df1.nl),
digits = max(2, digits-1), ...)
} else {
print(cbind(df1.nl = object@df1.nl,
df2.nl = object@df2.nl),
digits = max(2, digits-1), ...)
}
invisible(object)
}
setMethod("show", "Coef.rrvgam", function(object)
show.Coef.rrvgam(object))
setMethod("coef", "rrvgam", function(object, ...)
Coef.rrvgam(object, ...))
setMethod("coefficients", "rrvgam", function(object, ...)
Coef.rrvgam(object, ...))
setMethod("Coef", "rrvgam", function(object, ...)
Coef.rrvgam(object, ...))
lvplot.rrvgam <- function(object,
add = FALSE, show.plot = TRUE, rugplot = TRUE, y = FALSE,
type = c("fitted.values", "predictors"),
xlab = paste("Latent Variable",
if (Rank == 1) "" else " 1", sep = ""),
ylab = if (Rank == 1) switch(type, predictors = "Predictors",
fitted.values = "Fitted values") else "Latent Variable 2",
pcex = par()$cex, pcol = par()$col, pch = par()$pch,
llty = par()$lty, lcol = par()$col, llwd = par()$lwd,
label.arg= FALSE, adj.arg=-0.5,
sites= FALSE, spch = NULL, scol = par()$col,
scex = par()$cex, sfont = par()$font,
which.species = NULL,
check.ok = TRUE, ...) {
type <- match.arg(type, c("fitted.values", "predictors"))[1]
if ((Rank <- object@control$Rank) > 2)
stop("can only handle 'Rank' = 1 or 2 models")
M <- if (any(slotNames(object) == "predictors") &&
is.matrix(object@predictors))
ncol(object@predictors) else
object@misc$M
NOS <- ncol(object@y)
MSratio <- M / NOS
n <- object@misc$n
colx2.index <- object@control$colx2.index
cx1i <- object@control$colx1.index
if (!length(which.species))
which.species <- 1:NOS
if (check.ok)
if (!(length(cx1i) == 1 && names(cx1i) == "(Intercept)"))
stop("latent variable plots allowable only ",
"for 'noRRR = ~ 1' models")
Coeflist <- Coef(object)
Cmat <- Coeflist@C
latvarmat <- Coeflist@latvar # n x Rank
if (!show.plot)
return(latvarmat)
r.curves <- slot(object, type)
if (MSratio != 1 && type == "predictors")
stop("can only plot the predictors if M == S")
MorS <- ncol(r.curves) # Actually, here, the value is S always.
if (!add) {
if (Rank == 1) {
matplot(latvarmat,
if ( y && type == "fitted.values")
object@y[, which.species, drop = FALSE] else
r.curves[, which.species, drop = FALSE],
type = "n", xlab = xlab, ylab = ylab, ...)
} else { # Rank == 2
matplot(c(Coeflist@Optimum[1, which.species], latvarmat[, 1]),
c(Coeflist@Optimum[2, which.species], latvarmat[, 2]),
type = "n", xlab = xlab, ylab = ylab, ...)
}
}
pch <- rep_len(pch, length(which.species))
pcol <- rep_len(pcol, length(which.species))
pcex <- rep_len(pcex, length(which.species))
llty <- rep_len(llty, length(which.species))
lcol <- rep_len(lcol, length(which.species))
llwd <- rep_len(llwd, length(which.species))
adj.arg <- rep_len(adj.arg, length(which.species))
sppnames <- if (type == "predictors") dimnames(r.curves)[[2]] else
dimnames(object@y)[[2]]
if (Rank == 1) {
for (sppno in seq_along(which.species)) {
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], sppnames) else which.species[sppno]
if (is.na(indexSpecies))
stop("mismatch found in 'which.species'")
xx <- latvarmat
yy <- r.curves[, indexSpecies]
ooo <- sort.list(xx)
xx <- xx[ooo]
yy <- yy[ooo]
lines(xx, yy, col = lcol[sppno],
lwd = llwd[sppno], lty = llty[sppno])
if (y && type == "fitted.values") {
ypts <- object@y
if (NCOL(ypts) == ncol(r.curves))
points(xx, ypts[ooo, sppno], col = pcol[sppno],
cex = pcex[sppno], pch = pch[sppno])
}
}
if (rugplot) rug(xx)
} else {
if (sites) {
text(latvarmat[,1], latvarmat[,2], adj = 0.5,
labels = if (is.null(spch)) dimnames(latvarmat)[[1]] else
rep_len(spch, nrow(latvarmat)),
col = scol, cex = scex, font=sfont)
}
for (sppno in seq_along(which.species)) {
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], sppnames) else
which.species[sppno]
if (is.na(indexSpecies))
stop("mismatch found in 'which.species'")
points(Coeflist@Optimum[1, indexSpecies],
Coeflist@Optimum[2, indexSpecies],
col = pcol[sppno], cex = pcex[sppno], pch = pch[sppno])
}
if (label.arg) {
for (sppno in seq_along(which.species)) {
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], sppnames) else
which.species[sppno]
text(Coeflist@Optimum[1, indexSpecies],
Coeflist@Optimum[2, indexSpecies],
labels = (dimnames(Coeflist@Optimum)[[2]])[indexSpecies],
adj = adj.arg[sppno], col = pcol[sppno],
cex = pcex[sppno])
}
}
}
invisible(latvarmat)
}
setMethod("lvplot", "rrvgam",
function(object, ...) {
invisible(lvplot.rrvgam(object, ...))})
predict.rrvgam <- function (object, newdata = NULL,
type = c("link", "response", "terms"),
deriv = 0, ...) {
type <- match.arg(type, c("link", "response", "terms"))[1]
if (type != "link" && deriv != 0)
stop("Setting deriv = <positive integer> requires type='link'")
na.act <- object@na.action
object@na.action <- list()
ocontrol <- object@control
nice21 <- (length(ocontrol$colx1.index) == 1) &&
(names(ocontrol$colx1.index) == "(Intercept)")
if (!nice21)
stop("Can only handle 'noRRR = ~ 1'")
if (!length(newdata) && type == "response" &&
length(object@fitted.values)) {
if (length(na.act)) {
return(napredict(na.act[[1]], object@fitted.values))
} else {
return(object@fitted.values)
}
}
if (!length(newdata)) {
X <- model.matrixvlm(object, type = "lm", ...)
offset <- object@offset
tt <- terms(object)
if (!length(object@x))
attr(X, "assign") <- attrassignlm(X, tt)
} else {
if (is.smart(object) && length(object@smart.prediction)) {
setup.smart("read", smart.prediction = object@smart.prediction)
}
tt <- terms(object) # 20030811; object@terms$terms
X <- model.matrix(delete.response(tt), newdata,
contrasts = if (length(object@contrasts))
object@contrasts else NULL,
xlev = object@xlevels)
if (nice21 && nrow(X) != nrow(newdata)) {
as.save <- attr(X, "assign")
X <- X[rep_len(1, nrow(newdata)),, drop = FALSE]
dimnames(X) <- list(dimnames(newdata)[[1]], "(Intercept)")
attr(X, "assign") <- as.save # Restored
}
offset <- if (!is.null(off.num <- attr(tt, "offset"))) {
eval(attr(tt, "variables")[[off.num+1]], newdata)
} else if (!is.null(object@offset))
eval(object@call$offset, newdata)
if (is.smart(object) && length(object@smart.prediction)) {
wrapup.smart()
}
attr(X, "assign") <- attrassigndefault(X, tt)
}
cancoefs <- concoef(object)
latvarmat <- X[, ocontrol$colx2.index, drop = FALSE] %*% cancoefs
Rank <- ocontrol$Rank
NOS <- ncol(object@y)
sppnames <- dimnames(object@y)[[2]]
modelno <- ocontrol$modelno # 1,2,3,5 or 0
M <- if (any(slotNames(object) == "predictors") &&
is.matrix(object@predictors))
ncol(object@predictors) else
object@misc$M
MSratio <- M / NOS # 1st value is g(mean) = quadratic form in latvar
if (type == "terms") {
terms.mat <- matrix(0, nrow(X), Rank * NOS)
interceptvector <- rep_len(0, NOS)
} else {
etamat <- matrix(0, nrow(X), M) # Could contain derivatives
}
ind8 <- 1:Rank
which.species <- 1:NOS # Do it all for all species
for (sppno in seq_along(which.species)) {
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], sppnames) else which.species[sppno]
if (is.na(indexSpecies))
stop("mismatch found in 'which.species'")
temp345 <-
predictrrvgam(object, grid = latvarmat, sppno = thisSpecies,
Rank = Rank, deriv = deriv, MSratio = MSratio,
type = ifelse(type == "response", "link", type))
if (MSratio == 2) {
if (any(type == c("link", "response"))) {
etamat[, 2*sppno-1] <- temp345$yvals
etamat[, 2*sppno ] <- temp345$eta2
} else {
terms.mat[, ind8] <- temp345
interceptvector[sppno] <- attr(temp345, "constant")
}
} else {
if (any(type == c("link", "response"))) {
etamat[, sppno] <- temp345$yvals
} else {
terms.mat[, ind8] <- temp345
interceptvector[sppno] <- attr(temp345, "constant")
}
}
ind8 <- ind8 + Rank
}
if (length(offset) && any(offset != 0))
etamat <- etamat + offset
if (type == "link") {
dimnames(etamat) <-
list(dimnames(X)[[1]],
if (deriv == 0)
object@misc$predictors.names else NULL)
return(etamat)
} else if (type == "response") {
fv <- object@family@linkinv(etamat, extra = object@extra)
dimnames(fv) <- list(dimnames(fv)[[1]],
dimnames(object@fitted.values)[[2]])
return(fv)
} else {
attr(terms.mat, "constant") <- interceptvector
terms.mat
}
}
setMethod("predict", "rrvgam", function(object, ...)
predict.rrvgam(object, ...))
predictrrvgam <- function(object, grid, sppno, Rank = 1,
deriv = 0, MSratio = 1, type = "link") {
if (type != "link" && type != "terms")
stop("'link' must be \"link\" or \"terms\"")
if (ncol(grid <- as.matrix(grid)) != Rank)
stop("'grid' must have ", Rank, " columns")
if (!is.Numeric(1 + deriv, length.arg = 1,
positive = TRUE, integer.valued = TRUE))
stop("'deriv' must be a non-negative integer")
if (type == "terms" && deriv != 0)
stop("'deriv' must be 0 when type=\"terms\"")
temp.b <- object@Bspline[[sppno]]
if (type == "terms") {
meanlatvar <- colMeans(grid)
answer <- matrix(0, nrow(grid), Rank)
} else {
nlfunvalues <- 0
}
for (rindex in 1:Rank) {
temp <- temp.b[[rindex]] # temp is of class "vsmooth.spline.fit"
nlpart <- predict(temp, grid[, rindex], deriv = deriv)
yvals <- nlpart$y
if (type == "terms") {
answer[, rindex] <- yvals
} else {
nlfunvalues <- nlfunvalues + yvals
}
}
lcoef <- object@coefficients # linear coefs; dont use coef() (== Coef)
llcoef <- lcoef[(1+(sppno-1)*(MSratio+Rank)):(sppno*(MSratio+Rank))]
if (type == "terms") {
interceptvector <- llcoef[1]
for (rindex in 1:Rank) {
answer[, rindex] <- answer[, rindex] + (grid[, rindex] -
meanlatvar[rindex]) * llcoef[MSratio+rindex]
interceptvector <- interceptvector +
meanlatvar[rindex] * llcoef[MSratio+rindex]
}
} else {
linpar <- if (deriv == 0) {
llcoef[1] + grid %*% llcoef[-(1:MSratio)]
} else {
if (deriv == 1) llcoef[MSratio + rindex] else 0
}
nlfunvalues <- nlfunvalues + linpar # Now complete
}
if (type == "terms") {
attr(answer, "constant") <- interceptvector
answer
} else {
list(xvals = grid,
yvals = c(nlfunvalues),
eta2 = if (MSratio == 2) llcoef[MSratio] else NULL)
}
}
plot.rrvgam <-
function(x,
xlab = param.names("Latent Variable", Rank, skip1 = TRUE),
ylab = NULL, residuals.arg = FALSE,
pcol = par()$col, pcex = par()$cex, pch = par()$pch,
lcol = par()$col, lwd = par()$lwd, lty = par()$lty,
add = FALSE,
main = NULL,
center.cf = Rank > 1,
WhichRank = 1:Rank,
which.species = NULL, # a numeric or character vector
rugplot = TRUE, se.arg = FALSE, deriv = 0,
scale = 0, ylim = NULL,
overlay = FALSE, ...) {
Rank <- x@control$Rank
if (!is.logical(center.cf) || length(center.cf) != 1)
stop("bad input for argument 'center.cf'")
if (Rank > 1 && !center.cf)
stop("center.cf = TRUE is needed for models with Rank > 1")
NOS <- ncol(x@y)
sppnames <- dimnames(x@y)[[2]]
modelno <- x@control$modelno # 1,2,3, or 0
M <- if (any(slotNames(x) == "predictors") &&
is.matrix(x@predictors)) ncol(x@predictors) else x@misc$M
if (all((MSratio <- M / NOS) != c(1,2)))
stop("bad value for 'MSratio'")
pcol <- rep_len(pcol, Rank*NOS)
pcex <- rep_len(pcex, Rank*NOS)
pch <- rep_len(pch, Rank*NOS)
lcol <- rep_len(lcol, Rank*NOS)
lwd <- rep_len(lwd, Rank*NOS)
lty <- rep_len(lty, Rank*NOS)
xlab <- rep_len(xlab, Rank)
if (!length(which.species)) which.species <- 1:NOS
if (length(ylab))
ylab <- rep_len(ylab, length(which.species)) # Too long if overlay
if (length(main))
main <- rep_len(main, length(which.species)) # Too long if overlay
latvarmat <- latvar(x)
nice21 <- length(x@control$colx1.index) == 1 &&
names(x@control$colx1.index) == "(Intercept)"
if (!nice21)
stop("can only handle intercept-only models")
counter <- 0
for (sppno in seq_along(which.species)) {
thisSpecies <- which.species[sppno]
indexSpecies <- if (is.character(which.species))
match(which.species[sppno], sppnames) else which.species[sppno]
if (is.na(indexSpecies))
stop("mismatch found in 'which.species'")
terms.mat <- predictrrvgam(object = x, grid = latvarmat,
type = "terms",
sppno = indexSpecies, Rank = Rank,
deriv = deriv, MSratio = MSratio)
for (rindex in WhichRank) {
xvals <- latvarmat[, rindex]
yvals <- terms.mat[, rindex]
ooo <- sort.list(xvals)
xvals <- xvals[ooo]
yvals <- yvals[ooo]
if (!center.cf)
yvals <- yvals + attr(terms.mat, "constant")
if (!add)
if (sppno == 1 || !overlay) {
ylim.use <- if (length(ylim)) ylim else
ylim.scale(range(yvals), scale)
matplot(xvals, yvals, type = "n",
xlab = xlab[rindex],
ylab = if (length(ylab)) ylab[sppno] else
ifelse(overlay, "Fitted functions",
"Fitted function"),
main = if (length(main)) main[sppno] else
ifelse(overlay, "", sppnames[thisSpecies]),
ylim = ylim.use,
...)
}
if (residuals.arg) {
stop("cannot handle residuals = TRUE yet")
}
counter <- counter + 1
lines(xvals, yvals,
col = lcol[counter], lwd = lwd[counter], lty = lty[counter])
if (rugplot) rug(xvals)
}
}
invisible(x)
}
setMethod("plot", "rrvgam",
function(x, y, ...) {
if (!missing(y)) stop("cannot process the 'y' argument")
invisible(plot.rrvgam(x, ...))})
persp.rrvgam <-
function(x,
show.plot = TRUE, # zlim ignored if Rank == 1:
xlim = NULL, ylim = NULL, zlim = NULL,
gridlength = if (Rank == 1) 301 else c(51, 51),
which.species = NULL,
xlab = if (Rank == 1) "Latent Variable" else "Latent Variable 1",
ylab = if (Rank == 1) "Expected Value" else "Latent Variable 2",
zlab = "Expected value",
labelSpecies = FALSE, # For Rank == 1 only
stretch = 1.05, # quick and dirty, Rank == 1 only
main = "",
ticktype = "detailed",
col = if (Rank == 1) par()$col else "white",
lty = par()$lty,
lwd = par()$lwd,
rugplot = FALSE,
...) {
object <- x # don't like x as the primary argument
coefobj <- Coef(object)
if ((Rank <- coefobj@Rank) > 2)
stop("object must be a rank-1 or rank-2 model")
fvmat <- fitted(object)
NOS <- ncol(fvmat) # Number of species
M <- if (any(slotNames(object) == "predictors") &&
is.matrix(object@predictors)) ncol(object@predictors) else
object@misc$M
MSratio <- M / NOS # First value is g(mean) = quadratic form in latvar
xlim <- if (length(xlim))
xlim else
range(coefobj@latvar[, 1])
if (!length(ylim.orig <- ylim)) {
ylim <- if (Rank == 1)
c(0, max(fvmat)*stretch) else
range(coefobj@latvar[,2])
}
xlim <- rep_len(xlim, 2)
ylim <- rep_len(ylim, 2)
gridlength <- rep_len(gridlength, Rank)
latvar1 <- seq(xlim[1], xlim[2], length = gridlength[1])
latvar2 <- if (Rank == 2)
seq(ylim[1], ylim[2], len = gridlength[2]) else
NULL
latvarmat <- if (Rank == 2)
expand.grid(latvar1, latvar2) else
cbind(latvar1)
sppNames <- dimnames(object@y)[[2]]
if (!length(which.species)) {
which.species <- sppNames[1:NOS]
which.species.numer <- 1:NOS
} else
if (is.numeric(which.species)) {
which.species.numer <- which.species
which.species <- sppNames[which.species.numer] # Convert to character
} else {
which.species.numer <- match(which.species, sppNames)
}
LP <- matrix(NA_real_, nrow(latvarmat), NOS)
for (sppno in 1:NOS) {
temp <- predictrrvgam(object = object, grid = latvarmat, sppno = sppno,
Rank = Rank, deriv = 0, MSratio = MSratio)
LP[, sppno] <- temp$yval
}
if (MSratio == 2) {
LP <- kronecker(LP, matrix(1:0, 1, 2)) # n x M
}
fitvals <- object@family@linkinv(LP, extra = object@extra) # n by NOS
dimnames(fitvals) <- list(NULL, dimnames(fvmat)[[2]])
if (Rank == 1) {
if (show.plot) {
if (!length(ylim.orig)) # A revision:
ylim <- c(0, max(fitvals[,which.species.numer]) * stretch)
col <- rep_len(col, length(which.species.numer))
lty <- rep_len(lty, length(which.species.numer))
lwd <- rep_len(lwd, length(which.species.numer))
matplot(latvar1, fitvals, xlab = xlab, ylab = ylab,
type = "n", main = main, xlim = xlim, ylim = ylim, ...)
if (rugplot) rug(latvar(object))
for (sppno in seq_along(which.species.numer)) {
ptr2 <- which.species.numer[sppno] # points to species column
lines(latvar1, fitvals[,ptr2], col = col[sppno],
lty = lty[sppno], lwd = lwd [sppno], ...)
if (labelSpecies) {
ptr1 <- (1:nrow(fitvals))[max(fitvals[, ptr2]) ==
fitvals[, ptr2]]
ptr1 <- ptr1[1]
text(latvar1[ptr1], fitvals[ptr1, ptr2] + (stretch-1) *
diff(range(ylim)), label = sppNames[sppno],
col = col[sppno], ...)
}
}
}
} else {
max.fitted <- matrix(fitvals[,which.species[1]],
length(latvar1), length(latvar2))
if (length(which.species) > 1)
for (sppno in which.species[-1]) {
max.fitted <- pmax(max.fitted,
matrix(fitvals[, sppno],
length(latvar1), length(latvar2)))
}
if (!length(zlim))
zlim <- range(max.fitted, na.rm = TRUE)
perspdefault <- getS3method("persp", "default")
if (show.plot)
perspdefault(latvar1, latvar2, max.fitted,
zlim = zlim,
xlab = xlab, ylab = ylab, zlab = zlab,
ticktype = ticktype, col = col, main = main, ...)
}
invisible(list(fitted = fitvals,
latvar1grid = latvar1,
latvar2grid = if (Rank == 2) latvar2 else NULL,
max.fitted = if (Rank == 2) max.fitted else NULL))
}
if (!isGeneric("persp"))
setGeneric("persp", function(x, ...) standardGeneric("persp"))
setMethod("persp", "rrvgam", function(x, ...) persp.rrvgam(x = x, ...))
latvar.rrvgam <- function(object, ...) {
Coef(object, ...)@latvar
}
if (!isGeneric("lv"))
setGeneric("lv",
function(object, ...) {
.Deprecated("latvar")
standardGeneric("lv")
},
package = "VGAM")
setMethod("lv", "rrvgam",
function(object, ...) latvar.rrvgam(object, ...))
if (!isGeneric("latvar"))
setGeneric("latvar",
function(object, ...) standardGeneric("latvar"))
setMethod("latvar", "rrvgam",
function(object, ...) latvar.rrvgam(object, ...))
setClass(Class = "summary.rrvgam",
representation("misc" = "list",
"call" = "call"),
contains = "Coef.rrvgam")
summary.rrvgam <- function(object, ...) {
answer <- Coef(object, ...)
answer <- as(answer, "summary.rrvgam")
answer@misc <- object@misc
answer@call <- object@call
answer
}
setMethod("summary", "rrvgam", function(object, ...)
summary.rrvgam(object, ...))
show.summary.rrvgam <- function(x, ...) {
cat("\nCall:\n")
dput(x@call)
show.Coef.rrvgam(x, ...)
cat("\nNumber of species: ", x@NOS, "\n")
if (length(x@misc$dispersion) == 1) {
cat("\nDispersion parameter(s): ", x@misc$dispersion, "\n")
} else if (is.Numeric(x@dispersion)) {
cat("\nDispersion parameter(s)\n")
print( x@dispersion, ... )
}
invisible(x)
}
setMethod("show", "summary.rrvgam",
function(object)
show.summary.rrvgam(object))
concoef.rrvgam <- function(object, ...) {
Coef(object, ...)@C
}
concoef.Coef.rrvgam <- function(object, ...) {
if (length(list(...)))
warning("Too late! Ignoring the extra arguments")
object@C
}
if (FALSE) {
if (!isGeneric("ccoef"))
setGeneric("ccoef", function(object, ...) {
.Deprecated("concoef")
standardGeneric("ccoef")
})
setMethod("ccoef", "rrvgam", function(object, ...)
concoef.rrvgam(object, ...))
setMethod("ccoef", "Coef.rrvgam", function(object, ...)
concoef.Coef.rrvgam(object, ...))
}
setMethod("concoef", "rrvgam", function(object, ...)
concoef.rrvgam(object, ...))
setMethod("concoef", "Coef.rrvgam", function(object, ...)
concoef.Coef.rrvgam(object, ...))
if (!isGeneric("calibrate"))
setGeneric("calibrate", function(object, ...)
standardGeneric("calibrate"))
setMethod("calibrate", "rrvgam", function(object, ...)
calibrate.qrrvglm(object, ...))
setMethod("calibrate", "qrrvglm", function(object, ...)
calibrate.qrrvglm(object, ...))
setMethod("calibrate", "rrvglm", function(object, ...)
calibrate.rrvglm(object, ...))
Tol.rrvgam <- function(object, ...) {
stop("The tolerance for a 'rrvgam' object is undefined")
}
if (!isGeneric("Tol"))
setGeneric("Tol", function(object, ...) standardGeneric("Tol"))
setMethod("Tol", "rrvgam", function(object, ...)
Tol.rrvgam(object, ...))
setMethod("show", "rrvgam", function(object) show.vgam(object))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.