######
#cgam#
######
cgam <- function(formula, nsim = 0, family = gaussian, cpar = 1.2, data = NULL, weights = NULL, sc_x = FALSE, sc_y = FALSE, pnt = TRUE, pen = 0, var.est = NULL)
{
cl <- match.call()
if (is.character(family))
family <- get(family, mode = "function", envir = parent.frame())
if (is.function(family))
family <- family()
if (is.null(family$family))
stop("'family' not recognized!")
if (family$family == "ordered") {
rslt <- cgam.polr(formula, data = data, weights = weights, family = family, nsim = nsim, cpar = cpar)
rslt$call <- cl
return (rslt)
} else {
labels <- NULL
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
ynm <- names(mf)[1]
mt <- attr(mf, "terms")
y <- model.response(mf, "any")
if (family$family == "binomial") {
if (class(y) == "factor") {
y <- ifelse(y == levels(y)[1], 0, 1)
}
}
#additive
shapes1_add <- NULL; shapes2_add <- NULL; shapes_add <- NULL
xmat_add <- NULL; xnms_add <- NULL
tr <- NULL; pl <- NULL; umb <- NULL
tree.delta <- NULL; umbrella.delta <- NULL
tid1 <- NULL; tid2 <- NULL; tpos2 <- 0
uid1 <- NULL; uid2 <- NULL; upos2 <- 0
nums_add <- NULL; ks_add <- list(); sps_add <- NULL; xid_add <- 1
zmat <- NULL; zid <- NULL; zid0 <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_param <- NULL; is_fac <- NULL; vals <- NULL; st <- 1; ed <- 1
ztb <- list(); iztb <- 1
iadd <- 0
varlist_add <- NULL
#warp
warp.delta <- NULL
nums_wp <- NULL; ks_wp <- list(); ks_wps <- list(); sps_wp <- NULL
xmat_wp <- NULL; x1_wp <- NULL; x2_wp <- NULL; xnms_wp <- NULL; nks0_wp <- NULL; ks0_wp <- NULL; dc <- NULL
x1_wps <- NULL; x2_wps <- NULL; dcss <- list()
#zmat <- NULL; zid <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_fac <- NULL; is_param <- NULL; vals <- NULL; st <- 1; ed <- 1
iwps <- 0
varlist_wps <- NULL
#tri
tri.delta <- NULL
nums_tri <- NULL; ks_tri <- list(); ks_tris <- list(); nks_tris <- list(); sps_tri <- NULL
xmat_tri <- NULL; x1_tri <- NULL; x2_tri <- NULL; xnms_tri <- NULL; nks0_tri <- NULL; ks0_tri <- NULL; cvs <- NULL
x1_tris <- NULL; x2_tris <- NULL; cvss <- list()
#zmat <- NULL; zid <- NULL; zid1 <- NULL; zid2 <- NULL; znms <- NULL; is_fac <- NULL; is_param <- NULL; vals <- NULL; st <- 1; ed <- 1
#print (dim(mf))
itri <- 0
#print (head(mf))
for (i in 2:ncol(mf)) {
#additive
#print (attributes(mf[,i])$class)
if (!is.null(attributes(mf[,i])$categ)) {
if (is.numeric(attributes(mf[,i])$shape) & (attributes(mf[,i])$categ == "additive")) {
labels <- c(labels, "additive")
shapes1_add <- c(shapes1_add, attributes(mf[,i])$shape)
xmat_add <- cbind(xmat_add, mf[,i])
xnms_add <- c(xnms_add, attributes(mf[,i])$nm)
nums_add <- c(nums_add, attributes(mf[,i])$numknots)
sps_add <- c(sps_add, attributes(mf[,i])$space)
ks_add[[xid_add]] <- attributes(mf[,i])$knots
xid_add <- xid_add + 1
}
if (is.character(attributes(mf[,i])$shape) & (attributes(mf[,i])$categ == "additive")) {
labels <- c(labels, "additive")
shapes2_add <- c(shapes2_add, attributes(mf[,i])$shape)
if (attributes(mf[,i])$shape == "tree") {
pl <- c(pl, attributes(mf[,i])$pl)
treei <- tree.fun(mf[,i], attributes(mf[,i])$pl)
tree.delta <- rbind(tree.delta, treei)
tpos1 <- tpos2 + 1
tpos2 <- tpos2 + nrow(treei)
tid1 <- c(tid1, tpos1)
tid2 <- c(tid2, tpos2)
tr <- cbind(tr, mf[,i])
}
if (attributes(mf[,i])$shape == "umbrella") {
umbi <- umbrella.fun(mf[,i])
umbrella.delta <- rbind(umbrella.delta, umbi)
upos1 <- upos2 + 1
upos2 <- upos2 + nrow(umbi)
uid1 <- c(uid1, upos1)
uid2 <- c(uid2, upos2)
umb <- cbind(umb, mf[,i])
}
}
#warp
if (is.character(attributes(mf[, i])$shape) & (attributes(mf[,i])$categ == "warp")) {
iwps <- iwps + 1
labels <- c(labels, rep(paste("warp", iwps, sep = "_"), 2))
sps_wp <- attributes(mf[, i])$space
dcs <- attributes(mf[, i])$decreasing
dcss[[iwps]] <- dcs
ks0_wp <- attributes(mf[, i])$knots
nks0_wp <- attributes(mf[, i])$numknots
x1_wp <- (mf[, i])[, 1]
x1_wps <- cbind(x1_wps, x1_wp)
x2_wp <- (mf[, i])[, 2]
x2_wps <- cbind(x2_wps, x2_wp)
xmat_wp <- cbind(xmat_wp, mf[, i])
xnms_wp <- c(xnms_wp, attributes(mf[, i])$name)
#print (xnms_wp)
#print (dim(xmat_wp))
ans_warp <- makedelta_wps(x1t = x1_wp, x2t = x2_wp, m1_0 = nks0_wp[1], m2_0 = nks0_wp[2], k1 = ks0_wp$k1, k2 = ks0_wp$k2, space = sps_wp, decreasing = dcs)
warp.delta0 <- ans_warp$delta
k1 <- ans_warp$k1
k2 <- ans_warp$k2
ks_wp[[1]] <- k1
ks_wp[[2]] <- k2
ks_wps[[iwps]] <- ks_wp
if (iwps > 1) {
warp.delta0 <- warp.delta0[,-1]
}
varlist_wps <- c(varlist_wps, 1:ncol(warp.delta0)*0 + iwps)
warp.delta <- cbind(warp.delta, warp.delta0)
#print (dim(warp.delta))
#print (ks_wps)
}
#tri
if (is.character(attributes(mf[, i])$shape) & (attributes(mf[,i])$categ == "tri")) {
itri <- itri + 1
labels <- c(labels, rep(paste("tri", itri, sep = "_"), 2))
sps_tri <- attributes(mf[, i])$space
cvs <- attributes(mf[, i])$cvs
cvss[[itri]] <- cvs
ks0_tri <- attributes(mf[, i])$knots
nks0_tri <- attributes(mf[, i])$numknots
ks_tris[[itri]] <- ks0_tri
nks_tris[[itri]] <- nks0_tri
x1_tri <- (mf[, i])[, 1]
x1_tris <- cbind(x1_tris, x1_tri)
x2_tri <- (mf[, i])[, 2]
x2_tris <- cbind(x2_tris, x2_tri)
#xmat_tri <- cbind(x1_tri, x2_tri)
xmat_tri <- cbind(xmat_tri, mf[, i])
xnms_tri <- c(xnms_tri, attributes(mf[, i])$name)
}
#print (class((mf[, i])))
}
if (is.null(attributes(mf[,i])$shape)) {
if (!is.null(names(mf)[i])) {
znms <- c(znms, names(mf)[i])
}
if (!is.matrix(mf[,i])) {
zid <- c(zid, i)
is_param <- c(is_param, TRUE)
if (is.factor(mf[,i])) {
is_fac <- c(is_fac, TRUE)
ch_char <- suppressWarnings(is.na(as.numeric(levels(mf[, i]))))
if (any(ch_char)) {
vals <- c(vals, unique(levels(mf[, i]))[-1])
} else {
vals <- c(vals, as.numeric(levels(mf[, i]))[-1])
}
nlvs <- length(attributes(mf[,i])$levels)
ed <- st + nlvs - 2
zid1 <- c(zid1, st)
zid2 <- c(zid2, ed)
st <- st + nlvs - 1
zmat0 <- model.matrix(~ mf[, i])[, -1, drop = FALSE]
zmat <- cbind(zmat, zmat0)
ztb[[iztb]] <- mf[,i]
iztb <- iztb + 1
} else {
is_fac <- c(is_fac, FALSE)
zmat <- cbind(zmat, mf[, i])
ztb[[iztb]] <- mf[,i]
iztb <- iztb + 1
ed <- st
zid1 <- c(zid1, st)
zid2 <- c(zid2, ed)
st <- st + 1
vals <- c(vals, "")
}
} else {
is_param <- c(is_param, FALSE)
is_fac <- c(is_fac, FALSE)
zmat0 <- mf[, i]
mat_cols <- ncol(zmat0)
mat_rm <- NULL
#rm_num <- 0
for (irm in 1:mat_cols) {
if (all(round(diff(zmat0[, irm]), 8) == 0)) {
mat_rm <- c(mat_rm, irm)
}
}
if (!is.null(mat_rm)) {
zmat0 <- zmat0[, -mat_rm, drop = FALSE]
#rm_num <- rm_num + length(mat_rm)
}
zmat <- cbind(zmat, zmat0)
ztb[[iztb]] <- mf[,i]
iztb <- iztb + 1
vals <- c(vals, 1)
zid <- c(zid, i)
nlvs <- ncol(zmat0) + 1
ed <- st + nlvs - 2
zid1 <- c(zid1, st)
zid2 <- c(zid2, ed)
st <- st + nlvs - 1
}
}
}
dimnames(zmat)[[2]] <- NULL
if (family$family == "binomial" | family$family == "poisson" | family$family == "Gamma") {
wt.iter = TRUE
} else {wt.iter = FALSE}
#attr(xmat, "shape") <- shapes1
#print (labels)
if (any(labels == "additive") | !is.null(zmat)) {
xmat0_add <- xmat_add; shapes0_add <- shapes1_add; nums0_add <- nums_add; ks0_add <- ks_add; sps0_add <- sps_add; xnms0_add <- xnms_add; idx_s <- NULL; idx <- NULL
if (any(shapes1_add == 17)) {
kshapes <- length(shapes1_add)
obs <- 1:kshapes
idx_s <- obs[which(shapes1_add == 17)]; idx <- obs[which(shapes1_add != 17)]
xmat0_add[ ,1:length(idx_s)] <- xmat_add[ ,idx_s]
shapes0_add[1:length(idx_s)] <- shapes1_add[idx_s]
nums0_add[1:length(idx_s)] <- nums_add[idx_s]
sps0_add[1:length(idx_s)] <- sps_add[idx_s]
ks0_add[1:length(idx_s)] <- ks_add[idx_s]
xnms0_add[1:length(idx_s)] <- xnms_add[idx_s]
if (length(idx) > 0) {
xmat0_add[ ,(1 + length(idx_s)):kshapes] <- xmat_add[ ,idx]
shapes0_add[(1 + length(idx_s)):kshapes] <- shapes1_add[idx]
nums0_add[(1 + length(idx_s)):kshapes] <- nums_add[idx]
sps0_add[(1 + length(idx_s)):kshapes] <- sps_add[idx]
ks0_add[(1 + length(idx_s)):kshapes] <- ks_add[idx]
xnms0_add[(1 +length(idx_s)):kshapes] <- xnms_add[idx]
}
#xmat <- xmat0; nums <- nums0; ks <- ks0; sps <- sps0; xnms <- xnms0
}
shapes_add <- c(shapes1_add, shapes2_add)
}
#print (labels)
xnms <- c(xnms_wp, xnms_add, xnms_tri)
xmat <- cbind(xmat_wp, xmat_add, xmat_tri)
colnames(xmat) <- xnms
#boolz <- is.null(labels) & !is.null(zmat)
if (all(labels == "additive")) {
if (is.null(shapes_add)) {
nsim <- 0
}
ans <- cgam.fit(y = y, xmat = xmat0_add, zmat = zmat, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add, nsim = nsim, family = family, cpar = cpar, wt.iter = wt.iter, umbrella.delta = umbrella.delta, tree.delta = tree.delta, weights = weights, zid = zid, zid1 = zid1, zid2 = zid2, sc_x = sc_x, sc_y = sc_y, idx_s = idx_s, idx = idx, var.est = var.est, data = data)
if (!is.null(uid1) & !is.null(uid2)) {
uid1 <- uid1 + ans$d0 + ans$capm
uid2 <- uid2 + ans$d0 + ans$capm
}
if (!is.null(tid1) & !is.null(tid2)) {
tid1 <- tid1 + ans$d0 + ans$capm + ans$capu
tid2 <- tid2 + ans$d0 + ans$capm + ans$capu
}
#new:
knots_add <- ans$knots
numknots_add <- ans$numknots
#new:
if (length(knots_add) > 0) {
names(knots_add) <- xnms_add
}
rslt <- list(etahat = ans$etahat, muhat = ans$muhat, vcoefs = ans$vcoefs, xcoefs = ans$xcoefs, zcoefs = ans$zcoefs, ucoefs = ans$ucoefs, tcoefs = ans$tcoefs, coefs = ans$coefs, cic = ans$cic, d0 = ans$d0, edf0 = ans$edf0, etacomps = ans$etacomps, y = y, xmat_add = xmat_add, zmat = zmat, ztb = ztb, tr = tr, umb = umb, tree.delta = tree.delta, umbrella.delta = umbrella.delta, bigmat = ans$bigmat, shapes = shapes_add, shapesx = shapes1_add, shapesx0 = shapes0_add, prior.w = weights, wt = ans$wt, wt.iter = ans$wt.iter, family = family, SSE0 = ans$sse0, SSE1 = ans$sse1, pvals.beta = ans$pvals.beta, se.beta = ans$se.beta, null_df = ans$df.null, df = ans$df, resid_df_obs = ans$resid_df_obs, edf = ans$df_obs, null_deviance = ans$dev.null, deviance = ans$dev, tms = mt, capm = ans$capm, capms = ans$capms, capk = ans$capk, capt = ans$capt, capu = ans$capu, xid1 = ans$xid1, xid2 = ans$xid2, tid1 = tid1, tid2 = tid2, uid1 = uid1, uid2 = uid2, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, nsim = nsim, xnms_add = xnms_add, xnms = xnms, ynm = ynm, znms = znms, is_param = is_param, is_fac = is_fac, knots = knots_add, numknots = numknots_add, sps = sps_add, ms = ans$ms, cpar = ans$cpar, pl = pl, idx_s = idx_s, idx = idx, xmat0 = ans$xmat2, knots0 = ans$knots2, numknots0 = ans$numknots2, sps0 = ans$sps2, ms0 = ans$ms2, phihat = ans$phihat, pvs = ans$pvs, s.edf = ans$s.edf, bstats = ans$bstats, pvsz = ans$pvsz, z.edf = ans$z.edf, fstats = ans$fstats, vh = ans$vh, kts.var = ans$kts.var)
rslt$call <- cl
class(rslt) <- "cgam"
} else if (any(grepl("warp", labels, fixed = TRUE)) & !any(grepl("tri", labels, fixed = TRUE))) {
nprs <- sum(grepl("warp", labels, fixed = TRUE)) / 2
delta_add <- NULL
np_add <- 0
pb <- 0
if (any(labels == "additive")) {
#delta_add <- NULL
if (length(shapes0_add) > 0) {
ans_add <- make_delta_add(xmat = xmat0_add, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add)
delta_add <- ans_add$bigmat
#np_add: smooth only, conv and conc
np_add <- ans_add$np
pb <- nrow(delta_add)
varlist_add <- ans_add$varlist
}
}
delta_ut <- NULL
if (!is.null(umbrella.delta)) {
delta_add <- rbind(delta_add, umbrella.delta)
delta_ut <- rbind(delta_ut, umbrella.delta)
}
if (!is.null(tree.delta)) {
delta_add <- rbind(delta_add, tree.delta)
delta_ut <- rbind(delta_ut, tree.delta)
}
ans <- wps.fit(x1t = x1_wps, x2t = x2_wps, y = y, zmat = zmat, xmat_add = xmat_add, delta_add = delta_add, delta_ut = delta_ut, varlist_add = varlist_add, shapes_add = shapes_add, np_add = np_add, w = weights, pen = pen, pnt = pnt, cpar = cpar, decrs = dcss, delta = warp.delta, kts = ks_wps, wt.iter = wt.iter, family = family, nsim = nsim, nprs = nprs, idx_s = idx_s, idx = idx)
rslt <- list(ks_wps = ans$kts, muhat = ans$muhat, SSE1 = ans$sse1, SSE0 = ans$sse0, edfc = ans$edf, edf0 = ans$edf0, delta = ans$delta, y = y, zmat = ans$zmat, zmat_0 = ans$zmat_0, xmat_wp = xmat_wp, xmat_add = xmat_add, xmat = xmat, shapes = shapes_add, coefs = ans$coefs, zcoefs = ans$zcoefs, pvals.beta = ans$pz, se.beta = ans$sez, gcv = ans$gcv, xnms_wp = xnms_wp, xnms_add = xnms_add, xnms = xnms, xid_add = xid_add, znms = znms, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, ynm = ynm, decrs = dcss, tms = mt, is_param = is_param, is_fac = is_fac, d0 = ans$d0, pb = pb, pen = ans$pen, cpar = ans$cpar, wt.iter = wt.iter, cic = ans$cic, family = family, nprs_wps = nprs, varlist_wps = varlist_wps, coef_wp = ans$coef_wp, varlist_add = varlist_add, np_add = np_add, coef_add = ans$coef_add, coef_ut = ans$coef_ut, etacomps = ans$etacomps)
class(rslt) <- "wps"
if (!is.null(delta_add)) {
class(rslt) <- c("wps", "cgam")
if (min(which(labels == "additive")) < min(which(grepl("warp", labels, fixed = TRUE)))) {
class(rslt) <- c("cgam", "wps")
}
}
} else if (any(grepl("tri", labels, fixed = TRUE))) {
nprs <- sum(grepl("tri", labels, fixed = TRUE)) / 2
nprs_wps <- sum(grepl("warp", labels, fixed = TRUE)) / 2
amat_wp <- NULL
dmat_wp <- NULL
delta_add <- NULL
np_add <- 0
pb <- 0
if (any(labels == "additive")) {
if (length(shapes0_add) > 0) {
ans_add <- make_delta_add(xmat = xmat0_add, shapes = shapes0_add, numknots = nums0_add, knots = ks0_add, space = sps0_add)
delta_add <- ans_add$bigmat
np_add <- ans_add$np
pb <- nrow(delta_add)
varlist_add <- ans_add$varlist
}
}
if (!is.null(umbrella.delta)) {
delta_add <- rbind(delta_add, umbrella.delta)
}
if (!is.null(tree.delta)) {
delta_add <- rbind(delta_add, tree.delta)
}
if (any(grepl("warp", labels, fixed = TRUE))) {
ans_wps <- makeamat_wps(ks_wps, nprs_wps)
amat_wp <- ans_wps$amat
dmat_wp <- ans_wps$dmat
valist_wps <- ans_wps$valist_wps
}
ans <- trispl.fit(x1t = x1_tris, x2t = x2_tris, y = y, zmat = zmat, xmat_add = xmat_add, delta_add = delta_add, varlist_add = varlist_add, shapes_add = shapes_add, np_add = np_add, xmat_wp = xmat_wp, delta_wp = warp.delta, varlist_wps = varlist_wps, amat_wp = amat_wp, dmat_wp = dmat_wp, w = weights, lambda = pen, pnt = TRUE, cpar = cpar, cvss = cvss, delta = tri.delta, kts = ks_tris, nkts = nks_tris, wt.iter = wt.iter, family = family, nsim = nsim, nprs = nprs)
rslt <- list(muhat = ans$muhat, etahat = ans$etahat, trimat = ans$trimat, y = y, zmat = ans$zmat, capk = ans$capk, xmat_tri = xmat_tri, xmat_add = xmat_add, xmat_wp = xmat_wp, xmat = xmat, shapes = shapes_add, coefs = ans$thhat, zcoefs = ans$zcoefs, se.beta = ans$se.beta, pvals.beta = ans$pvals.beta, gcv = ans$gcv, edf = ans$edf, edf0 = ans$edf0, cic = ans$cic, sse0 = ans$sse0, sse1 = ans$sse1, family = family, nprs = nprs, nprs_wps = nprs_wps, varlist_wps = varlist_wps, varlist_add = varlist_add, varlist_tri = ans$varlist, xnms_tri = xnms_tri, xnms_add = xnms_add, xnms_wp = xnms_wp, xnms = xnms, znms = znms, zid = zid, vals = vals, zid1 = zid1, zid2 = zid2, ynm = ynm, decrs = dcss, cvss = cvss, tms = mt, is_param = is_param, is_fac = is_fac, d0 = ans$d0, pen = ans$pen, cpar = cpar, wt.iter = wt.iter, np_add = np_add, pb = pb, coef_add = ans$coef_add, coef_tri = ans$coef_tri, coef_wp = ans$coef_wp, etacomps = ans$etacomps, ks_wps = ks_wps, knots_lst = ans$knots_lst, trimat_lst = ans$trimat_lst, bmat_lst = ans$bmat_lst, capk_lst = ans$capk_lst)
class(rslt) <- "trispl"
if (!is.null(delta_add) & is.null(warp.delta)) {
if (labels[1] == "additive") {
class(rslt) <- c("cgam", "trispl")
} else {class(rslt) <- c("trispl", "cgam")}
}
if (!is.null(warp.delta) & is.null(delta_add)) {
#if (labels[1] == "warp") {
if (grepl("warp", labels[1], fixed = TRUE)) {
class(rslt) <- c("wps", "trispl")
} else {class(rslt) <- c("trispl", "wps")}
}
if (!is.null(warp.delta) & !is.null(delta_add)) {
if (labels[1] == "additive") {
class(rslt) <- c("cgam", "trispl", "wps")
} else if (grepl("warp", labels[1], fixed = TRUE)) {
class(rslt) <- c( "wps", "trispl", "cgam")
} else {class(rslt) <- c("trispl", "cgam", "wps")}
}
}
}
rslt$call <- cl
rslt$labels <- labels
return (rslt)
}
###############
#amat function#
###############
amat.fun <- function(x)
{
obs <- 1:length(x)
amat <- NULL
xu <- unique(x)
nx <- sort(xu)
hd <- head(nx, 1)
tl <- nx
while (length(tl) > 1) {
hd <- head(tl, 1)
tl <- tl[-1]
paired <- 0
for (i in 1:length(tl)) {
a1 <- 1:length(x)*0
if (hd * tl[i] > 0) {
if (hd < 0) {
if (tl[i] > hd) {
a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1
} else {
a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
}
}
if (hd > 0) {
if (tl[i] < hd) {
a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1
} else {
a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
}
}
#if (!all(a1 == 0) & paired == 0 ) {amat <- rbind(amat, a1); paired <- 1}
}
if (hd * tl[i] == 0) {
if (hd == 0) {
#if (tl[i] > 0) {
a1[min(obs[which(x == hd)])] <- 1; a1[min(obs[which(x == tl[i])])] <- -1
#}
} else {
a1[min(obs[which(x == hd)])] <- -1; a1[min(obs[which(x == tl[i])])] <- 1
}
}
if (!all(a1 == 0) & paired == 0 ) {amat <- rbind(amat, a1); paired <- 1}
}
}
dimnames(amat) <- NULL
amat
}
###############
#bmat function#
###############
bmat.fun <- function(x)
{
obs <- 1:length(x)
bmat <- NULL
hd <- head(x, 1)
tl <- x
j <- 0
while (length(tl) > 1) {
hd <- head(tl, 1)
tl <- tl[-1]
paired <- 0
j <- j + 1
for (i in 1:length(tl)) {
b1 <- 1:length(x)*0
if (hd == tl[i]) {b1[j] <- -1; b1[j + i] <- 1}
if (!all(b1 == 0) & paired == 0 ) {bmat <- rbind(bmat, b1); paired <- 1}
}
}
dimnames(bmat) <- NULL
bmat
}
##########
#cgam.fit#
##########
cgam.fit <- function(y, xmat, zmat, shapes, numknots, knots, space, nsim, family = gaussian(), cpar = 1.2, wt.iter = FALSE, umbrella.delta = NULL, tree.delta = NULL, weights = NULL, zid = zid, zid1 = zid1, zid2 = zid2, sc_x = FALSE, sc_y = FALSE, idx_s = NULL, idx = NULL, var.est = NULL, data = NULL) {
cicfamily <- CicFamily(family)
llh.fun <- cicfamily$llh.fun
#new: use log link in gamma
linkfun <- cicfamily$linkfun
etahat.fun <- cicfamily$etahat.fun
gr.fun <- cicfamily$gr.fun
wt.fun <- cicfamily$wt.fun
zvec.fun <- cicfamily$zvec.fun
muhat.fun <- cicfamily$muhat.fun
ysim.fun <- cicfamily$ysim.fun
deriv.fun <- cicfamily$deriv.fun
dev.fun <- cicfamily$dev.fun
n <- length(y)
sm <- 1e-7
#sm <- 1e-5
capl <- length(xmat) / n
if (capl < 1) {capl <- 0}
if (round(capl, 8) != round(capl, 1)) {stop ("Incompatible dimensions for xmat!")}
#new:
if (capl > 0 & sc_x) {
for (i in 1:capl) {xmat[,i] <- (xmat[,i] - min(xmat[,i])) / (max(xmat[,i]) - min(xmat[,i]))}
#for (i in 1:capl) {xmat[,i] <- (xmat[,i] - mean(xmat[,i])) / sd(xmat[,i])}
#for (i in 1:capl) {xmat[,i] <- xmat[,i] / sd(xmat[,i])}
}
#new:
if (sc_y) {
sc <- sd(y)
y <- y / sc
}
capk <- length(zmat) / n
if (capk < 1) {capk <- 0}
if (round(capk, 8) != round(capk, 1)) {stop ("Incompatible dimensions for zmat!")}
#new:
capls <- sum(shapes == 17)
####################################################
#get basis functions for the constrained components#
####################################################
delta <- NULL
varlist <- NULL
xid1 <- NULL; xid2 <- NULL; xpos2 <- 0
knotsuse <- list(); numknotsuse <- NULL
mslst <- list()
#new:
capm <- 0
capms <- 0
if (capl - capls > 0) {
del1_ans <- makedelta(xmat[, 1], shapes[1], numknots[1], knots[[1]], space = space[1])
del1 <- del1_ans$amat
knotsuse[[1]] <- del1_ans$knots
mslst[[1]] <- del1_ans$ms
numknotsuse <- c(numknotsuse, length(del1_ans$knots))
m1 <- length(del1) / n
#new code: record the number of columns of del1 if shapes0[1] == 17:
if (shapes[1] == 17) {capms <- capms + m1}
var1 <- 1:m1*0 + 1
xpos1 <- xpos2 + 1
xpos2 <- xpos2 + m1
xid1 <- c(xid1, xpos1)
xid2 <- c(xid2, xpos2)
if (capl == 1) {
delta <- del1
varlist <- var1
} else {
for (i in 2:capl) {
#new code:
del2_ans <- makedelta(xmat[,i], shapes[i], numknots[i], knots[[i]], space = space[i])
del2 <- del2_ans$amat
knotsuse[[i]] <- del2_ans$knots
mslst[[i]] <- del2_ans$ms
numknotsuse <- c(numknotsuse, length(del2_ans$knots))
m2 <- length(del2) / n
#new code: record the number of columns of del2 if shapes0[i] == 17:
if (shapes[i] == 17) {capms <- capms + m2}
xpos1 <- xpos2 + 1
xpos2 <- xpos2 + m2
xid1 <- c(xid1, xpos1)
xid2 <- c(xid2, xpos2)
delta <- rbind(del1, del2)
varlist <- 1:(m1 + m2)*0
varlist[1:m1] <- var1
varlist[(m1 + 1):(m1 + m2)] <- (1:m2)*0 + i
var1 <- varlist
m1 <- m1 + m2
del1 <- delta
}
}
if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0 & capk > 0) {
bigmat <- rbind(1:n*0 + 1, t(zmat), t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
np <- 1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) + capms
} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0 & capk == 0) {
bigmat <- rbind(1:n*0 + 1, t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
np <- 1 + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) + capms
} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0 & capk > 0) {
bigmat <- rbind(1:n*0 + 1, t(zmat), delta)
np <- 1 + capk + capms
} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0 & capk == 0) {
bigmat <- rbind(1:n*0 + 1, delta)
np <- 1 + capms
} else {
print ("error in capk, shapes!")
}
#new:
capm <- length(delta) / n - capms
} else {
if (capk + capls > 0) {
#new:
if (capls < 1 & capk > 0) {
bigmat <- rbind(1:n*0 + 1, t(zmat))
np <- 1 + capk
} else if (capls > 0) {
delta <- NULL; varlist <- NULL
del1_ans <- makedelta(xmat[,1], 17, numknots[1], knots[[1]], space = space[1])
del1 <- del1_ans$amat
knotsuse[[1]] <- del1_ans$knots
mslst[[1]] <- del1_ans$ms
numknotsuse <- c(numknotsuse, length(del1_ans$knots))
m1 <- length(del1) / n
var1 <- 1:m1*0 + 1
xpos1 <- xpos2 + 1
xpos2 <- xpos2 + m1
xid1 <- c(xid1, xpos1)
xid2 <- c(xid2, xpos2)
if (capls == 1) {
delta <- del1
varlist <- var1
} else {
for (i in 2:capls) {
del2_ans <- makedelta(xmat[,i], 17, numknots[i], knots[[i]], space = space[i])
del2 <- del2_ans$amat
knotsuse[[i]] <- del2_ans$knots
mslst[[i]] <- del2_ans$ms
numknotsuse <- c(numknotsuse, length(del2_ans$knots))
m2 <- length(del2) / n
xpos1 <- xpos2 + 1
xpos2 <- xpos2 + m2
xid1 <- c(xid1, xpos1)
xid2 <- c(xid2, xpos2)
delta <- rbind(del1, del2)
varlist <- 1:(m1 + m2)*0
varlist[1:m1] <- var1
varlist[(m1 + 1):(m1 + m2)] <- (1:m2)*0 + i
var1 <- varlist
m1 <- m1 + m2
del1 <- delta
}
}
if (capk < 1){
bigmat <- rbind(1:n*0 + 1, delta)
capms <- length(delta) / n
np <- 1 + capms
} else {
bigmat <- rbind(1:n*0 + 1, t(zmat), delta)
capms <- length(delta) / n
np <- 1 + capk + capms
}
}
} else {bigmat <- matrix(1:n*0 + 1, nrow = 1); capm <- 0; capms <- 0; np <- 1}
}
if (!is.null(umbrella.delta)) {
bigmat <- rbind(bigmat, umbrella.delta)
capu <- length(umbrella.delta) / n
} else {capu <- 0}
if (!is.null(tree.delta)) {
bigmat <- rbind(bigmat, tree.delta)
capt <- length(tree.delta) / n
} else {capt <- 0}
if (!is.null(umbrella.delta) | !is.null(tree.delta))
delta_ut <- rbind(umbrella.delta, tree.delta)
if (capl + capk + capu + capt > 0) {
# if (capl + capu + capt > 0) {
#new:
if (capl - capls + capu + capt > 0) {
#new:initialize cvec
cvec <- NULL
#new: initialize face
face <- NULL
if (wt.iter) {
etahat <- etahat.fun(n, y, fml = family$family)
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
} else {wt <- wt.fun(y, etahat, n, weights, fml = family$family)}
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
# gmat <- t(bigmat %*% sqrt(diag(wt)))
#new: avoid memory allocation error
gmat <- t(bigmat)
for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
zsend <- gmat[, 1:np, drop = FALSE]
#ans <- coneB(zvec, t(dsend), zsend)
ans <- coneB(zvec, dsend, zsend)
face <- ans$face
etahat <- t(bigmat) %*% ans$coefs
#new: for monotinic variance estimation
vh <- NULL
kts.var <- NULL
if (!is.null(var.est)) {
#x.var <- var.est
#test more:
if (!is.null(data)) {
nms <- names(data)
nm <- attributes(var.est)$nm
if (nm %in% nms) {
x.var <- data[[nm]]
attributes(x.var) <- attributes(var.est)
}
} else {x.var <- var.est}
db.exp <- attributes(x.var)$db.exp
kts.var <- attributes(x.var)$var.knots
shape.var <- attributes(x.var)$shape
iter <- 0
diff.var <- 100
oldeta <- etahat
evec0 <- y - etahat
#bigwt <- 10*max(evec0)
bigwt <- 1000
while(diff.var > 1e-8 & iter < 10) {
iter <- iter + 1
evec <- y - etahat
fit.var <- varest(evec, x.var, shape=shape.var, var.knots=kts.var, db.exp=db.exp)
v1 <- fit.var$vhat
wt0 <- 1/v1
wt0[wt0 > bigwt] <- bigwt
wt <- wt.fun(y, etahat, n, weights = wt0, fml = family$family)
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
gmat <- t(bigmat)
for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
zsend <- gmat[, 1:np, drop = FALSE]
ans <- coneB(zvec, dsend, zsend)
etahat <- t(bigmat) %*% ans$coefs
diff.var <- mean((etahat - oldeta)^2)
oldeta <- etahat
}
kts.var <- fit.var$var.knots
vh <- v1
}
if (wt.iter) {
muhat <- muhat.fun(etahat, fml = family$family)
diff <- 1
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
nrep <- 0
##########
#iterate!#
##########
while (diff > sm & mdiff & nrep < n^2){
oldmu <- muhat
nrep <- nrep + 1
gr <- gr.fun(y, etahat, weights, fml = family$family)
wt <- wt.fun(y, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
#zvec <- cvec / sqrt(wt)
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
# gmat <- t(bigmat %*% sqrt(diag(wt)))
gmat <- t(bigmat)
for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
zsend <- gmat[, 1:np, drop = FALSE]
#ans <- coneB(zvec, t(dsend), zsend)
ans <- coneB(zvec, dsend, zsend, face = face)
etahat <- t(bigmat) %*% ans$coefs
muhat <- muhat.fun(etahat, fml = family$family)
diff <- mean((muhat - oldmu)^2)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
}
}
yhat <- ans$yhat
coefskeep <- ans$coefs
########################
#if capk >= 0, we have:#
########################
zcoefs <- coefskeep[1:(capk + 1)]
######################
#we will always have:#
######################
vcoefs <- coefskeep[1:np]
#######################
#if capm > 0, we have:#
#######################
xcoefs <- NULL
#if (capm > 0) {
# xcoefs <- coefskeep[(np + 1):(np + capm)]
#}
#new:
if (capl > 0) {
xcoefs <- coefskeep[(np - capms + 1):(np + capm)]
}
#######################
#if capu > 0, we have:#
#######################
ucoefs <- NULL
if (capu > 0) {
ucoefs <- coefskeep[(np + 1 + capm):(np + capm + capu)]
}
#######################
#if capt > 0, we have:#
#######################
tcoefs <- NULL
if (capt > 0) {
tcoefs <- coefskeep[(np + 1 + capm + capu):(np + capm + capu + capt)]
}
#########################################################
#if we have at least one constrained predictor, we have:#
#########################################################
thvecs <- NULL
if (capl > 0) {
#new code:
dcoefs <- coefskeep[(np - capms + 1):(np + capm)]
#dcoefs <- coefskeep[(np + 1):(np + capm)]
#####################################################
#thvecs is f(x), where x has one of the eight shapes#
#####################################################
thvecs <- matrix(nrow = capl, ncol = n)
ncon <- 1
for (i in 1:capl) {
thvecs[i,] <- t(delta[varlist == i,]) %*% dcoefs[varlist == i]
#new:
#if (shapes[i] > 2 & shapes[i] < 5) {
if (shapes[i] > 2 & shapes[i] < 5 | shapes[i] > 10 & shapes[i] < 13) {
ncon <- ncon + 1
#thvecs[i,] <- thvecs[i,] + zcoefs[ncon] * xmat[,i]
thvecs[i,] <- thvecs[i,] + vcoefs[capk + ncon] * xmat[,i]
}
}
}
#new:order thvecs back
#print (idx_s)
#print (idx)
#if (!is.null(idx_s)) {
if (length(idx_s) > 0) {
thvecs0 <- thvecs
thvecs0[idx_s,] <- thvecs[1:length(idx_s), ]
#if (!is.null(idx)) {
if (length(idx) > 0) {
thvecs0[idx,] <- thvecs[(1+length(idx_s)):capl, ]
}
thvecs <- thvecs0
}
thvecs_ut <- NULL
if (capu + capt > 0) {
thvecs_ut <- t(delta_ut) %*% coefskeep[(np + 1 + capm):(np + capm + capu + capt)]
}
if (!is.null(thvecs_ut)) {
thvecs <- rbind(thvecs, t(thvecs_ut))
}
#new: problem when not gaussian
if (sc_y) {
y <- y*sc
etahat <- etahat*sc
for (i in 1:nrow(thvecs)) {
thvecs[i,] <- thvecs[i,] * sc
}
}
etakeep <- etahat
muhatkeep <- muhat.fun(etakeep, fml = family$family)
wtkeep <- wt
df_obs <- sum(abs(coefskeep) > 0)
#llh <- llh.fun(y, muhat, etahat, n, weights, fml = family$family)
if (family$family == "Gamma") {
if ((n - cpar * df_obs) <= 0) {
phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / df_obs
} else {
phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / (n - cpar * df_obs)
}
} else {phihat <- NULL}
#print (phihat)
llh <- llh.fun(y, muhatkeep, etakeep, phihat, n, weights, fml = family$family)
if (family$family == "poisson") {
mu0 <- mean(y)
eta0 <- log(mu0)
} else {mu0 <- NULL}
#new: get p-values for smooth components
pvs <- NULL
s.edf <- NULL
bstats <- NULL
if (is.numeric(shapes)) {
#changed to include shape==17
if (all(shapes >= 9 & shapes <= 17)) {
for (i in 1:capl) {
ansi <- cgam.pv(y=y, xmat=xmat, zmat=zmat, shapes=shapes, delta=delta, np=np, capms=capms, numknotsuse=numknotsuse, varlist=varlist, family=family, weights=weights, test_id=i, nsims=100)
pvi <- ansi$pv
edfi <- 1.5*sum(xcoefs[varlist == i] > 1e-8)
bstati <- ansi$bstat
pvs <- c(pvs, pvi)
s.edf <- c(s.edf, edfi)
bstats <- c(bstats, bstati)
}
}
}
#new: get p-values for categorical predictors, not for each level
pvsz <- NULL
z.edf <- NULL
fstats <- NULL
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
prior.w <- weights
if (capk > 0) {
sse1 <- sum(prior.w * (y - muhatkeep)^2)
ansi <- cgam.pvz(y=y, bigmat=bigmat, df_obs=df_obs, sse1=sse1, np=np, zid=zid, zid1=zid1, zid2=zid2, muhat = muhatkeep, etahat = etakeep, coefskeep = coefskeep, wt.iter=wt.iter, family=family, weights=weights)
pvsz <- ansi$pvs
z.edf <- ansi$edfs
fstats <- ansi$fstats
}
} else if (capk + capls > 0 & capl - capls + capt + capu == 0) {
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
prior.w <- weights
vmat <- t(bigmat[1:np, , drop = FALSE])
if (wt.iter) {
nrep <- 0
muhat <- mean(y) + 1:n*0
etahat <- linkfun(muhat)
diff <- 1
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
while (diff > sm & mdiff & nrep < n^2) {
nrep <- nrep + 1
oldmu <- muhat
zhat <- etahat + (y - muhat) * deriv.fun(muhat, fml = family$family)
#w <- diag(as.vector(prior.w / deriv.fun(muhat)))
#w <- diag(as.vector(prior.w * (deriv.fun(muhat, fml = family$family))^(-1)))
#b <- solve(t(vmat) %*% w %*% vmat) %*% t(vmat) %*% w %*% zhat
w <- as.vector(prior.w * (deriv.fun(muhat, fml = family$family))^(-1))
tvmat <- t(vmat)
for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
#print (tvmat)
b <- solve(tvmat %*% vmat) %*% tvmat %*% zhat
etahat <- vmat %*% b
muhat <- muhat.fun(etahat, fml = family$family)
diff <- mean((muhat - oldmu)^2)
mdiff <- abs(max(muhat) - 1)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
}
zcoefs <- b[1:(capk + 1)]
#se.beta <- sqrt(diag(solve(t(vmat) %*% w %*% vmat)))[1:(capk + 1)]
se.beta <- sqrt(diag(solve(tvmat %*% vmat)))[1:(capk + 1)]
zstat <- zcoefs / se.beta
pvals.beta <- 1 - pchisq(zstat^2, df = 1)
#new: gamma only
if (family$family == "Gamma") {
phihat <- sum(((y - muhatkeep) / muhatkeep)^2) / (n - np)
} else {phihat <- NULL}
} else {
#w <- diag(prior.w)
#b <- solve(t(vmat) %*% w %*% vmat) %*% t(vmat) %*% w %*% y
w <- prior.w
tvmat <- t(vmat)
for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
b <- solve(tvmat %*% vmat) %*% tvmat %*% y
etahat <- vmat %*% b
muhat <- muhat.fun(etahat, fml = family$family)
sdhat2 <- sum(prior.w * (y - muhat)^2) / (n - np)
zcoefs <- b[1:(capk + 1)]
se.beta <- sqrt(diag(solve(tvmat %*% vmat) * sdhat2))[1:(capk + 1)]
tstat <- zcoefs / se.beta
pvals.beta <- (1 - pt(abs(tstat), df = n - np)) * 2
#new: for var estimation
vh <- NULL
kts.var <- NULL
if (!is.null(var.est)) {
x.var <- var.est
db.exp <- attributes(x.var)$db.exp
kts.var <- attributes(x.var)$var.knots
shape.var <- attributes(x.var)$shape
iter <- 0
diff.var <- 100
oldeta <- etahat
evec0 <- y - oldeta
#bigwt <- 10*max(evec0)
bigwt <- 1000
while(diff.var > 1e-8 & iter < 10) {
iter <- iter + 1
evec <- y - etahat
fit.var <- varest(evec, x.var, shape=shape.var, var.knots=kts.var, db.exp=db.exp)
v1 <- fit.var$vhat
#w <- 1/v1
w0 <- 1/v1
w0[w0 > bigwt] <- bigwt
w <- w0
tvmat <- t(vmat)
for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
b <- solve(tvmat %*% vmat) %*% tvmat %*% y
etahat <- vmat %*% b
diff.var <- mean((etahat - oldeta)^2)
oldeta <- etahat
#print (diff.var)
}
kts.var <- fit.var$var.knots
vh <- v1
}
}
#add thvecs if capls > 0:
thvecs <- NULL
if (capls > 0) {
thvecs <- matrix(nrow = capls, ncol = n)
dcoefs <- b[(capk + 2):np]
for (i in 1:capls) {
thvecs[i,] <- t(delta[varlist == i,]) %*% dcoefs[varlist == i]
}
}
#new:
if (sc_y) {
y <- y*sc
etahat <- etahat*sc
muhat <- muhat.fun(etahat, fml = family$family)
for (i in 1:nrow(thvecs)) {
thvecs[i,] = thvecs[i,] * sc
}
}
llh <- llh.fun(y, muhat, etahat, phihat, n, weights, fml = family$family)
df_obs <- np
dfmean <- np
rslt <- new.env()
rslt$family <- family
rslt$wt.iter <- wt.iter
#rslt$wt <- diag(w)
rslt$wt <- w
rslt$bigmat <- bigmat
rslt$etahat <- etahat
rslt$muhat <- muhat
rslt$d0 <- np
rslt$capm <- 0
rslt$capms <- capms
rslt$capk <- capk
rslt$capu <- capu
rslt$capt <- capt
rslt$xid1 <- xid1 + np - capms
rslt$xid2 <- xid2 + np - capms
rslt$dfmean <- dfmean
rslt$edf0 <- dfmean
#rslt$llh <- llh
#if (nsim > 0) {
if (family$family == "binomial") {
nobs <- sum(prior.w)
} else {nobs <- n}
if ((nobs - np - cpar * (dfmean - np)) <= 0) {
rslt$cic <- llh + log(1 + 2 * dfmean / (dfmean - np))
} else {
#new:
if (n<=200){cpar=1.5}
rslt$cic <- llh + log(1 + 2 * dfmean / (nobs - np - cpar * (dfmean - np)))
}
#rslt$cic <- llh + log(1 + 2 * dfmean / (n - np - 1.5 * (dfmean - np)))
#}
rslt$zcoefs <- zcoefs
rslt$coefs <- b
rslt$vcoefs <- b
rslt$xcoefs <- b[(capk + 2):np]
rslt$se.beta <- se.beta
rslt$pvals.beta <- pvals.beta
rslt$dev <- dev.fun(y, muhat, etahat, weights, fml = family$family)$dev
rslt$dev.null <- dev.fun(y, muhat, etahat, weights, fml = family$family)$dev.null
rslt$df <- n - np
rslt$df.null <- n - 1
#rslt$resid_df_obs <- n - np - 1.5 * (df_obs - np)
rslt$resid_df_obs <- n - cpar * df_obs
rslt$vhat <- etahat
rslt$vmat <- vmat
rslt$etacomps <- thvecs
rslt$knots <- knotsuse
rslt$numknots <- numknotsuse
rslt$ms <- mslst
#new: used for predict when there are only shape == 17 x's
#rslt$xmat <- xmat
rslt$xmat2 <- xmat
rslt$knots2 <- knotsuse
rslt$numknots2 <- numknotsuse
rslt$ms2 <- mslst
rslt$vh <- vh
rslt$kts.var <- kts.var
return (rslt)
}
##########
#get cic#
##########
if (capl - capls + capu + capt > 0 & nsim > 0) {
dfs <- 1:nsim
#new: initialize face
face <- NULL
for (isim in 1:nsim) {
#set.seed(123)
ysim <- ysim.fun(n, mu0, fml = family$family)
if (wt.iter) {
etahat <- etahat.fun(n, ysim, fml = family$family)
gr <- gr.fun(ysim, etahat, weights, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
} else {wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)}
zvec <- zvec.fun(cvec, wt, ysim, fml = family$family)
# gmat <- t(bigmat %*% sqrt(diag(wt)))
gmat <- t(bigmat)
for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
zsend <- gmat[, 1:np, drop = FALSE]
#ans <- try(coneB(zvec, t(dsend), zsend))
ans <- try(coneB(zvec, dsend, zsend))
face <- ans$face
if (class(ans) == "try-error") next
if (wt.iter) {
etahat <- t(bigmat) %*% ans$coefs
muhat <- muhat.fun(etahat, fml = family$family)
diff <- 1
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
##########
#iterate!#
##########
nrep <- 0
while (diff > sm & nrep < n^2 & mdiff > sm) {
nrep <- nrep + 1
oldmu <- muhat
gr <- gr.fun(ysim, etahat, weights, fml = family$family)
wt <- wt.fun(ysim, etahat, n, weights, fml = family$family)
cvec <- wt * etahat - gr
#zvec <- cvec / sqrt(wt)
zvec <- zvec.fun(cvec, wt, y, fml = family$family)
# gmat <- t(bigmat %*% sqrt(diag(wt)))
gmat <- t(bigmat)
for (i in 1:n) {gmat[i,] <- bigmat[,i] * sqrt(wt[i])}
dsend <- gmat[, (np + 1):(np + capm + capu + capt), drop = FALSE]
zsend <- gmat[, 1:np, drop = FALSE]
#ans <- try(coneB(zvec, t(dsend), zsend))
ans <- try(coneB(zvec, dsend, zsend, face=face))
if (class(ans) == "try-error") next
etahat <- t(bigmat) %*% ans$coefs
muhat <- muhat.fun(etahat, fml = family$family)
diff <- mean((muhat - oldmu)^2)
if (family$family == "binomial") {
mdiff <- abs(max(muhat) - 1) > sm
} else {mdiff <- TRUE}
}
}
dfs[isim] <- sum(abs(ans$coefs) > 0)
}
dfmean <- mean(dfs)
#print (dfmean)
} else if (capl - capls + capu + capt > 0 & nsim == 0) {
dfmean <- NULL
}
###################################################
#if the user does not give any predictor, we have:#
###################################################
} else {
rslt <- new.env()
rslt$family <- family
rslt$wt.iter <- wt.iter
rslt$muhat <- 1:n*0 + mean(y)
rslt$etahat <- linkfun(muhat)
rslt$dfmean <- 1
print ("No predictor is provided")
return (rslt)
}
if (capl - capls + capu + capt > 0) {
#new:
#exclude the case we only have z or unrestricted smooth
xid1 <- xid1 + np - capms
xid2 <- xid2 + np - capms
#xid1 <- xid1 + np
#xid2 <- xid2 + np
rslt <- new.env()
rslt$phihat <- phihat
rslt$family <- family
rslt$wt.iter <- wt.iter
rslt$wt <- wtkeep
rslt$bigmat <- bigmat
rslt$etahat <- etakeep
rslt$muhat <- muhatkeep
rslt$d0 <- np
rslt$capm <- capm
rslt$capms <- capms
rslt$capk <- capk
rslt$capu <- capu
rslt$capt <- capt
rslt$xid1 <- xid1
rslt$xid2 <- xid2
rslt$coefs <- coefskeep
rslt$vcoefs <- vcoefs
rslt$xcoefs <- xcoefs
rslt$zcoefs <- zcoefs
rslt$ucoefs <- ucoefs
rslt$tcoefs <- tcoefs
rslt$dfmean <- dfmean
#rslt$llh <- llh
#if (nsim > 0) {
if (!is.null(dfmean)) {
#check!
#new: use dfmean - np if n - np - 1.5 * (dfmean - np) < 0
if ((n - np - cpar * (dfmean - np)) <= 0) {
rslt$cic <- llh + log(1 + 2 * dfmean / (dfmean - np))
} else {
#new:
if (n<=200){cpar=1.5}
rslt$cic <- llh + log(1 + 2 * dfmean / (n - np - cpar * (dfmean - np)))
}
} else {rslt$cic <- NULL}
rslt$edf0 <- dfmean
rslt$etacomps <- thvecs
vmat <- t(bigmat[1:np, , drop = FALSE])
rslt$vmat <- vmat
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
prior.w <- weights
#w <- diag(as.vector(prior.w / deriv.fun(muhatkeep, fml = family$family)))
w <- as.vector(prior.w / deriv.fun(muhatkeep, fml = family$family))
###############################################
#the case capk = 0 and capk >= 0 are combined:#
###############################################
#debugged: vhat -> muhat.fun(vhat)
vhat <- vmat %*% vcoefs
muvhat <- muhat.fun(vhat, fml = family$family)
#debugged: yhat -> muhatkeep
sse1 <- sum(prior.w * (y - muhatkeep)^2)
sse0 <- sum(prior.w * (y - muvhat)^2)
#new:
#if ((n - np - cpar * df_obs) <= 0) {
if ((n - cpar * df_obs) <= 0) {
#sdhat2 <- sse1 / (df_obs - np)
sdhat2 <- sse1 / df_obs
} else {
#sdhat2 <- sse1 / (n - np - 1.5 * (df_obs - np))
#sdhat2 <- sse1 / (n - np - cpar * df_obs)
sdhat2 <- sse1 / (n - cpar * df_obs)
}
#debugged: vmat -> vmat and duse
#new: coefskeep include zcoefs; bigmat include vmat
pmat <- vmat
capbm <- length(bigmat) / n
bigmat_nv <- bigmat[(np + 1):capbm, , drop = FALSE]
coefs_nv <- coefskeep[(np + 1):capbm]
duse <- coefs_nv > 1e-8
if (sum(duse) >= 1) {
pmat <- cbind(vmat, t(bigmat_nv[duse, , drop = FALSE]))
}
if (wt.iter) {
#new:
tpmat <- t(pmat)
for (i in 1:n) {tpmat[,i] <- tpmat[,i] * w[i]}
se2 <- solve(tpmat %*% pmat)
} else {
# se2 <- solve(t(pmat) %*% diag(prior.w) %*% pmat) * sdhat2
tpmat <- t(pmat)
for (i in 1:n) {tpmat[,i] <- tpmat[,i] * prior.w[i]}
se2 <- solve(tpmat %*% pmat) * sdhat2
}
se.beta <- 1:(capk + 1)*0
tstat <- 1:(capk + 1)*0
pvals.beta <- 1:(capk + 1)*0
rslt$zcoefs <- zcoefs
for (i in 1:(capk + 1)) {
se.beta[i] <- sqrt(se2[i,i])
tstat[i] <- zcoefs[i] / se.beta[i]
#new code: n - np - 1.5 * (df_obs - np) must be positive
#if ((n - np - cpar * df_obs) <= 0) {
if ((n - cpar * df_obs) <= 0) {
pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]), df_obs))
warning ('Effective degrees of freedom is close to the number of observations! Inference about parametric covariates is not reliable!')
} else {
#pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]), n - np - cpar * df_obs))
pvals.beta[i] <- 2 * (1 - pt(abs(tstat[i]), n - cpar * df_obs))
}
}
rslt$se.beta <- se.beta
rslt$pvals.beta <- pvals.beta
rslt$sse1 <- sse1
rslt$sse0 <- sse0
rslt$dev <- dev.fun(y, muhatkeep, etakeep, weights, fml = family$family)$dev
rslt$dev.null <- dev.fun(y, muhatkeep, etakeep, weights, fml = family$family)$dev.null
rslt$df <- n - np
rslt$df.null <- n - 1
#rslt$resid_df_obs <- n - np - cpar * df_obs
rslt$resid_df_obs <- n - cpar * df_obs
rslt$df_obs <- df_obs
rslt$vhat <- vhat
if (length(knotsuse) == 0) {
knotsuse <- NULL
}
#new: order back
knotsuse2 <- knotsuse
numknotsuse2 <- numknotsuse
mslst2 <- mslst
xmat2 <- xmat
#if (!is.null(idx_s)) {
if (length(idx_s) > 0) {
knotsuse0 <- knotsuse
numknotsuse0 <- numknotsuse
mslst0 <- mslst
knotsuse0[idx_s] <- knotsuse[1:length(idx_s)]
numknotsuse0[idx_s] <- numknotsuse[1:length(idx_s)]
mslst0[idx_s] <- mslst[1:length(idx_s)]
#if (!is.null(idx)) {
if (length(idx) > 0) {
knotsuse0[idx] <- knotsuse[(1+length(idx_s)):capl]
numknotsuse0[idx] <- numknotsuse[(1+length(idx_s)):capl]
mslst0[idx] <- mslst[(1+length(idx_s)):capl]
}
knotsuse <- knotsuse0
numknotsuse <- numknotsuse0
mslst <- mslst0
}
# the following has shape = 17 at the beginning: knots2 etc
rslt$knots2 <- knotsuse2
rslt$numknots2 <- numknotsuse2
rslt$ms2 <- mslst2
rslt$xmat2 <- xmat2
rslt$knots <- knotsuse
rslt$numknots <- numknotsuse
rslt$ms <- mslst
rslt$capms <- capms
rslt$cpar <- cpar
rslt$pvs <- pvs
rslt$s.edf <- s.edf
rslt$bstats <- bstats
#new:
rslt$pvsz <- pvsz
rslt$z.edf <- z.edf
rslt$fstats <- fstats
rslt$vh <- vh
rslt$kts.var <- kts.var
#rslt$sdhat2 <- sdhat2
return (rslt)
}
}
###########
#CicFamily#
###########
CicFamily <- function(object,...)UseMethod("CicFamily")
CicFamily <- function(object) {
#temp:
#if (object$family == "ordered") {
# object$family = "gaussian"
#}
#new:
if (object$family == "gaussian") {
linkfun <- function (mu) mu
}
if (object$family == "poisson") {
linkfun <- function (mu) log(mu)
}
if (object$family == "binomial") {
linkfun <- binomial()$linkfun
}
if (object$family == "Gamma") {
#linkfun <- function (mu) log(mu)
linkfun <- Gamma(link="log")$linkfun
}
llh.fun <- function(y, muhat = NULL, etahat = NULL, phihat = NULL, n = NULL, weights = NULL, fml = object$family){
sm <- 1e-7
#sm <- 1e-5
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
w <- weights
#new: avoid Inf
if (fml == "poisson") {
llh <- 2 * sum(w[w!=0] * (muhat[w!=0] - y[w!=0] * etahat[w!=0])) / n
}
if (fml == "binomial") {
llh <- 0
if (all(0 <= y) & all(y <= 1)) {
for (i in 1:n) {
if (muhat[i] > 0 & muhat[i] < 1) {
llh <- llh + w[i] * (y[i] * log(muhat[i]) + (1 - y[i]) * log(1 - muhat[i]))
}
}
llh <- (-2/n) * llh
} else {
stop ("y values must be 0 <= y <= 1!")
}
}
if (fml == "gaussian") {
if (all(w == 1)) {
llh <- log(sum((y - etahat)^2))
} else {
llh <- log(sum(w[w!=0] * (y[w!=0] - etahat[w!=0])^2)) - sum(log(w[w!=0])) / n
}
}
if (fml == "Gamma") {
vuhat <- 1 / phihat
#print (vuhat)
#llh <- 2 * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) / n
llh <- 2 / n * (vuhat * sum(w[w!=0] * (etahat[w!=0] + y[w!=0] * exp(-etahat[w!=0]))) + n * (log(gamma(vuhat)) - vuhat * log(vuhat)) - (vuhat-1) * sum(log(y)))
#print (llh)
}
llh
}
etahat.fun <- function(n, y, fml = object$family){
if (fml == "poisson") {
etahat <- 1:n*0 + log(mean(y))
}
if (fml == "binomial") {
etahat <- 1:n*0
}
if (fml == "Gamma") {
etahat <- 1:n*0 + log(mean(y))
}
etahat
}
gr.fun <- function(y, etahat = NULL, weights = NULL, fml = object$family){
n <- length(y)
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
w <- weights
if (fml == "poisson") {
gr <- w * (exp(etahat) - y)
}
if (fml == "binomial") {
if (all(etahat == 0)) {
gr <- w * (1/2 - y)
} else {
gr <- 1:n*0
for (i in 1:n) {
if (etahat[i] > 100) {
gr[i] <- w[i] * (1 - y[i])
} else {gr[i] <- w[i] * (exp(etahat[i]) / (1 + exp(etahat[i])) - y[i])}
}
}
}
if (fml == "Gamma") {
gr <- w * (1 - y * exp(-etahat))
}
gr
}
wt.fun <- function(y, etahat = NULL, n = NULL, weights = NULL, fml = object$family){
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
w <- weights
if (fml == "poisson") {
wt <- w * exp(etahat)
}
if (fml == "binomial") {
if (all(etahat == 0)){
#wt <- 1:n*0 + 1/4
wt <- w * (1:n*0 + 1/4)
} else {
wt <- 1:n*0
for (i in 1:n) {
if (etahat[i] > 100) {
wt[i] <- 0
} else {
wt[i] <- w[i] * exp(etahat[i]) / ((1 + exp(etahat[i]))^2)
}
}
}
}
if (fml == "gaussian") {
wt <- w # (1:n*0 + 1) / w
}
if (fml == "Gamma") {
wt <- w * y * exp(-etahat)
}
wt <- as.vector(wt)
wt
}
zvec.fun <- function(cvec = NULL, wt = NULL, y, sm = 1e-7, fml = object$family) {
n <- length(y)
if (fml == "gaussian") {
#zvec <- y
zvec <- wt^(1/2) * y
}
if (fml == "poisson") {
#zvec <- cvec / wt
zvec <- cvec / sqrt(wt)
}
if (fml == "binomial") {
zvec = 1:n*0
zvec[wt == 0] <- 1 / sm
zvec[wt > 0] <- cvec[wt > 0] / sqrt(wt[wt > 0])
}
if (fml == "Gamma") {
zvec <- cvec / sqrt(wt)
}
zvec
}
muhat.fun <- function(etahat, wt = NULL, fml = object$family){
n <- length(etahat)
if (fml == "poisson") {
muhat <- exp(etahat)
}
if (fml == "binomial") {
muhat <- 1:n*0
id1 = which(etahat > 100)
id2 = which(etahat <= 100)
muhat[id1] = 1
muhat[id2] = exp(etahat[id2]) / (1 + exp(etahat[id2]))
#muhat[wt == 0] <- 1
#for (i in 1:n) {
# if (etahat[i] > 100) {
# muhat[i] <- 1
# } else {
# muhat[i] <- exp(etahat[i]) / (1 + exp(etahat[i]))
# }
#}
}
if (fml == "gaussian") {
muhat <- etahat
}
if (fml == "Gamma") {
muhat <- exp(etahat)
}
muhat
}
ysim.fun <- function(n, mu0 = NULL, fml = object$family, shp0 = NULL) {
if (fml == "binomial") {
ysim <- 1:n*0
ysim[runif(n) < .5] <- 1
}
if (fml == "poisson") {
if (!is.null(mu0)) {
ysim <- rpois(n, mu0)
}
}
if (fml == "gaussian") {
ysim <- rnorm(n)
}
if (fml == "Gamma") {
ysim <- rgamma(n, shape=1)
}
ysim
}
deriv.fun <- function(muhat, fml = object$family) {
if (fml == "binomial") {
deriv <- 1 / (muhat * (1 - muhat))
}
if (fml == "poisson") {
deriv <- 1 / muhat
}
if (fml == "gaussian") {
deriv <- 1
}
if (fml == "Gamma") {
deriv <- 1 / muhat
}
deriv
}
dev.fun <- function(y, muhat, etahat, weights, fml = object$family){
n <- length(y)
sm <- 1e-7
#sm <- 1e-5
if (is.null(weights)) {
weights <- 1:n*0 + 1
}
w <- weights
vmat <- matrix(1:n*0 + 1, ncol = 1)
if (fml == "poisson") {
#dev <- 2 * sum(w * (y * log(y / muhat) - y + muhat))
dev <- 0
for (i in 1:n) {
if (y[i] == 0) {
dev <- dev + 2 * w[i] * muhat[i]
} else {
dev <- dev + 2 * w[i] * (y[i] * log(y[i] / muhat[i]) - y[i] + muhat[i])
}
}
}
if (fml == "binomial") {
dev <- 0
for (i in 1:n) {
if (y[i] == 0 & w[i] != 0) {
dev <- dev + 2 * w[i] * log(w[i] / (w[i] - w[i] * muhat[i]))
} else if (y[i] == 1 & w[i] != 0) {
dev <- dev + 2 * w[i] * log(w[i] / (w[i] * muhat[i]))
} else if (0 < y[i] & y[i] < 1 & w[i] != 0) {
dev <- dev + 2 * w[i] * y[i] * log(w[i] * y[i] / (w[i] * muhat[i])) + 2 * (w[i] - w[i] * y[i]) * log((w[i] - w[i] * y[i]) / (w[i] - w[i] * muhat[i]))
} else {
stop ("y values must be 0 <= y <= 1!")
}
}
}
if (fml == "gaussian") {
dev <- sum(w * (y - muhat)^2)
}
if (fml == "Gamma") {
dev <- 0
for (i in 1:n) {
if (y[i] == 0) {
sm <- 1e-4
dev <- dev + 2 * w[i] * ((sm - muhat[i]) / muhat[i] - log(sm / muhat[i]))
} else {
dev <- dev + 2 * w[i] * ((y[i] - muhat[i]) / muhat[i] - log(y[i] / muhat[i]))
}
}
}
###################
#get null deviance#
###################
if (fml == "binomial" | fml == "poisson" | fml == "Gamma") {
diff <- 1
muhat0 <- mean(y) + 1:n*0
if (fml == "poisson") {
etahat0 <- log(muhat0)
}
if (fml == "binomial") {
etahat0 <- log(muhat0 / (1 - muhat0))
}
if (fml == "Gamma") {
etahat0 <- log(muhat0)
}
while (diff > sm) {
oldmu <- muhat0
zhat <- etahat0 + (y - muhat0) * deriv.fun(muhat0, fml = fml)
# wmat <- diag(as.vector(w / deriv.fun(muhat0, fml = fml)))
# b <- solve(t(vmat) %*% wmat %*% vmat) %*% t(vmat) %*% wmat %*% zhat
wm <- as.vector(w / deriv.fun(muhat0, fml = fml))
tvmat <- t(vmat)
for (i in 1:n) {tvmat[,i] <- tvmat[,i] * wm[i]}
b <- solve(tvmat %*% vmat) %*% tvmat %*% zhat
etahat0 <- vmat %*% b
muhat0 <- muhat.fun(etahat0, fml = fml)
diff <- mean((muhat0 - oldmu)^2)
}
if (fml == "poisson") {
#dev.null <- 2 * sum(w * (y * log(y / muhat0) - y + muhat0))
dev.null <- 0
for (i in 1:n) {
if (y[i] == 0) {
dev.null <- dev.null + 2 * w[i] * muhat0[i]
} else {
dev.null <- dev.null + 2 * w[i] * (y[i] * log(y[i] / muhat0[i]) - y[i] + muhat0[i])
}
}
}
if (fml == "binomial") {
dev.null <- 0
for (i in 1:n) {
if (y[i] == 0 & w[i] != 0) {
dev.null <- dev.null + 2 * w[i] * log(w[i] / (w[i] - w[i] * muhat0[i]))
} else if (y[i] == 1 & w[i] != 0) {
dev.null <- dev.null + 2 * w[i] * log(w[i] / (w[i] * muhat0[i]))
} else if (0 < y[i] & y[i] < 1 & w[i] != 0) {
dev.null <- dev.null + 2 * w[i] * y[i] * log(w[i] * y[i] / (w[i] * muhat0[i])) + 2 * (w[i] - w[i] * y[i]) * log((w[i] - w[i] * y[i]) / (w[i] - w[i] * muhat0[i]))
} else {
stop ("y values must be 0 <= y <= 1!")
}
}
}
if (fml == "Gamma") {
dev.null <- 0
for (i in 1:n) {
if (y[i] == 0) {
sm <- 1e-4
dev.null <- dev.null + 2 * w[i] * ((sm - muhat0[i]) / muhat0[i] - log(sm / muhat0[i]))
} else {
dev.null <- dev.null + 2 * w[i] * ((y[i] - muhat0[i]) / muhat0[i] - log(y[i] / muhat0[i]))
}
}
}
}
if (fml == "gaussian") {
# wmat <- diag(w)
# b <- solve(t(vmat) %*% wmat %*% vmat) %*% t(vmat) %*% wmat %*% y
tvmat <- t(vmat)
for (i in 1:n) {tvmat[,i] <- tvmat[,i] * w[i]}
b <- solve(tvmat %*% vmat) %*% tvmat %*% y
etahat0 <- vmat %*% b
muhat0 <- muhat.fun(etahat0, fml = fml)
dev.null <- sum(w * (y - muhat0)^2)
}
rslt <- new.env()
rslt$dev <- dev
rslt$dev.null <- dev.null
rslt
}
ans <- list(llh.fun = llh.fun, etahat.fun = etahat.fun, gr.fun = gr.fun, wt.fun = wt.fun, zvec.fun = zvec.fun, muhat.fun = muhat.fun, ysim.fun = ysim.fun, deriv.fun = deriv.fun, dev.fun = dev.fun, linkfun = linkfun)
class(ans) <- "CicFamily"
return (ans)
}
#######################
#eight shape functions#
#######################
incr <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 1
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
decr <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 2
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 3
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 4
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
incr.conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 5
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
decr.conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 6
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
incr.conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 7
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
decr.conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 8
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
#s.incr <- function(x, numknots = 0, knots = 0, space = "E")
#{
# cl <- match.call()
# pars <- match.call()[-1]
# attr(x, "nm") <- deparse(pars$x)
# attr(x, "shape") <- 9
# attr(x, "numknots") <- numknots
# attr(x, "knots") <- knots
# attr(x, "space") <- space
# attr(x, "categ") <- "additive"
# #class(x) <- "additive"
# x
#}
#s.decr <- function(x, numknots = 0, knots = 0, space = "E")
#{
# cl <- match.call()
# pars <- match.call()[-1]
# attr(x, "nm") <- deparse(pars$x)
# attr(x, "shape") <- 10
# attr(x, "numknots") <- numknots
# attr(x, "knots") <- knots
# attr(x, "space") <- space
# attr(x, "categ") <- "additive"
# #class(x) <- "additive"
# x
#}
s.incr <- function(x, numknots = 0, knots = 0, var.knots = 0, space = "E", db.exp = FALSE)
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 9
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
attr(x, "db.exp") <- db.exp
attr(x, "var.knots") <- var.knots
#class(x) <- "additive"
x
}
s.decr <- function(x, numknots = 0, knots = 0, var.knots = 0, space = "E", db.exp = FALSE)
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 10
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
attr(x, "db.exp") <- db.exp
attr(x, "var.knots") <- var.knots
#class(x) <- "additive"
x
}
s.conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 11
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s.conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 12
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s.incr.conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 13
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s.incr.conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 14
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s.decr.conv <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 15
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s.decr.conc <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 16
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
s <- function(x, numknots = 0, knots = 0, space = "E")
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- 17
attr(x, "numknots") <- numknots
attr(x, "knots") <- knots
attr(x, "space") <- space
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
######################################################
#tree function: give the tree shape to x and return x#
######################################################
tree <- function(x, pl = NULL)
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- "tree"
#stop (print (x))
#print (class(x))
if (is.null(pl)) {
if (is.numeric(x)) {
if (0%in%x) {
pl <- 0
} else {pl <- min(x)}
} else {
xu <- unique(x)
pl <- xu[1]
}
} else {
if (!(pl%in%x)) {
stop ("placebo level is not a level of the tree variable!")
}
}
attr(x, "pl") <- pl
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
##################################################
#tree.fun: make delta to a tree ordering variable#
##################################################
#tree.fun <- function(x)
#{
# if (min(x) != 0) {
# stop ("A tree ordering variable must have its placebo equal to 0!")
# }
# if (!all(round(x, 0) == x)) {
# stop ("All elements of a tree ordering variable must be integers!")
# }
# if (any(x < 0))
# stop ("All elements of a tree ordering variable must be positive!")
# nx <- x
# obs <- 1:length(x)
# delta <- matrix(0, nrow = length(attributes(factor(x))$levels) - 1, ncol = length(x))
# pl <- min(nx)
# for (i in 1:nrow(delta)) {
# nx <- nx[which(nx != pl)]
# pl <- min(nx)
# index <- obs[x == pl]
# delta[i, index] <- 1
# }
# attr(delta, "shape") <- "tree"
# delta
#}
tree.fun <- function(x, pl = NULL)
{
#if (pl %in% x) {
if (is.null(pl)) {
if (is.numeric(x)) {
if (0%in%x) {
pl <- 0
} else {pl <- min(x)}
} else {
xu <- unique(x)
pl <- xu[1]
}
} else {
if (!(pl%in%x)) {
stop ("placebo level is not a level of the tree variable!")
}
}
pos <- x %in% pl
xu <- unique(x)
nx <- xu[xu != pl]
delta <- matrix(0, nrow = (length(xu) - 1), ncol = length(x))
for (i in 1:nrow(delta)) {
xi = nx[i]
delta[i, !pos & x %in% xi] <- 1
pos <- pos | x %in% xi
}
#} else {stop ("placebo level is not a level of the tree variable!")}
attr(delta, "shape") <- "tree"
delta
}
##############################################################
#umbrella function: give the umbrella shape to x and return x#
##############################################################
umbrella <- function(x)
{
cl <- match.call()
pars <- match.call()[-1]
attr(x, "nm") <- deparse(pars$x)
attr(x, "shape") <- "umbrella"
attr(x, "categ") <- "additive"
#class(x) <- "additive"
x
}
###########################################################
#umbrella.fun: make delta to an umbrella ordering variable#
###########################################################
umbrella.fun <- function(x)
{
amat <- amat.fun(x)
bmat <- bmat.fun(x)
constr <- t(rbind(amat, bmat))
vmat <- qr.Q(qr(constr), complete = TRUE)[, -(1:(qr(constr)$rank)), drop = FALSE]
if (!is.null(bmat)) {
wperp <- t(rbind(t(vmat), bmat))
wmat <- qr.Q(qr(wperp), complete = TRUE)[, -(1:(qr(wperp)$rank)), drop = FALSE]
atil <- amat %*% wmat
delta <- t(wmat %*% t(atil) %*% solve(atil %*% t(atil)))
} else {
delta <- t(t(amat) %*% solve(amat %*% t(amat)))
}
attr(delta, "shape") <- "umbrella"
delta
}
###################################################
#find delta for a specific predictor x and a shape#
###################################################
###################################################
#find delta for a specific predictor x and a shape#
###################################################
makedelta = function(x, sh, numknots = 0, knots = 0, space = "E", suppre = FALSE, interp = FALSE) {
#if (!interp) {
#x = (x - min(x)) / (max(x) - min(x))
#}
n = length(x)
# find unique x values
#round(x,8) will make 0 edge in amat!
#xu = sort(unique(round(x, 8)))
#new: center and scale to avoid numerical instabillity
xu = sort(unique(x))
n1 = length(xu)
sm = 1e-7
ms = NULL
# increasing or decreasing
if (sh < 3) {
amat = matrix(0, nrow = n1 - 1, ncol = n)
for (i in 1: (n1 - 1)) {
amat[i, x > xu[i]] = 1
}
if (sh == 2) {amat = -amat}
if (!interp) {
for (i in 1:(n1 - 1)) {
#new: use ms in predict.cgam
ms = c(ms, mean(amat[i, ]))
amat[i, ] = amat[i, ] - mean(amat[i, ])
}
}
} else if (sh == 3 | sh == 4) {
# convex or concave
amat = matrix(0, nrow = n1 - 2 ,ncol = n)
#for (i in 1: (n1 - 2)) {
# amat[i, x > xu[i]] = x[x > xu[i]] - xu[i]
#}
for (i in 1: (n1 - 2)) {
amat[i, x > xu[i+1]] = x[x > xu[i+1]] - xu[i+1]
}
if (sh == 4) {amat = -amat}
xm = cbind(1:n*0+1,x)
xpx = solve(t(xm) %*% xm)
pm = xm %*% xpx %*% t(xm)
#new: use ms in predict.cgam
if (!interp) {
ms = amat %*% t(pm)
#amat = amat - amat %*% t(pm)
amat = amat - ms
}
} else if (sh > 4 & sh < 9) {
amat = matrix(0, nrow = n1 - 1, ncol = n)
if (sh == 5) { ### increasing convex
for (i in 1:(n1 - 1)) {
amat[i, x > xu[i]] = (x[x > xu[i]] - xu[i]) / (max(x) - xu[i])
}
if (!interp) {
for (i in 1:(n1 - 1)) {
ms = c(ms, mean(amat[i, ]))
amat[i,] = amat[i,] - mean(amat[i,])
}
}
} else if (sh == 6) { ## decreasing convex
for (i in 1:(n1 - 1)) {
amat[i, x < xu[i + 1]] = (x[x < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
}
if (!interp) {
for (i in 1:(n1 - 1)) {
ms = c(ms, mean(amat[i, ]))
amat[i,] = amat[i,] - mean(amat[i, ])
}
}
#print (ms)
} else if (sh == 7) { ## increasing concave
for (i in 1:(n1 - 1)) {
amat[i, x < xu[i + 1]] = (x[x < xu[i + 1]] - xu[i + 1]) / (min(x) - xu[i + 1])
}
if (!interp) {
for (i in 1:(n1 - 1)) {
ms = c(ms, mean(amat[i, ]))
amat[i,] = -amat[i,] + mean(amat[i,])
}
}
} else if (sh == 8) {## decreasing concave
for (i in 1:(n1 - 1)) {
amat[i, x > xu[i]] = (x[x > xu[i]] - xu[i]) / (max(x) - xu[i])
}
if (!interp) {
for (i in 1:(n1 - 1)) {
ms = c(ms, mean(amat[i, ]))
amat[i,] = -amat[i,] + mean(amat[i,])
}
}
}
} else if (sh > 8 & sh < 18) {
#new: add two knots
#if (all(knots == 0) & numknots == 0) {
if (length(knots) < 2 & numknots == 0) {
if (sh == 9 | sh == 10) {#1 2
#k = trunc(n1^(1/5)) + 4
if (n1 <= 50) {
k = 5
} else if (n1>50 && n1<100) {
k = 6
} else if (n1>= 100 && n1<200) {
k = 7
} else {
k = trunc(n1^(1/5)) + 6
}
} else {
#k = trunc(n1^(1/7) + 4)
if (n1 <= 50) {
k = 5
} else if (n1>50 && n1<100) {
k = 6
} else if (n1>= 100 && n1<200) {
k = 7
} else {
k = trunc(n1^(1/7)) + 6
}
}
if (space == "Q") {
t = quantile(xu, probs = seq(0, 1, length = k), names = FALSE)
}
if (space == "E") {
#t = 0:k / k * (max(x) - min(x)) + min(x)
t = 0:(k-1) / (k-1) * (max(x) - min(x)) + min(x)
}
#} else if (any(knots != 0) & numknots == 0) {
} else if (length(knots) >= 2 & numknots == 0) {
t = knots
#} else if (all(knots == 0) & numknots != 0) {
} else if (length(knots) < 2 & numknots != 0) {
if (space == "Q") {
t = quantile(xu, probs = seq(0, 1, length = numknots), names = FALSE)
}
if (space == "E") {
#k = numknots
#new: numknots should be the # of all knots
k = numknots - 1
#if (sh == 9 | sh == 10) {#1 2
# k = trunc(n1^(1/5)) + 4
#} else {k = trunc(n1^(1/7) + 4)}
t = 0:k / k * (max(x) - min(x)) + min(x)
}
#} else if (any(knots != 0) & numknots != 0) {
} else if (length(knots) >= 2 & numknots != 0) {
#t0 = quantile(xu, probs = seq(0, 1, length = numknots), names = FALSE)
t = knots
if (!suppre) {
print("'knots' is used! 'numknots' is not used!")
}
#print ("'knots' is used!")
#if (numknots != length(knots)) {
# if (!suppre) {
# print("length(knots) is not equal to 'numknots'! 'knots' is used!")
# }
#} else if (any(t0 != knots)) {
# if (!suppre) {
# print("equal x-quantiles knots != 'knots'! 'knots' is used! ")
# }
#}
}
if (sh == 9) {#1
amat_ans = monincr(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 10) {#2
amat_ans = mondecr(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 11) {#3
amat_ans = convex(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 12) {#4
amat_ans = concave(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 13) {#5
amat_ans = incconvex(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 14) {#6
amat_ans = incconcave(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
} else if (sh == 15) {#7
#amat_ans = -incconcave(x, t, interp)
amat_ans = incconcave(x, t, interp)
amat = -amat_ans$sigma
if (!interp) {
ms = -amat_ans$ms
}
} else if (sh == 16) {#8
#amat_ans = -incconvex(x, t, interp)
amat_ans = incconvex(x, t, interp)
amat = -amat_ans$sigma
if (!interp) {
ms = -amat_ans$ms
}
} else if (sh == 17) {#unconstrained
amat_ans = incconvex(x, t, interp)
amat = amat_ans$sigma
ms = amat_ans$ms
#amat = -incconcave(x, t)
#amat = rbind(x, t(bcspl(x, m = length(t), knots = t)$bmat))
#amat = rbind(x, convex(x, t))
}
}
#if (sh < 9) {
# rslt = list(amat = amat, knots = 0, ms = ms)
#} else {
# rslt = list(amat = amat, knots = t, ms = ms)
#}
if (sh < 9) {t = 0}
rslt = list(amat = amat, knots = t, ms = ms)
rslt
}
#####
#make basis for additive components
make_delta_add = function(xmat, shapes, numknots, knots, space) {
capl = ncol(xmat)
if (capl < 1) {capl = 0}
if (round(capl, 8) != round(capl, 1)) {stop ("Incompatible dimensions for xmat!")}
capls = sum(shapes == 17)
delta = NULL
varlist = NULL
knotsuse = list(); numknotsuse = NULL
mslst = list()
capm = 0
capms = 0
del1_ans = makedelta(xmat[, 1], shapes[1], numknots[1], knots[[1]], space = space[1], interp = FALSE)
del1 = del1_ans$amat
knotsuse[[1]] = del1_ans$knots
mslst[[1]] = del1_ans$ms
numknotsuse = c(numknotsuse, length(del1_ans$knots))
m1 = nrow(del1)
#new code: record the number of columns of del1 if shapes0[1] == 17:
if (shapes[1] == 17) {capms = capms + m1}
var1 = 1:m1*0 + 1
if (capl == 1) {
delta = del1
varlist = var1
} else {
for (i in 2:capl) {
del2_ans = makedelta(xmat[,i], shapes[i], numknots[i], knots[[i]], space = space[i], interp = FALSE)
del2 = del2_ans$amat
knotsuse[[i]] = del2_ans$knots
mslst[[i]] = del2_ans$ms
numknotsuse = c(numknotsuse, length(del2_ans$knots))
m2 = nrow(del2)
#new code: record the number of columns of del2 if shapes0[i] == 17:
if (shapes[i] == 17) {capms = capms + m2}
delta = rbind(del1, del2)
varlist = 1:(m1 + m2)*0
varlist[1:m1] = var1
varlist[(m1 + 1):(m1 + m2)] = (1:m2)*0 + i
var1 = varlist
m1 = m1 + m2
del1 = delta
}
}
if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) > 0) {
bigmat = rbind(t(xmat[, shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13]), delta)
np = sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) + capms
} else if (sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13) == 0) {
bigmat = delta
np = capms
}
rslt = list(bigmat=bigmat, np=np, varlist=varlist)
return(rslt)
}
# Monotone increasing
monincr = function(xs, t, interp = FALSE) {
n = length(xs)
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
k = length(t) - 2
m = k + 2
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
knt = 1:m
for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
for (j in 1:(k-1)) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = 0
index = x > t[j] & x <= t[j+1]
sigma[j, index] = (x[index] - t[j])^2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = 1 - (x[index] - t[j+2])^2 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+2] #& x <= t[m]
sigma[j, index] = 1
}
index = x >= t[1] & x <= t[k]
sigma[k, index] = 0
index = x > t[k] & x <= t[k+1]
sigma[k, index] = (x[index] - t[k])^2 / (t[k+2] - t[k]) / (t[k+1] - t[k])
index = x > t[k+1] & x <= t[k+2]
sigma[k, index] = 1 - (x[index] - t[k+2])^2 / (t[k+2] - t[k+1]) / (t[k+2] - t[k])
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = 1 - (t[2] - x[index])^2 / (t[2] - t[1])^2
index = x > t[2]
sigma[k+1, index] = 1
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = 0
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = (x[index] - t[k+1])^2 / (t[k+2] - t[k+1])^2
#new:
ms = NULL
if (!interp) {
ms = apply(sigma, 1, mean)
for (i in 1:m) {
sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
########################################################
# Monotone decreasing
mondecr = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
n = length(x)
k = length(t) - 2
m = k + 2
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
#knt = 1:m
#for (i in 1:(k + 2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
#t = x[knt]
for (j in 1:(k - 1)) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = 1
index = x > t[j] & x <= t[j+1]
sigma[j, index] = 1 - (x[index] - t[j])^2 / (t[j+2] - t[j]) / (t[j+1] - t[j])
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = (x[index] - t[j+2])^2 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+2]
sigma[j, index] = 0
}
index = x >= t[1] & x <= t[k]
sigma[k, index] = 1
index = x > t[k] & x <= t[k+1]
sigma[k, index] = 1 - (x[index] - t[k])^2 / (t[k+2] - t[k]) / (t[k+1] - t[k])
index = x > t[k+1] & x <= t[k+2]
sigma[k, index] = (x[index] - t[k+2])^2 / (t[k+2] - t[k+1]) / (t[k+2] - t[k])
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = (t[2] - x[index])^2 / (t[2] - t[1])^2
index = x > t[2]
sigma[k+1, index] = 0
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = 1
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = 1 - (x[index] - t[k+1])^2 / (t[k+2] - t[k+1])^2
ms = NULL
if (!interp) {
ms = apply(sigma, 1, mean)
for (i in 1:m) {
sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
########################################################
# Convex
convex = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
n = length(x)
k = length(t) - 2
m = k + 2
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
#knt = 1:m
#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
for (j in 1:(k-1)) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = 0
index = x > t[j] & x <= t[j+1]
sigma[j, index] = (x[index] - t[j])^3 / (t[j+2] - t[j]) / (t[j+1] - t[j]) / 3
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = x[index] - t[j+1] - (x[index] - t[j+2])^3 / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) / 3 + (t[j+1] - t[j])^2 / 3 /(t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
index = x > t[j+2]
sigma[j, index] = (x[index] - t[j+1]) + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
}
index = x >= t[1] & x <= t[k]
sigma[k, index] = 0
index = x > t[k] & x <= t[k+1]
sigma[k, index] = (x[index] - t[k])^3 / (t[k+2] - t[k]) / (t[k+1] - t[k]) / 3
index = x > t[k+1] & x <= t[k+2]
sigma[k, index] = x[index] - t[k+1] - (x[index] - t[k+2])^3 / (t[k+2] - t[k]) / (t[k+2] - t[k+1]) / 3 + (t[k+1] - t[k])^2 / 3 / (t[k+2] -t[k]) - (t[k+2] - t[k+1])^2 / 3 / (t[k+2] - t[k])
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = x[index] - t[1] + (t[2] - x[index])^3 / (t[2] - t[1])^2 / 3 - (t[2] - t[1]) / 3 #-(t[2]-t[1])^3/(t[2]-t[1])^2/3 #
index = x > t[2]
sigma[k+1, index] = x[index] - t[1] - (t[2] - t[1]) / 3 #-(t[2]-t[1])^3/(t[2]-t[1])^2/3#
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = 0
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = (x[index] - t[k+1])^3 / (t[k+2] - t[k+1])^2 / 3
ms = NULL
if (!interp) {
xm = cbind(1:n*0+1, x)
pm = xm %*% solve(t(xm) %*% xm) %*% t(xm)
ms = matrix(0, nrow = nrow(sigma), ncol = ncol(sigma))
for (i in 1:m) {
ms[i,] = pm %*% sigma[i,]
ms[i,] = ms[i, rank(xs)]
#rng=max(sigma[i,])
#sigma[i,]=sigma[i,]/rng
sigma[i,] = sigma[i,] - pm %*% sigma[i,]
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - pm %*% sigma[i,]
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
########################################################
# Concave
concave = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
n = length(x)
k = length(t) - 2
m = k + 2
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
#knt = 1:m
#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
#t = x[knt]
for (j in 1:k) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = x[index] - t[1]
index = x > t[j] & x <= t[j+1]
sigma[j, index] = t[j] - t[1] + ((t[j+1] - t[j])^3 - (t[j+1] - x[index])^3) / 3 / (t[j+1] - t[j]) / (t[j+2] - t[j]) + (x[index] - t[j]) * (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + ((t[j+2] - t[j+1])^3 - (t[j+2] - x[index])^3) / 3 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+2]
sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
}
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = -(t[2] - x[index])^3 / 3 / (t[2] - t[1])^2
index = x > t[2]
sigma[k+1, index] = 0
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = x[index] - t[1]
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = t[k+1] - t[1] + ((t[k+2] - t[k+1])^2 * (x[index] - t[k+1]) - (x[index] - t[k+1])^3 / 3) / (t[k+2] - t[k+1])^2
ms = NULL
if (!interp) {
xm = cbind(1:n*0+1, x)
pm = xm %*% solve(t(xm) %*% xm) %*% t(xm)
ms = matrix(0, nrow = nrow(sigma), ncol = ncol(sigma))
for (i in 1:m) {
ms[i,] = pm %*% sigma[i,]
ms[i,] = ms[i, rank(xs)]
sigma[i,] = sigma[i,] - pm %*% sigma[i,]
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - pm %*% sigma[i,]
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
########################################################
# Increasing and Convex
incconvex = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
n = length(x)
k = length(t) - 2
m = k + 3
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
#knt = 1:(k+2)
#for (i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
for (j in 1:(k-1)) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = 0
index = x > t[j] & x <= t[j+1]
sigma[j, index] = (x[index] - t[j])^3 / (t[j+2] - t[j]) / (t[j+1] - t[j]) / 3
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = x[index] - t[j+1] - (x[index] - t[j+2])^3 / (t[j+2] - t[j]) / (t[j+2] - t[j+1]) / 3 + (t[j+1] - t[j])^2 / 3 /(t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
index = x > t[j+2]
sigma[j, index] = (x[index] - t[j+1]) + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) - (t[j+2] - t[j+1])^2 / 3 / (t[j+2] - t[j])
}
index = x >= t[1] & x <= t[k]
sigma[k, index] = 0
index = x > t[k] & x <= t[k+1]
sigma[k, index] = (x[index] - t[k])^3 / (t[k+2] - t[k]) / (t[k+1] - t[k]) / 3
index = x > t[k+1] & x <= t[k+2]
sigma[k, index] = x[index] - t[k+1] - (x[index] - t[k+2])^3 / (t[k+2] - t[k]) / (t[k+2] - t[k+1]) / 3 + (t[k+1] - t[k])^2 / 3 / (t[k+2] - t[k]) - (t[k+2] - t[k+1])^2 / 3 / (t[k+2] - t[k])
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = x[index] - t[1] + (t[2] - x[index])^3 / (t[2] - t[1])^2 / 3 - (t[2] - t[1]) / 3
index = x > t[2]
sigma[k+1, index] = x[index] - t[1] - (t[2] - t[1]) / 3
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = 0
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = (x[index] - t[k+1])^3 / (t[k+2] - t[k+1])^2 / 3
sigma[k+3,] = x
ms = NULL
if (!interp) {
ms = apply(sigma, 1, mean)
for (i in 1:m) {
sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
########################################################
# Increasing and Concave
incconcave = function(xs, t, interp = FALSE) {
#xs = (xs - min(xs)) / (max(xs) - min(xs))
x = sort(xs)
n = length(x)
k = length(t) - 2
m = k + 3
sigma = matrix(0, nrow = m, ncol = n)
obs = 1:n
#knt = 1:(k+2)
#for(i in 1:(k+2)) {knt[i] = min(obs[abs(x - t[i]) == min(abs(x - t[i]))])}
for (j in 1:k) {
index = x >= t[1] & x <= t[j]
sigma[j, index] = x[index] - t[1]
index = x > t[j] & x <= t[j+1]
sigma[j, index] = t[j] - t[1] + ((t[j+1] - t[j])^3 - (t[j+1] - x[index])^3) / 3 / (t[j+1] - t[j]) / (t[j+2] - t[j]) + (x[index] - t[j]) * (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+1] & x <= t[j+2]
sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + ((t[j+2] - t[j+1])^3 - (t[j+2] - x[index])^3) / 3 / (t[j+2] - t[j+1]) / (t[j+2] - t[j])
index = x > t[j+2]
sigma[j, index] = t[j] - t[1] + (t[j+1] - t[j])^2 / 3 / (t[j+2] - t[j]) + (t[j+2] - t[j+1]) * (t[j+1] - t[j]) / (t[j+2] - t[j]) + (t[j+2] - t[j+1])^2 / 3 /(t[j+2] - t[j])
}
index = x >= t[1] & x <= t[2]
sigma[k+1, index] = -(t[2] - x[index])^3 / 3 / (t[2] - t[1])^2
index = x > t[2]
sigma[k+1, index] = 0
index = x >= t[1] & x <= t[k+1]
sigma[k+2, index] = x[index] - t[1]
index = x > t[k+1] & x <= t[k+2]
sigma[k+2, index] = t[k+1] - t[1] + ((t[k+2] - t[k+1])^2 * (x[index] - t[k+1]) - (x[index] - t[k+1])^3 / 3) / (t[k+2] - t[k+1])^2
sigma[k+3, ] = x
ms = NULL
if (!interp) {
ms = apply(sigma, 1, mean)
for (i in 1:m) {
sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
} else {
for (i in 1:m) {
#sigma[i,] = sigma[i,] - mean(sigma[i,])
sigma[i,] = sigma[i, rank(xs)]
}
}
rslt = list(sigma = sigma, ms = ms)
rslt
}
##############
#summary.cgam#
##############
summary.cgam <- function(object,...) {
if (!is.null(object$zcoefs)) {
family <- object$family
resid_df_obs <- object$resid_df_obs
wt.iter <- object$wt.iter
coefs <- object$zcoefs
se <- object$se.beta
tval <- coefs / se
pvalbeta <- object$pvals.beta
n <- length(coefs)
sse0 <- object$SSE0
sse1 <- object$SSE1
cic <- object$cic
deviance <- object$deviance
null_deviance <- object$null_deviance
df <- object$df
null_df <- object$null_df
zid <- object$zid
#zid1 <- object$zid1 - 1 - length(shapes)
#zid2 <- object$zid2 - 1 - length(shapes)
#new: zid1, zid2 just index zmat not bigmat
zid1 <- object$zid1
zid2 <- object$zid2
tms <- object$tms
is_param <- object$is_param
is_fac <- object$is_fac
vals <- object$vals
pvs <- object$pvs
s.edf <- object$s.edf
bstats <- object$bstats
if (wt.iter) {
rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "z.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
rownames(rslt1)[1] <- "(Intercept)"
if (n > 1) {
lzid <- length(zid1)
for (i in 1:lzid) {
pos1 <- zid1[i]; pos2 <- zid2[i]
for (j in pos1:pos2) {
if (!is_param[i]) {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")
} else {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")
}
}
}
}
rslt1 <- as.matrix(rslt1)
#new:
rslt2 <- NULL
if (!is.null(pvs)) {
rslt2 <- data.frame("edf" = round(s.edf, 4), "mixture of Beta" = round(bstats, 4), "p.value" = round(pvs, 4))
#rownames(rslt2) <- attributes(tms)$term.labels
#debugged: check more
if (!is.null(zid)) {
rownames(rslt2) <- (attributes(tms)$term.labels)[-(zid-1)]
} else {
rownames(rslt2) <- (attributes(tms)$term.labels)
}
}
} else {
rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "t.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
rownames(rslt1)[1] <- "(Intercept)"
if (n > 1) {
lzid <- length(zid1)
for (i in 1:lzid) {
pos1 <- zid1[i]; pos2 <- zid2[i];
for (j in pos1:pos2) {
if (!is_param[i]) {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")
} else {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")
}
}
}
}
rslt1 <- as.matrix(rslt1)
rslt2 <- NULL
if (!is.null(pvs)) {
rslt2 <- data.frame("edf" = round(s.edf, 4), "mixture of Beta" = round(bstats, 4), "p.value" = round(pvs, 4))
#debugged: check more
if (!is.null(zid)) {
rownames(rslt2) <- (attributes(tms)$term.labels)[-(zid-1)]
} else {
rownames(rslt2) <- (attributes(tms)$term.labels)
}
}
}
#if (!is.null(sse0) & !is.null(sse1)) {
#rslt2 <- cbind(SSE.Linear = sse0, SSE.Full = sse1)
#new:
# rslt2 <- data.frame("SSE.Linear" = sse0, "SSE.Full" = sse1)
# rownames(rslt2)[1] <- ""
# ans <- list(call = object$call, coefficients = rslt1, residuals = rslt2, zcoefs = coefs, cic = cic, null_deviance = null_deviance, null_df = null_df, deviance = deviance, df = df, resid_df_obs = resid_df_obs, family = family)
# class(ans) <- "summary.cgam"
# ans
#} else {
ans <- list(call = object$call, coefficients = rslt1, coefficients2 = rslt2, zcoefs = coefs, cic = cic, null_deviance = null_deviance, null_df = null_df, deviance = deviance, df = df, resid_df_obs = resid_df_obs, family = family)
class(ans) <- "summary.cgam"
ans
#}
} else {
ans <- list(zcoefs = object$zcoefs)
class(ans) <- "summary.cgam"
ans
}
}
#############
#summary.wps#
#############
summary.wps <- function(object,...) {
if (!is.null(object$zcoefs)) {
coefs <- object$zcoefs
se <- object$se.beta
#tval <- object$tz
pvalbeta <- object$pvals.beta
tval <- coefs / se
n <- length(coefs)
#sse0 <- object$SSE0
#sse1 <- object$SSE1
zid <- object$zid
#new: zid1, zid2 just index zmat not bigmat
zid1 <- object$zid1
zid2 <- object$zid2
#tms <- object$tms
#zmat <- object$zmat
#is_mat <- object$is_mat
is_param <- object$is_param
is_fac <- object$is_fac
vals <- object$vals
tms <- object$tms
#new:
cic <- object$cic
rslt1 <- data.frame("Estimate" = round(coefs, 4), "StdErr" = round(se, 4), "t.value" = round(tval, 4), "p.value" = round(pvalbeta, 4))
rownames(rslt1)[1] <- "(Intercept)"
if (n > 1) {
lzid <- length(zid1)
for (i in 1:lzid) {
pos1 <- zid1[i]; pos2 <- zid2[i]
for (j in pos1:pos2) {
if (!is_param[i]) {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], rownames(rslt1)[j + 1], sep = "")
} else {
rownames(rslt1)[j + 1] <- paste(attributes(tms)$term.labels[zid[i] - 1], vals[j], sep = "")
}
}
}
}
rslt1 <- as.matrix(rslt1)
#if (!is.null(sse0) & !is.null(sse1)) {
# rslt2 <- data.frame("SSE.Linear" = sse0, "SSE.Full" = sse1)
#new:
# rownames(rslt2)[1] <- ""
#rslt2 <- as.matrix(rslt2)
# ans <- list(call = object$call, coefficients = rslt1, residuals = rslt2, zcoefs = coefs)
# class(ans) <- "summary.wps"
# ans
#} else {
ans <- list(call = object$call, coefficients = rslt1, zcoefs = coefs, cic = cic)
class(ans) <- "summary.wps"
ans
#}
} else {
ans <- list(zcoefs = object$zcoefs)
class(ans) <- "summary.wps"
ans
}
}
#####################
#print.summary.cgam #
#####################
print.summary.cgam <- function(x,...) {
if (!is.null(x$zcoefs)) {
#if (!is.null(x$se.beta)) {
cat("Call:\n")
print(x$call)
cat("\n")
cat("Coefficients:")
cat("\n")
printCoefmat(x$coefficients, P.values = TRUE, has.Pvalue = TRUE)
cat("\n")
if (x$family$family == "binomial") {
cat("(Dispersion parameter for binomial family taken to be 1)", "\n")
}
if (x$family$family == "poisson") {
cat("(Dispersion parameter for poisson family taken to be 1)", "\n")
}
if (x$family$family == "gaussian") {
cat("(Dispersion parameter for gaussian family taken to be ", round(x$deviance/x$df,4),")","\n", sep="")
}
cat("\n")
cat("Null deviance: ",round(x$null_deviance,4), "", "on", x$null_df, "", "degrees of freedom", "\n")
cat("Residual deviance: ",round(x$deviance,4), "", "on", x$resid_df_obs, "", "observed degrees of freedom", "\n")
#if (is.null(x$cic)) {
# message("CIC value is not available when there is no shape-restricted predictor")
#} else {message("CIC: ", round(x$cic,4))}
if (!is.null(x$coefficients2)) {
cat("\n")
cat("Approximate significance of smooth terms: \n")
printCoefmat(x$coefficients2, P.values = TRUE, has.Pvalue = TRUE)
}
if (!is.null(x$cic)) {
cat("CIC: ", round(x$cic,4))
}
#if (!is.null(x$residuals)) {
# cat("==============================================================", "\n")
# cat("Call:\n")
# print(x$call)
# cat("\n")
# printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
#}
} else {
print ("No linear predictor is defined")
#print ("Residual degree of freedom is negative!")
}
}
####################
#print.summary.wps #
####################
print.summary.wps <- function(x,...) {
if (!is.null(x$zcoefs)) {
#if (!is.null(x$se.beta)) {
cat("Call:\n")
print(x$call)
cat("\n")
cat("Coefficients:")
cat("\n")
printCoefmat(x$coefficients, P.values = TRUE, has.Pvalue = TRUE)
#if (!is.null(x$residuals)) {
#cat("==============================================================", "\n")
# cat("+--------------------------------------+\n")
# cat("| SSE.Linear vs SSE.Full |\n")
#printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
# cat("+--------------------------------------+\n")
#cat("Call:\n")
#print(x$call)
#cat("\n")
# printCoefmat(x$residuals, P.values = TRUE, has.Pvalue = TRUE)
#}
cat("\n")
if (!is.null(x$cic)) {
cat("CIC: ", round(x$cic,4), "\n", sep = "")
}
} else {
print ("No linear predictor is defined")
}
}
##############
#predict.cgam#
##############
predict.cgam = function(object, newData, interval = c("none", "confidence", "prediction"), type = c("response", "link"), level = 0.95, n.mix = 500,...) {
#print (is.data.frame(newData))
#print (newData)
#new:
family = object$family
cicfamily = CicFamily(family)
muhat.fun = cicfamily$muhat.fun
if (!inherits(object, "cgam")) {
warning("calling predict.cgam(<fake-cgam-object>) ...")
}
if (missing(newData) || is.null(newData)) {
#if (missing(newData)) {
etahat = object$etahat
muhat = muhat.fun(etahat, fml = family$family)
#ans = list(fit = muhat, etahat = etahat, newbigmat = object$bigmat)
ans = list(fit = muhat)
return (ans)
}
if (!is.data.frame(newData)) {
#newData = as.data.frame(newData)
stop ("newData must be a data frame!")
}
#shapes = object$shapes
#new: used for ci
prior.w = object$prior.w
vh = object$vh
y = object$y
#new: only use shapes for x != umb or tree
shapes = object$shapesx
np = object$d0; capm = object$capm; capk = object$capk; capt = object$capt; capu = object$capu
#new:
xid10 = object$xid1; xid20 = object$xid2;
uid1 = object$uid1; uid2 = object$uid2; tid1 = object$tid1; tid2 = object$tid2
#new:
xmat0 = object$xmat0; knots0 = object$knots0; numknots0 = object$numknots0; sps0 = object$sps0; ms0 = object$ms0
zmat = object$zmat; umb = object$umb; tr = object$tr
#new:
ztb = object$ztb; zid1 = object$zid1; zid2 = object$zid2; iz = 1
bigmat = object$bigmat; umbrella.delta = object$umbrella.delta; tree.delta = object$tree.delta
coefs = object$coefs; zcoefs = object$zcoefs; vcoefs = object$vcoefs; xcoefs0 = object$xcoefs; ucoefs = object$ucoefs; tcoefs = object$tcoefs
tt = object$tms
Terms = delete.response(tt)
#model.frame will re-organize newData in the original order in formula
m = model.frame(Terms, newData)
#print (m)
newdata = m
#print (head(newdata))
#new:
newx0 = NULL; newxv = NULL; newx = NULL; newx_s = NULL; newu = NULL; newt = NULL; newz = NULL; newv = NULL
#newz = list(); iz = 1;
rn = nrow(newdata)
#print (rn)
#new:
newetahat = 0; newmuhat = 0
#new: newx_sbasis
newxbasis = NULL; newx_sbasis = NULL; newubasis = NULL; newtbasis = NULL; newbigmat = NULL
#######################
#local helper function#
#######################
my_line = function(xp = NULL, y, x, end, start) {
slope = NULL
intercept = NULL
yp = NULL
slope = (y[end] - y[start]) / (x[end] - x[start])
intercept = y[end] - slope * x[end]
yp = intercept + slope * xp
ans = new.env()
ans$slope = slope
ans$intercept = intercept
ans$yp = yp
ans
}
#get shape attributes and elem out of newdata
for (i in 1:ncol(newdata)) {
if (is.null(attributes(newdata[,i])$shape)) {
if (is.factor(newdata[,i])) {
lvli = levels(newdata[,i])
ztbi = levels(ztb[[iz]])
newdatai = NULL
if (!any(lvli %in% ztbi)) {
stop ("new factor level must be among factor levels in the fit!")
} else {
id1 = which(ztbi %in% lvli)
#if (any(id1 > 1)) {
#delete the base level
klvls = length(ztbi)
if (klvls > 1) {
newimat = matrix(0, nrow = rn, ncol = klvls-1)
for (i1 in 1:rn) {
if (newdata[i1,i] != ztbi[1]) {
id_col = which(ztbi %in% newdata[i1,i]) - 1
newimat[i1,id_col] = 1
}
}
#for (i1 in 1: klvls) {
# which(levels(newdatai) == ztbi[i1])
# for (i2 in 1:rn) {
# if (levels(newdatai[i2, i1]) == lvli[i2]) {
# if (lvli[i2] != ztbi[1]) {
# newimat[i2, i1] = 1
# }
# }
# }
#}
#ztb_use = ztbi[-1]
#nc = length(id1[id1 > 1])
#newdatai = matrix(0, nrow = rn, ncol = nc)
#for (ic in 1:nc) {
# id2 = which(newdata[,i] == ztb_use[ic])
# newdatai[id2, ic] = 1
#}
newdatai = newimat
}
}
#if (length(levels(newdata[,i])) > 2) {
# klvls = length(levels(newdata[,i]))
# vals = as.numeric(levels(newdata[,i]))
# newimat = matrix(0, nrow = rn, ncol = klvls - 1)
# for (i1 in 1: (klvls - 1)) {
# for (i2 in 1:rn) {
# if (newdatai[i2] == vals[i1 + 1]) {
# newimat[i2, i1] = 1
# }
# }
# }
# newdatai = newimat
#}
} else {
newdatai = newdata[,i]
}
newz = cbind(newz, newdatai)
iz = iz + 1
#print (head(newz))
}
if (is.numeric(attributes(newdata[,i])$shape)) {
newx0 = cbind(newx0, newdata[,i])
if ((attributes(newdata[,i])$shape > 2 & attributes(newdata[,i])$shape < 5) | (attributes(newdata[,i])$shape > 10 & attributes(newdata[,i])$shape < 13)) {
newxv = cbind(newxv, newdata[,i])
}
}
if (is.character(attributes(newdata[,i])$shape)) {
if (attributes(newdata[,i])$shape == "tree") {
newt = cbind(newt, newdata[,i])
}
if (attributes(newdata[,i])$shape == "umbrella") {
newu = cbind(newu, newdata[,i])
}
}
}
#print (head(newt))
#print (head(newu))
#print (head(newx0))
#print (head(newxv))
#print (head(newz))
#print (head(newv))
#new: separate x and x_s, move shape 17 to the beginning
if (!is.null(shapes)) {
if (any(shapes == 17)) {
kshapes <- length(shapes)
obs <- 1:kshapes
idx_s <- obs[which(shapes == 17)]; idx <- obs[which(shapes != 17)]
newx1 <- newx0
shapes0 <- 1:kshapes*0
newx1[,1:length(idx_s)] <- newx0[,idx_s]
shapes0[1:length(idx_s)] <- shapes[idx_s]
if (length(idx) > 0) {
newx1[,(1 + length(idx_s)):kshapes] <- newx0[,idx]
shapes0[(1 + length(idx_s)):kshapes] <- shapes[idx]
}
newx0 <- newx1; shapes <- shapes0
}
#new:xmat is ordinal, xmat_s is smooth
xmat = xmat_s = NULL
newx = newx_s = NULL
knots = NULL
ms_x = ms = NULL
sh_x = sh = NULL
if (all(shapes < 9)) {
newx = newx0
xid1 = xid10; xid2 = xid20
xmat = xmat0
#new:
sh_x = shapes
ms_x = ms0
} else if (all(shapes > 8)) {
newx_s = newx0
xid1_s = xid10; xid2_s = xid20
xmat_s = xmat0
numknots = numknots0
knots = knots0
sps = sps0
sh = shapes
ms = ms0
} else if (any(shapes > 8) & any(shapes < 9)) {
newx = newx0[, shapes < 9, drop = FALSE]
xmat = xmat0[, shapes < 9, drop = FALSE]
#new:
sh_x = shapes[shapes < 9]
ms_x = ms0[shapes < 9]
xid1 = xid10[shapes < 9]; xid2 = xid20[shapes < 9]
newx_s = newx0[, shapes > 8, drop = FALSE]
xmat_s = xmat0[, shapes > 8, drop = FALSE]
sh = shapes[shapes > 8]
ms = ms0[shapes > 8]
xid1_s = xid10[shapes > 8]; xid2_s = xid20[shapes > 8]
numknots = numknots0[shapes > 8]
knots = knots0[shapes > 8]
sps = sps0[shapes > 8]
}
}
#if (!is.null(shapes) & any(shapes == 17)) {
if (!is.null(shapes)) {
vcoefs = vcoefs[1:(1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))]
} else {vcoefs = vcoefs[1:(1 + capk)]}
if (capk > 0) {
vcoefs_nz = vcoefs[-c(2:(1 + capk))]
} else {vcoefs_nz = vcoefs}
newv = cbind(1:rn*0 + 1, newz, newxv)
newv_nz = cbind(1:rn*0 + 1, newxv)
#print (newv_nz)
#print (vcoefs)
#etahat_v = as.vector(newv %*% vcoefs)
etahat_v_nz = as.vector(newv_nz %*% vcoefs_nz)
#print (etahat_v_nz)
#new:
etahat_s = 1:rn*0; newx_sbasis = NULL; xs_coefs = NULL; var_xs = NULL
if (!is.null(newx_s)) {
ks = ncol(newx_s)
del = NULL
for (i in 1:ks) {
xi = xmat_s[,i]
nxi = newx_s[,i]
if (any(nxi > max(xi)) | any(nxi < min(xi))) {
stop ("No extrapolation is allowed in cgam prediction!")
}
#new: scale accordingly
#nxi = (nxi - min(xi)) / (max(xi) - min(xi))
pos1 = xid1_s[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
pos2 = xid2_s[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
xs_coefs = c(xs_coefs, xcoefs0[pos1:pos2])
msi = ms[[i]]
deli_ans = makedelta(nxi, sh[i], numknots[i], knots[[i]], space = sps[i], suppre = TRUE, interp = TRUE)
deli = deli_ans$amat
if (sh[i] > 10 & sh[i] < 13) {
#x = xmat_s[,i]
xs = sort(xi)
ord = order(xi)
nx = length(xi)
obs = 1:nx
nr = nrow(deli)
nc = length(nxi)
ms2 = matrix(0, nrow = nr, ncol = nc)
for (i1 in 1:nc) {
for (i2 in 1:nr) {
ms2[i2, i1] = my_line(xp = nxi[i1], y = msi[i2, ][ord], x = xs, end = nx, start = 1)$yp
}
}
deli = deli - ms2
} else {
deli = deli - msi
}
#new:
var_xs = c(var_xs, 1:nrow(deli)*0 + i)
del = rbind(del, deli)
}
newx_sbasis = t(del)
etahat_s = as.vector(newx_sbasis %*% xs_coefs)
}
etahat_x = 1:rn*0; newxbasis = NULL; xcoefs = NULL; var_x = NULL
if (!is.null(newx)) {
kx = ncol(xmat)
del = NULL
for (i in 1:kx) {
xi = xmat[,i]
#new:
nxi = newx[,i]
if (any(nxi > max(xi)) | any(nxi < min(xi))) {
stop ("No extrapolation is allowed in cgam prediction!")
}
#new: scale accordingly
#nxi = (nxi - min(xi)) / (max(xi) - min(xi))
shi = sh_x[i]
pos1 = xid1[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
pos2 = xid2[i] - (1 + capk + sum(shapes > 2 & shapes < 5 | shapes > 10 & shapes < 13))
xcoefs = c(xcoefs, xcoefs0[pos1:pos2])
deli = pred_del(xi, shi, nxi, ms_x[[i]])
#new:
var_x = c(var_x, 1:nrow(deli)*0 + i)
del = rbind(del, deli)
}
#new:
newxbasis = t(del)
etahat_x = as.vector(newxbasis %*% xcoefs)
}
#new:
etahat_u = 1:rn*0; newuedge = NULL
if (!is.null(newu)) {
for (j in 1:ncol(umb)) {
u = umb[,j]
nu = length(u)
us = sort(u)
ord = order(u)
if (any(newu[,j] > max(u)) | any(newu[,j] < min(u))) {
stop ("no extrapolation is allowed in cgam!")
}
#u_edges = t(umbrella.fun(u))
pos1 = uid1[j]; pos2 = uid2[j]
u_edges = t(bigmat[pos1:pos2, , drop = FALSE])
nuedge = ncol(u_edges)
newuedge0 = matrix(0, nrow = rn, ncol = nuedge)
for (i in 1:nuedge) {
ue = u_edges[,i]; udist = round(diff(ue), 4); s = sign(udist)
ueord = ue[ord]; udistord = round(diff(ueord), 4); sord = sign(udistord)
obs = 1:(nu-1)
posin = obs[sord != 0] + 1
pos = unique(c(1, posin, nu))
uk = us[pos]
uek = ueord[pos]
npos = length(pos)
nd = length(newu[,j])
for (l in 1:(npos-1)) {
if (any(newu[,j] == uk[npos])) {
ids = which(newu[,j] == uk[npos])
newuedge0[ids,i] = uek[npos]
}
if (any(newu[,j] >= uk[l]) & any(newu[,j] < uk[l+1])) {
ids = which(newu[,j] >= uk[l] & newu[,j] < uk[l+1])
newuedge0[ids,i] = uek[l]
}
}
}
newuedge = cbind(newuedge, newuedge0)
}
newubasis = newuedge
etahat_u = as.vector(newubasis %*% ucoefs)
}
# we don't have a newtbasis
# estimate tree
etahat_t = 1:rn*0
#print (newt)
if (!is.null(newt)) {
#each column of tru is a "table" of all levels of a tree var
tru = unique(tr)#; placebo = min(tru)
#tr_etahat = 0
for (i in 1: ncol(newt)) {
pos1 = tid1[i] - np - capm - capu
pos2 = tid2[i] - np - capm - capu
newtu = unique(newt[,i])
#tcoefi = tcoefs[pos1:pos2]
#check! add 0 to be the coef for placebo
tcoefi = c(0, tcoefs[pos1:pos2])
if (!all(newtu %in% tru[,i])) {
stop ("new tree ordering factor must be among the old tree ordering factors!")
}
#placebo = min(tru[,i])
#new:
placebo = object$pl[i]
if (any(newtu != placebo)) {
#id_cf = which(tru[,i] %in% newtu) - 1 #no coef for placebo
#tr_etahat = tr_etahat + tcoefi[id_cf]
#id_cf = which(newtu > placebo)
id_cf = which(newtu != placebo)
t_use = newtu[id_cf]
nt = length(t_use)
for (it in 1:nt) {
etahat_t[newt[,i] == t_use[it]] = etahat_t[newt[,i] == t_use[it]] + tcoefi[tru[,i] == t_use[it]]
}
}
}
#etahat_t = tr_etahat
}
#print (etahat_v)
#print (etahat_s)
#print (etahat_x)
#print (etahat_u)
#print (etahat_t)
etahat_z = 1:rn*0; zcf = NULL
#print (newt)
#print (newdata)
if (!is.null(newz)) {
#delete the coef for 1
zcoefs = zcoefs[-1]
etahat_z = newz %*% zcoefs
#ztbu = unique(ztb)#; placebo = min(ztbu)
#lz = ncol(newz)
#for (i in 1:lz) {
# pos1 = zid1[i]
# pos2 = zid2[i]
# zcoefi = zcoefs[pos1:pos2]
# zlvi = zcoefs[pos1:pos2]
# etahat_z[newz[,i] == 1] = etahat_z[newz[,i] == 1] + zcoefi
#}
}
#print (etahat_z)
etahat_v = etahat_v_nz + etahat_z
newetahat = etahat_v + etahat_s + etahat_x + etahat_u + etahat_t
newmuhat = as.vector(muhat.fun(newetahat, fml = family$family))
if (!is.null(newt)) {
newtbasis = t(tree.delta[,1:nrow(newData),drop = FALSE])
}
#print (newv)
#newbigmat = t(cbind(newv, newx_sbasis, newxbasis, newubasis, newtbasis))
#ans = list(v = newv, xbs = newxbasis, x_sbs = newx_sbasis, ubs = newubasis, tbs = newtbasis, bigmat = newbigmat, vcoefs = vcoefs, xcoefs = xcoefs, xs_coefs = xs_coefs, ucoefs = ucoefs, etahat = newetahat, muhat = newmuhat)
#ans = list(muhat = newmuhat)
#muhat = newmuhat