Nothing
R/family.censored.R
In VGAM: Vector Generalized Linear and Additive Models
Defines functions
truncweibull
pgamma.deriv.unscaled
is.na.SurvS4
is.SurvS4
weibullR
weibull.mean
cens.rayleigh
cens.normal
ecens.exponential
cens.poisson
Documented in
cens.normal
cens.poisson
cens.poisson
cens.rayleigh
is.na.SurvS4
is.SurvS4
pgamma.deriv.unscaled
pgamma.deriv.unscaled
truncweibull
weibull.mean
weibullR
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
cens.poisson <-
function(link = "loglink", imu = NULL,
biglambda = 10,
smallno = 1e-10) {
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Censored Poisson distribution\n\n",
"Link: ", namesof("mu", link, earg = earg), "\n",
"Variance: mu"),
infos = eval(substitute(function(...) {
list(M1 = 1,
expected = FALSE,
multipleResponses = FALSE,
parameters.names = c("mu"),
link = .link ,
biglambda = .biglambda ,
smallno = .smallno )
}, list( .link = link,
.biglambda = biglambda,
.smallno = smallno ))),
initialize = eval(substitute(expression({
if (anyNA(y))
stop("NAs are not allowed in the response")
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 3,
Is.integer.y = TRUE)
centype <- attr(y, "type")
if (centype == "right") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- rep_len(FALSE, n)
extra$interval <- rep_len(FALSE, n)
init.mu <- pmax(y[, 1], 1/8)
} else
if (centype == "left") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- rep_len(FALSE, n)
extra$leftcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
init.mu <- pmax(y[, 1], 1/8)
} else
if (centype == "interval" ||
centype == "interval2") {
temp <- y[, 3]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- ifelse(temp == 2, TRUE, FALSE)
extra$intervalcensored <- ifelse(temp == 3, TRUE, FALSE)
init.mu <- pmax((y[, 1] + y[, 2])/2, 1/8) # intervalcensored
if (any(extra$uncensored))
init.mu[extra$uncensored] <-
pmax(y[extra$uncensored, 1], 1/8)
if (any(extra$rightcensored))
init.mu[extra$rightcensored] <-
pmax(y[extra$rightcensored, 1], 1/8)
if (any(extra$leftcensored))
init.mu[extra$leftcensored] <-
pmax(y[extra$leftcensored, 1], 1/8)
} else
if (centype == "counting") {
stop("type == 'counting' not compatible with cens.poisson()")
init.mu <- pmax(y[, 1], 1/8)
stop("currently not working")
} else
stop("response have to be in a class of SurvS4")
if (length( .imu )) init.mu <- 0 * y[, 1] + .imu
predictors.names <-
namesof("mu", .link, earg = .earg, short = TRUE)
if (!length(etastart))
etastart <- theta2eta(init.mu, link = .link , earg = .earg )
}), list( .link = link, .earg = earg, .imu = imu))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu <- eta2theta(eta, link = .link , earg = .earg )
as.vector(mu)
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c("mu" = .link )
misc$earg <- list("mu" = .earg )
}), list( .link = link, .earg = earg ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
loglikelihood = function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$intervalcensored
if (residuals){
stop("loglikelihood residuals not implemented yet")
} else {
sum(w[cen0] * dpois(y[cen0, 1], mu[cen0], log = TRUE)) +
sum(w[cenU] * log1p(-ppois(y[cenU, 1] - 1, mu[cenU]))) +
sum(w[cenL] * ppois(y[cenL, 1] - 1, mu[cenL], log.p = TRUE)) +
sum(w[cenI] * log(ppois(y[cenI, 2], mu[cenI]) -
ppois(y[cenI, 1], mu[cenI])))
}
},
vfamily = "cens.poisson",
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta, link = .link , earg = .earg )
okay1 <- all(is.finite(lambda)) && all(0 < lambda)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$intervalcensored
lambda <- eta2theta(eta, link = .link , earg = .earg )
dl.dlambda <- y[, 1] / lambda - 1 # uncensored
yllim <- yulim <- y[, 1] # uncensored
if (any(cenU)) {
yllim[cenU] <- y[cenU, 1]
densm1 <- dpois(yllim-1, lambda)
queue <- ppois(yllim-1, lambda, lower.tail = FALSE)
dl.dlambda[cenU] <- densm1[cenU] / queue[cenU]
if (any(fix.up <- (is.na(dl.dlambda) |
queue < .smallno ) &
cenU &
lambda > .biglambda )) {
zedd <- (yllim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
dl.dlambda[fix.up] <-
mills.ratio(-zedd) / sqrt(lambda[fix.up])
}
} # cenU
if (any(cenL)) {
yulim[cenL] <- y[cenL, 1] - 1
densm0 <- dpois(yulim, lambda)
Queue <- ppois(yulim, lambda) # Left tail probability
dl.dlambda[cenL] <- -densm0[cenL] / Queue[cenL]
if (any(fix.up <- (is.na(dl.dlambda) |
Queue < .smallno ) &
cenL &
lambda > .biglambda )) {
dl.dlambda[fix.up] <-
-mills.ratio(yulim[fix.up] / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])) / sqrt(lambda[fix.up])
}
} # cenL
if (any(cenI)) {
yllim[cenI] <- y[cenI, 1] + 1
yulim[cenI] <- y[cenI, 2]
Queue1 <- ppois(yllim-1, lambda)
Queue2 <- ppois(yulim , lambda)
densm02 <- dpois(yulim , lambda)
densm12 <- dpois(yllim-1, lambda)
dl.dlambda[cenI] <-
(-densm02[cenI]+densm12[cenI]) / (Queue2[cenI]-Queue1[cenI])
} # cenI
dlambda.deta <- dtheta.deta(theta = lambda,
link = .link , earg = .earg )
der1 <- c(w) * dl.dlambda * dlambda.deta
der1
}), list( .link = link, .earg = earg,
.biglambda = biglambda ,
.smallno = smallno ))),
weight = eval(substitute(expression({
d2lambda.deta2 <- d2theta.deta2(theta = lambda,
link = .link , earg = .earg )
d2l.dlambda2 <- 1 / lambda # uncensored; Fisher scoring
if (any(cenU)) {
densm2 <- dpois(yllim-2, lambda)
d2l.dlambda2[cenU] <- (dl.dlambda[cenU])^2 -
(densm2[cenU] - densm1[cenU]) / queue[cenU]
if (any(fix.up <- is.na(d2l.dlambda2) | fix.up)) {
warning("incorrect as one of the queue uses yulim, not yulim-1")
zeddm1 <- (yllim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
zeddm2 <- (yllim[fix.up] - 2) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
d2l.dlambda2[fix.up] <- (dl.dlambda[fix.up])^2 -
(mills.ratio(-zeddm2) -
mills.ratio(-zeddm1)) / sqrt(lambda[fix.up])
}
} # cenU
if (any(cenL)) {
densm1 <- dpois(yulim-1, lambda)
d2l.dlambda2[cenL] <- (dl.dlambda[cenL])^2 -
(densm0[cenL] - densm1[cenL]) / Queue[cenL]
if (any(fix.up <- is.na(d2l.dlambda2) | fix.up)) {
d2l.dlambda2[fix.up] <- (dl.dlambda[fix.up])^2 -
(mills.ratio(yulim[fix.up] / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])) -
mills.ratio((yulim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up]))) / sqrt(lambda[fix.up])
}
} # cenL
if (any(cenI)) {
densm03 <- dpois(yulim-1, lambda)
densm13 <- dpois(yllim-2, lambda)
d2l.dlambda2[cenI] <- (dl.dlambda[cenI])^2 -
(densm13[cenI]-densm12[cenI]-densm03[cenI] +
densm02[cenI]) / (Queue2[cenI]-Queue1[cenI])
} # cenI
wz <- c(w) * (dlambda.deta^2) * d2l.dlambda2
wz
}), list( .link = link, .earg = earg,
.biglambda = biglambda ,
.smallno = smallno ))))
} # cens.poisson
if (FALSE)
cens.exponential <-
ecens.exponential <- function(link = "loglink", location = 0) {
if (!is.Numeric(location, length.arg = 1))
stop("bad input for 'location'")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Censored exponential distribution\n\n",
"Link: ", namesof("rate", link, tag = TRUE), "\n",
"Mean: ", "mu = ", location, " + 1 / ",
namesof("rate", link, tag = FALSE), "\n",
"Variance: ",
if (location == 0) "Exponential: mu^2" else
paste("(mu-", location, ")^2", sep = "")),
initialize = eval(substitute(expression({
extra$location <- .location
if (any(y[, 1] <= extra$location))
stop("all responses must be greater than ", extra$location)
predictors.names <- namesof("rate", .link , .earg , tag = FALSE)
type <- attr(y, "type")
if (type == "right" || type == "left"){
mu <- y[, 1] + (abs(y[, 1] - extra$location) < 0.001) / 8
} else if (type == "interval") {
temp <- y[, 3]
mu <- ifelse(temp == 3,
y[, 2] + (abs(y[, 2] - extra$location) < 0.001) / 8,
y[, 1] + (abs(y[, 1] - extra$location) < 0.001) / 8)
}
if (!length(etastart))
etastart <- theta2eta(1/(mu-extra$location), .link , .earg )
if (type == "right") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- rep_len(FALSE, n)
extra$interval <- rep_len(FALSE, n)
} else
if (type == "left") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- rep_len(FALSE, n)
extra$leftcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
} else
if (type == "counting") {
stop("type == 'counting' not recognized")
extra$uncensored <- rep(temp == 1, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
extra$leftcensored <- rep_len(FALSE, n)
extra$rightcensored <- rep_len(FALSE, n)
extra$counting <- ifelse(temp == 0, TRUE, FALSE)
} else
if (type == "interval") {
temp <- y[, 3]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- ifelse(temp == 2, TRUE, FALSE)
extra$interval <- ifelse(temp == 3, TRUE, FALSE)
} else
stop("'type' not recognized")
}), list( .location = location, .link = link ))),
linkinv = eval(substitute(function(eta, extra = NULL)
extra$location + 1 / eta2theta(eta, .link , .earg ),
list( .link = link ) )),
last = eval(substitute(expression({
misc$location <- extra$location
misc$link <- c("rate" = .link)
misc$multipleResponses <- FALSE
}), list( .link = link ))),
link = eval(substitute(function(mu, extra = NULL)
theta2eta(1 / (mu - extra$location), .link , .earg ),
list( .link = link ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
rate <- 1 / (mu - extra$location)
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$interval
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cenL] * log1p(-exp(-rate[cenL] *
(y[cenL, 1] - extra$location)))) +
sum(w[cenU] * (-rate[cenU]*(y[cenU, 1]-extra$location))) +
sum(w[cen0] * (log(rate[cen0]) -
rate[cen0]*(y[cen0, 1]-extra$location))) +
sum(w[cenI] * log(-exp(-rate[cenI]*(y[cenI, 2]-extra$location))+
exp(-rate[cenI]*(y[cenI, 1]-extra$location))))
}, list( .link = link ))),
vfamily = c("ecens.exponential"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
rate <- 1 / (mu - extra$location)
okay1 <- all(is.finite(rate)) && all(0 < rate)
okay1
}, list( .link = link ))),
deriv = eval(substitute(expression({
rate <- 1 / (mu - extra$location)
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$interval
dl.drate <- 1/rate - (y[, 1]-extra$location) # uncensored
tmp200 <- exp(-rate*(y[, 1]-extra$location))
tmp200b <- exp(-rate*(y[, 2]-extra$location)) # interval censored
if (any(cenL))
dl.drate[cenL] <- (y[cenL, 1]-extra$location) *
tmp200[cenL] / (1 - tmp200[cenL])
if (any(cenU))
dl.drate[cenU] <- -(y[cenU, 1]-extra$location)
if (any(cenI))
dl.drate[cenI] <- ((y[cenI, 2] - extra$location) *
tmp200b[cenI] - (y[cenI, 1] - extra$location) *
tmp200[cenI]) / (-tmp200b[cenI] + tmp200[cenI])
drate.deta <- dtheta.deta(rate, .link , .earg )
c(w) * dl.drate * drate.deta
}), list( .link = link ) )),
weight = eval(substitute(expression({
A123 <- ((mu-extra$location)^2) # uncensored d2l.drate2
Lowpt <- ifelse(cenL, y[, 1], extra$location)
Lowpt <- ifelse(cenI, y[, 1], Lowpt) #interval censored
Upppt <- ifelse(cenU, y[, 1], Inf)
Upppt <- ifelse(cenI, y[, 2], Upppt) #interval censored
tmp300 <- exp(-rate*(Lowpt - extra$location))
d2l.drate2 <- 0 * y[, 1]
ind50 <- Lowpt > extra$location
d2l.drate2[ind50] <- (Lowpt[ind50]-extra$location)^2 *
tmp300[ind50] / (1-tmp300[ind50])
d2l.drate2 <- d2l.drate2 + (exp(-rate*(Lowpt-extra$location)) -
exp(-rate*(Upppt-extra$location))) * A123
wz <- c(w) * (drate.deta^2) * d2l.drate2
wz
}), list( .link = link ))))
} # cens.exponential
cennormal <-
cens.normal <- function(lmu = "identitylink", lsd = "loglink",
imethod = 1, zero = "sd") {
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsd <- as.list(substitute(lsd))
esd <- link2list(lsd)
lsd <- attr(esd, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
new("vglmff",
blurb = c("Censored univariate normal\n\n",
"Links: ", namesof("mu", lmu, tag = TRUE), "; ",
namesof("sd", lsd, tag = TRUE), "\n",
"Conditional variance: sd^2"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
zero = .zero ,
multiple.responses = FALSE,
parameters.names = c("mu", "sd"),
expected = TRUE )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
if (!length(extra$leftcensored))
extra$leftcensored <- rep_len(FALSE, n)
if (!length(extra$rightcensored))
extra$rightcensored <- rep_len(FALSE, n)
if (any(extra$rightcensored & extra$leftcensored))
stop("some observations are both right and left censored!")
predictors.names <-
c(namesof("mu", .lmu , earg = .emu , tag = FALSE),
namesof("sd", .lsd , earg = .esd , tag = FALSE))
if (!length(etastart)) {
anyc <- extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # can be all data
junk <- lm.wfit(x = cbind(x[!i11, ]),
y = y[!i11], w = w[!i11])
sd.y.est <- sqrt(sum(w[!i11] * junk$resid^2) / junk$df.residual)
etastart <- cbind(mu = y,
rep_len(theta2eta(sd.y.est, .lsd), n))
if (any(anyc))
etastart[anyc, 1] <- x[anyc, , drop = FALSE] %*% junk$coeff
}
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd,
.imethod = imethod ))),
linkinv = eval(substitute( function(eta, extra = NULL) {
eta2theta(eta[, 1], .lmu , earg = .emu )
}, list( .lmu = lmu, .emu = emu ))),
last = eval(substitute(expression({
misc$link <- c("mu" = .lmu , "sd" = .lsd )
misc$earg <- list("mu" = .emu , "sd" = .esd )
misc$expected <- TRUE
misc$multipleResponses <- FALSE
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cen0 <- !cenL & !cenU # uncensored obsns
mum <- eta2theta(eta[, 1], .lmu , earg = .emu )
sdv <- eta2theta(eta[, 2], .lsd , earg = .esd )
Lower <- ifelse(cenL, y, -Inf)
Upper <- ifelse(cenU, y, Inf)
ell1 <- -log(sdv[cen0]) -
0.5 * ((y[cen0] - mum[cen0])/sdv[cen0])^2
ell2 <- log1p(-pnorm((mum[cenL] - Lower[cenL]) / sdv[cenL]))
ell3 <- log1p(-pnorm(( Upper[cenU] - mum[cenU]) / sdv[cenU]))
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cen0] * ell1) +
sum(w[cenL] * ell2) +
sum(w[cenU] * ell3)
}, list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
vfamily = c("cens.normal"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mum <- eta2theta(eta[, 1], .lmu )
sdv <- eta2theta(eta[, 2], .lsd )
okay1 <- all(is.finite(mum)) &&
all(is.finite(sdv)) && all(0 < sdv)
okay1
}, list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
deriv = eval(substitute(expression({
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cen0 <- !cenL & !cenU # uncensored obsns
Lower <- ifelse(cenL, y, -Inf)
Upper <- ifelse(cenU, y, Inf)
mum <- eta2theta(eta[, 1], .lmu )
sdv <- eta2theta(eta[, 2], .lsd )
dl.dmu <- (y-mum) / sdv^2
dl.dsd <- (((y-mum)/sdv)^2 - 1) / sdv
dmu.deta <- dtheta.deta(mum, .lmu , earg = .emu )
dsd.deta <- dtheta.deta(sdv, .lsd , earg = .esd )
if (any(cenL)) {
mumL <- mum - Lower
temp21L <- mumL[cenL] / sdv[cenL]
PhiL <- pnorm(temp21L)
phiL <- dnorm(temp21L)
fred21 <- phiL / (1 - PhiL)
dl.dmu[cenL] <- -fred21 / sdv[cenL]
dl.dsd[cenL] <- mumL[cenL] * fred21 / sdv[cenL]^2
rm(fred21)
}
if (any(cenU)) {
mumU <- Upper - mum
temp21U <- mumU[cenU] / sdv[cenU]
PhiU <- pnorm(temp21U)
phiU <- dnorm(temp21U)
fred21 <- phiU / (1 - PhiU)
dl.dmu[cenU] <- fred21 / sdv[cenU] # Negated
dl.dsd[cenU] <- mumU[cenU] * fred21 / sdv[cenU]^2
rm(fred21)
}
c(w) * cbind(dl.dmu * dmu.deta,
dl.dsd * dsd.deta)
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
weight = eval(substitute(expression({
A1 <- 1 - pnorm((mum - Lower) / sdv) # Lower
A3 <- 1 - pnorm((Upper - mum) / sdv) # Upper
A2 <- 1 - A1 - A3 # Middle; uncensored
wz <- matrix(0, n, 3)
wz[, iam(1, 1,M)] <- A2 * 1 / sdv^2 # ed2l.dmu2
wz[, iam(2, 2,M)] <- A2 * 2 / sdv^2 # ed2l.dsd2
mumL <- mum - Lower
temp21L <- mumL / sdv
PhiL <- pnorm(temp21L)
phiL <- dnorm(temp21L)
temp31L <- ((1-PhiL) * sdv)^2
wz.cenL11 <- phiL * (phiL - (1-PhiL)*temp21L) / temp31L
wz.cenL22 <- mumL * phiL * ((1-PhiL) * (2 - temp21L^2) +
mumL * phiL / sdv) / (sdv * temp31L)
wz.cenL12 <- phiL * ((1-PhiL)*(temp21L^2 - 1) -
temp21L*phiL) / temp31L
wz.cenL11[!is.finite(wz.cenL11)] <- 0
wz.cenL22[!is.finite(wz.cenL22)] <- 0
wz.cenL12[!is.finite(wz.cenL12)] <- 0
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + A1 * wz.cenL11
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] + A1 * wz.cenL22
wz[, iam(1, 2, M)] <- A1 * wz.cenL12
mumU <- Upper - mum # often Inf
temp21U <- mumU / sdv # often Inf
PhiU <- pnorm(temp21U) # often 1
phiU <- dnorm(temp21U) # often 0
temp31U <- ((1-PhiU) * sdv)^2 # often 0
tmp8 <- (1-PhiU)*temp21U
wzcenU11 <- phiU * (phiU - tmp8) / temp31U
tmp9 <- (1-PhiU) * (2 - temp21U^2)
wzcenU22 <- mumU * phiU * (tmp9 +
mumU * phiU / sdv) / (sdv * temp31U)
wzcenU12 <- -phiU * ((1-PhiU)*(temp21U^2 - 1) -
temp21U*phiU) / temp31U
wzcenU11[!is.finite(wzcenU11)] <- 0 # Needed when Upper==Inf
wzcenU22[!is.finite(wzcenU22)] <- 0 # Needed when Upper==Inf
wzcenU12[!is.finite(wzcenU12)] <- 0 # Needed when Upper==Inf
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + A3 * wzcenU11
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] + A3 * wzcenU22
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] + A3 * wzcenU12
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dmu.deta^2
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] * dsd.deta^2
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] * dmu.deta * dsd.deta
c(w) * wz
}), list( .lmu = lmu, .lsd = lsd ))))
} # cens.normal
cens.rayleigh <- function(lscale = "loglink",
oim = TRUE) {
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.logical(oim) || length(oim) != 1)
stop("bad input for argument 'oim'")
new("vglmff",
blurb = c("Censored Rayleigh distribution\n\n",
"f(y) = y*exp(-0.5*(y/scale)^2)/scale^2, y>0, scale>0\n",
"Link: ",
namesof("scale", lscale, earg = escale ), "\n", "\n",
"Mean: scale * sqrt(pi / 2)"),
initialize = eval(substitute(expression({
if (NCOL(y) != 1)
stop("response must be a vector or a one-column matrix")
if (length(extra$leftcensored))
stop("cannot handle left-censored data")
if (!length(extra$rightcensored))
extra$rightcensored <- rep_len(FALSE, n)
predictors.names <-
namesof("scale", .lscale , earg = .escale , tag = FALSE)
if (!length(etastart)) {
a.init <- (y+1/8) / sqrt(pi/2)
etastart <- theta2eta(a.init, .lscale , earg = .escale )
}
}), list( .lscale = lscale, .escale = escale ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
Scale * sqrt(pi/2)
}, list( .lscale = lscale, .escale = escale ))),
last = eval(substitute(expression({
misc$link <- c("scale" = .lscale )
misc$earg <- list("scale" = .escale )
misc$oim <- .oim
}), list( .lscale = lscale, .escale = escale,
.oim = oim ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
cen0 <- !extra$rightcensored # uncensored obsns
cenU <- extra$rightcensored
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cen0] * (log(y[cen0]) - 2*log(Scale[cen0]) -
0.5*(y[cen0]/Scale[cen0])^2)) -
sum(w[cenU] * (y[cenU]/Scale[cenU])^2) * 0.5
}, list( .lscale = lscale, .escale = escale ))),
vfamily = c("cens.rayleigh"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
okay1 <- all(is.finite(Scale)) && all(0 < Scale)
okay1
}, list( .lscale = lscale, .escale = escale ))),
deriv = eval(substitute(expression({
cen0 <- !extra$rightcensored # uncensored obsns
cenU <- extra$rightcensored
Scale <- eta2theta(eta, .lscale , earg = .escale )
dl.dScale <- ((y/Scale)^2 - 2) / Scale
dScale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
dl.dScale[cenU] <- y[cenU]^2 / Scale[cenU]^3
c(w) * dl.dScale * dScale.deta
}), list( .lscale = lscale,
.escale = escale ))),
weight = eval(substitute(expression({
ned2l.dScale2 <- 4 / Scale^2
wz <- dScale.deta^2 * ned2l.dScale2
if ( .oim ) {
d2l.dScale2 <- 3 * (y[cenU])^2 / (Scale[cenU])^4
d2Scale.deta2 <- d2theta.deta2(Scale[cenU], .lscale ,
earg = .escale )
wz[cenU] <- (dScale.deta[cenU])^2 * d2l.dScale2 -
dl.dScale[cenU] * d2Scale.deta2
} else {
ned2l.dScale2[cenU] <- 6 / (Scale[cenU])^2
wz[cenU] <- (dScale.deta[cenU])^2 * ned2l.dScale2[cenU]
}
c(w) * wz
}), list( .lscale = lscale, .escale = escale,
.oim = oim ))))
} # cens.rayleigh
weibull.mean <-
function(lmean = "loglink", lshape = "loglink",
imean = NULL, ishape = NULL,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1, zero = "shape") {
imeann <- imean
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
lmeann <- as.list(substitute(lmean))
emeann <- link2list(lmeann)
lmeann <- attr(emeann, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (length(ishape))
if (!is.Numeric(ishape, positive = TRUE))
stop("argument 'ishape' values must be positive")
if (length(imeann))
if (!is.Numeric(imeann, positive = TRUE))
stop("argument 'imean' values must be positive")
blurb.vec <- c(namesof("mean", lmeann, earg = emeann),
namesof("shape", lshape, earg = eshape))
new("vglmff",
blurb = c("Weibull distribution (parameterized by the mean)\n\n",
"Links: ",
blurb.vec[1], ", ",
blurb.vec[2], "\n",
"Mean: mean\n",
"Variance: mean^2 * (gamma(1 + 2/shape) / ",
"gamma(1 + 1/shape)^2 - 1)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero,
.lmeann = lmeann ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("mean", "shape"),
lmean = .lmeann ,
lshape = .lshape ,
zero = .zero )
}, list( .zero = zero,
.lmeann = lmeann, .lshape = lshape ))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
mynames1 <- param.names("mean" , ncoly, skip1 = TRUE)
mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lmeann , earg = .emeann , tag = FALSE),
namesof(mynames2, .lshape , earg = .eshape , tag = FALSE))
predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)]
Meann.init <- matrix(if (length( .imeann ))
.imeann else 0.5 * colMeans(y),
n, ncoly, byrow = TRUE) + 0.5 * y
Shape.init <- matrix(if (length( .ishape )) .ishape else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .ishape ) ||
!length( .imeann )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
if (!is.Numeric(Shape.init[, ilocal]))
Shape.init[, ilocal] <- 1 / fit0$coef["X"]
} # ilocal
etastart <-
cbind(theta2eta(Meann.init, .lmeann , earg = .emeann ),
theta2eta(Shape.init, .lshape , earg = .eshape ))[,
interleave.VGAM(M, M1 = M1)]
}
}
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape,
.imeann = imeann, .ishape = ishape,
.probs.y = probs.y,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Meann
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
last = eval(substitute(expression({
regnotok <- any(Shape <= 2)
if (any(Shape <= 1)) {
warning("MLE regularity conditions are violated",
"(shape <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < Shape & Shape < 2)) {
warning("MLE regularity conditions are violated",
"(1 < shape < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == Shape)) {
warning("MLE regularity conditions are violated",
"(shape == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than ",
"when shape > 2")
}
M1 <- extra$M1
avector <- c(rep_len( .lmeann , ncoly),
rep_len( .lshape , ncoly))
misc$link <- avector[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .emeann
misc$earg[[M1*ii ]] <- .eshape
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$RegCondOK <- !regnotok # Save this for later
misc$expected <- TRUE # all(cen0)
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape,
.imethod = imethod ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * dweibull(x = y, shape = Shape,
scale = Meann / gamma(1 + 1/Shape),
log = TRUE))
}
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
vfamily = c("weibull.mean"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
okay1 <- all(is.finite(Meann)) && all(0 < Meann) &&
all(is.finite(Shape)) && all(0 < Shape)
okay1
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
deriv = eval(substitute(expression({
M1 <- 2
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
if (FALSE) {
} else {
EulerM <- -digamma(1.0)
AA <- (EulerM - 1)^2 + (pi^2) / 6
BB <- digamma(1 + 1/Shape)
CC <- y * gamma(1 + 1/Shape) / Meann
dl.dmeann <- (CC^Shape - 1) * Shape / Meann # Agrees
dl.dshape <- 1/Shape -
(log(y/Meann) + lgamma(1 + 1/Shape)) * (CC^Shape - 1) +
(BB / Shape) * (CC^Shape - 1)
}
dmeann.deta <- dtheta.deta(Meann, .lmeann , earg = .emeann )
dshape.deta <- dtheta.deta(Shape, .lshape , earg = .eshape )
myderiv <- c(w) * cbind(dl.dmeann * dmeann.deta,
dl.dshape * dshape.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
weight = eval(substitute(expression({
if (FALSE) {
} else {
ned2l.dmeann <- (Shape / Meann)^2 #
ned2l.dshape <- AA / Shape^2 # Unchanged
ned2l.dshapemeann <- (EulerM - 1 + BB) / Meann
}
wz <- array(c(c(w) * ned2l.dmeann * dmeann.deta^2,
c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dshapemeann * dmeann.deta * dshape.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .eshape = eshape ))))
} # weibull.mean
weibullR <-
function(lscale = "loglink", lshape = "loglink",
iscale = NULL, ishape = NULL,
lss = TRUE,
nrfs = 1,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1, zero = "shape") {
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (!is.Numeric(nrfs, length.arg = 1) ||
nrfs < 0 ||
nrfs > 1)
stop("bad input for argument 'nrfs'")
if (length(ishape))
if (!is.Numeric(ishape, positive = TRUE))
stop("argument 'ishape' values must be positive")
if (length(iscale))
if (!is.Numeric(iscale, positive = TRUE))
stop("argument 'iscale' values must be positive")
scale.TF <- if (lss) c(TRUE, FALSE) else c(FALSE, TRUE)
scale.12 <- if (lss) 1:2 else 2:1
blurb.vec <- c(namesof("scale", lscale, earg = escale),
namesof("shape", lshape, earg = eshape))
blurb.vec <- blurb.vec[scale.12]
new("vglmff",
blurb = c("Weibull distribution\n\n",
"Links: ",
blurb.vec[1], ", ",
blurb.vec[2], "\n",
"Mean: scale * gamma(1 + 1/shape)\n",
"Variance: scale^2 * (gamma(1 + 2/shape) - ",
"gamma(1 + 1/shape)^2)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = if ( .lss )
c("scale", "shape") else
c("shape", "scale"),
lss = .lss ,
lscale = .lscale ,
lshape = .lshape ,
zero = .zero )
}, list( .zero = zero, .scale.12 = scale.12, .scale.TF = scale.TF,
.lscale = lscale ,
.lshape = lshape ,
.lss = lss
))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
if ( .lss ) {
mynames1 <- param.names("scale", ncoly, skip1 = TRUE)
mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lscale , earg = .escale , tag = FALSE),
namesof(mynames2, .lshape , earg = .eshape , tag = FALSE))
} else {
mynames1 <- param.names("shape", ncoly, skip1 = TRUE)
mynames2 <- param.names("scale", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lshape , earg = .eshape , tag = FALSE),
namesof(mynames2, .lscale , earg = .escale , tag = FALSE))
}
predictors.names <- predictors.names[
interleave.VGAM(M, M1 = M1)]
Shape.init <- matrix(if (length( .ishape )) .ishape else 0 + NA,
n, ncoly, byrow = TRUE)
Scale.init <- matrix(if (length( .iscale )) .iscale else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .ishape ) ||
!length( .iscale )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
if (!is.Numeric(Shape.init[, ilocal]))
Shape.init[, ilocal] <- 1 / fit0$coef["X"]
if (!is.Numeric(Scale.init[, ilocal]))
Scale.init[, ilocal] <- exp(fit0$coef["Intercept"])
} # ilocal
etastart <- if ( .lss )
cbind(theta2eta(Scale.init, .lscale , earg = .escale ),
theta2eta(Shape.init, .lshape , earg = .eshape ))[,
interleave.VGAM(M, M1 = M1)] else
cbind(theta2eta(Shape.init, .lshape , earg = .eshape ),
theta2eta(Scale.init, .lscale , earg = .escale ))[,
interleave.VGAM(M, M1 = M1)]
}
}
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.iscale = iscale, .ishape = ishape,
.probs.y = probs.y,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
Scale * gamma(1 + 1 / Shape)
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
last = eval(substitute(expression({
regnotok <- any(Shape <= 2)
if (any(Shape <= 1)) {
warning("MLE regularity conditions are violated",
"(shape <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < Shape & Shape < 2)) {
warning("MLE regularity conditions are violated",
"(1 < shape < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == Shape)) {
warning("MLE regularity conditions are violated",
"(shape == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than when ",
"shape > 2")
}
M1 <- extra$M1
avector <- if ( .lss ) c(rep_len( .lscale , ncoly),
rep_len( .lshape , ncoly)) else
c(rep_len( .lshape , ncoly),
rep_len( .lscale , ncoly))
misc$link <- avector[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- if ( .lss ) .escale else .eshape
misc$earg[[M1*ii ]] <- if ( .lss ) .eshape else .escale
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$nrfs <- .nrfs
misc$RegCondOK <- !regnotok # Save this for later
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.imethod = imethod,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss,
.nrfs = nrfs ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * dweibull(y, shape = Shape, scale = Scale, log = TRUE))
}
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
vfamily = c("weibullR"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
all(is.finite(Shape)) && all(0 < Shape)
okay1
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
deriv = eval(substitute(expression({
M1 <- 2
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
dl.dshape <- 1 / Shape + log(y / Scale) -
log(y / Scale) * (y / Scale)^Shape
dl.dscale <- (Shape / Scale) * (-1.0 + (y / Scale)^Shape)
dshape.deta <- dtheta.deta(Shape, .lshape, earg = .eshape )
dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
myderiv <- if ( .lss )
c(w) * cbind(dl.dscale * dscale.deta,
dl.dshape * dshape.deta) else
c(w) * cbind(dl.dshape * dshape.deta,
dl.dscale * dscale.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
weight = eval(substitute(expression({
EulerM <- -digamma(1.0)
ned2l.dscale <- (Shape / Scale)^2
ned2l.dshape <- (6*(EulerM - 1)^2 + pi^2)/(6*Shape^2) # KK (2003)
ned2l.dshapescale <- (EulerM-1) / Scale
wz <- if ( .lss )
array(c(c(w) * ned2l.dscale * dscale.deta^2,
c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dshapescale * dscale.deta * dshape.deta),
dim = c(n, M / M1, 3)) else
array(c(c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dscale * dscale.deta^2,
c(w) * ned2l.dshapescale * dscale.deta * dshape.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .eshape = eshape, .nrfs = nrfs,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ))))
} # weibullR
setOldClass(c("SurvS4", "Surv"))
SurvS4 <-
function (time, time2, event, type = c("right", "left", "interval",
"counting", "interval2"), origin = 0) {
nn <- length(time)
ng <- nargs()
if (missing(type)) {
if (ng == 1 || ng == 2)
type <- "right" else if (ng == 3)
type <- "counting" else stop("Invalid number of arguments")
} else {
type <- match.arg(type)
ng <- ng - 1
if (ng != 3 && (type == "interval" || type == "counting"))
stop("Wrong number of args for this type of survival data")
if (ng != 2 && (type == "right" || type == "left" ||
type == "interval2"))
stop("Wrong number of args for this type of survival data")
}
who <- !is.na(time)
if (ng == 1) {
if (!is.numeric(time))
stop("Time variable is not numeric")
ss <- cbind(time, 1)
dimnames(ss) <- list(NULL, c("time", "status"))
} else if (type == "right" || type == "left") {
if (!is.numeric(time))
stop("Time variable is not numeric")
if (length(time2) != nn)
stop("Time and status are different lengths")
if (is.logical(time2))
status <- 1 * time2 else if (is.numeric(time2)) {
who2 <- !is.na(time2)
if (max(time2[who2]) == 2)
status <- time2 - 1 else status <- time2
if (any(status[who2] != 0 & status[who2] != 1))
stop("Invalid status value")
} else stop("Invalid status value")
ss <- cbind(time, status)
dimnames(ss) <- list(NULL, c("time", "status"))
} else if (type == "counting") {
if (length(time2) != nn)
stop("Start and stop are different lengths")
if (length(event) != nn)
stop("Start and event are different lengths")
if (!is.numeric(time))
stop("Start time is not numeric")
if (!is.numeric(time2))
stop("Stop time is not numeric")
who3 <- who & !is.na(time2)
if (any(time[who3] >= time2[who3]))
stop("Stop time must be > start time")
if (is.logical(event))
status <- 1 * event else if (is.numeric(event)) {
who2 <- !is.na(event)
if (max(event[who2]) == 2)
status <- event - 1 else status <- event
if (any(status[who2] != 0 & status[who2] != 1))
stop("Invalid status value")
} else stop("Invalid status value")
ss <- cbind(time - origin, time2 - origin, status)
} else {
if (type == "interval2") {
event <- ifelse(is.na(time), 2, ifelse(is.na(time2),
0, ifelse(time == time2, 1, 3)))
if (any(time[event == 3] > time2[event == 3]))
stop("Invalid interval: start > stop")
time <- ifelse(event != 2, time, time2)
type <- "interval"
} else {
temp <- event[!is.na(event)]
if (!is.numeric(temp))
stop("Status indicator must be numeric")
if (length(temp) > 0 && any(temp != floor(temp) |
temp < 0 | temp > 3))
stop("Status indicator must be 0, 1, 2 or 3")
}
status <- event
ss <- cbind(time, ifelse(!is.na(event) & event == 3,
time2, 1), status)
}
attr(ss, "type") <- type
class(ss) <- "SurvS4"
ss
} # SurvS4
is.SurvS4 <- function(x) inherits(x, "SurvS4")
setIs(class1 = "SurvS4",
class2 = "matrix") # Forces vglm()@y to be a matrix
as.character.SurvS4 <- function (x, ...) {
class(x) <- NULL
type <- attr(x, "type")
if (type == "right") {
temp <- x[, 2]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "+", " "))
paste(format(x[, 1]), temp, sep = "")
} else if (type == "counting") {
temp <- x[, 3]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "+", " "))
paste("(", format(x[, 1]), ",", format(x[, 2]), temp, "]", sep = "")
} else if (type == "left") {
temp <- x[, 2]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "<", " "))
paste(temp, format(x[, 1]), sep = "")
} else {
stat <- x[, 3]
temp <- c("+", "", "-", "]")[stat + 1]
temp2 <- ifelse(stat == 3, paste("(", format(x[, 1]),
", ", format(x[, 2]), sep = ""), format(x[, 1]))
ifelse(is.na(stat), as.character(NA), paste(temp2, temp, sep = ""))
}
} # as.character.SurvS4
"[.SurvS4" <- function(x, i, j, drop = FALSE) {
if (missing(j)) {
temp <- class(x)
type <- attr(x, "type")
class(x) <- NULL
x <- x[i, , drop = FALSE]
class(x) <- temp
attr(x, "type") <- type
x
} else {
class(x) <- NULL
NextMethod("[")
}
} # "[.SurvS4"
is.na.SurvS4 <- function(x) {
as.vector( (1* is.na(unclass(x)))%*% rep(1, ncol(x)) >0)
}
show.SurvS4 <- function (object)
print.default(as.character.SurvS4(object), quote = FALSE)
setMethod("show", "SurvS4",
function(object)
show.SurvS4(object))
pgamma.deriv.unscaled <- function(q, shape) {
gam0 <- exp(lgamma(shape) + pgamma(q = q, shape = shape,
log.p = TRUE))
I.sq <- pgamma(q = q, shape = shape)
alld <- pgamma.deriv(q = q, shape = shape) # 6-coln matrix
tmp3 <- alld[, 3] / I.sq # RHS of eqn (4.5) of \cite{wing:1989}
G1s <- digamma(shape) + tmp3 # eqn (4.9)
gam1 <- gam0 * G1s
dG1s <- trigamma(shape) + alld[, 4] / I.sq - tmp3^2 # eqn (4.13)
G2s <- dG1s + G1s^2 # eqn (4.12)
gam2 <- gam0 * G2s
cbind("0" = gam0,
"1" = gam1,
"2" = gam2)
} # pgamma.deriv.unscaled
truncweibull <-
function(lower.limit = 1e-5,
lAlpha = "loglink", lBetaa = "loglink",
iAlpha = NULL, iBetaa = NULL,
nrfs = 1,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1,
zero = "Betaa") {
lAlpha <- as.list(substitute(lAlpha))
eAlpha <- link2list(lAlpha)
lAlpha <- attr(eAlpha, "function.name")
lBetaa <- as.list(substitute(lBetaa))
eBetaa <- link2list(lBetaa)
lBetaa <- attr(eBetaa, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (!is.Numeric(nrfs, length.arg = 1) ||
nrfs < 0 ||
nrfs > 1)
stop("bad input for argument 'nrfs'")
if (length(iAlpha))
if (!is.Numeric(iAlpha, positive = TRUE))
stop("argument 'iAlpha' values must be positive")
if (length(iBetaa))
if (!is.Numeric(iBetaa, positive = TRUE))
stop("argument 'iBetaa' values must be positive")
new("vglmff",
blurb = c("Truncated weibull distribution\n\n",
"Links: ",
namesof("Alpha", lAlpha, earg = eAlpha), ", ",
namesof("Betaa", lBetaa, earg = eBetaa), "\n",
if (length( lower.limit ) < 5)
paste("Truncation point(s): ",
lower.limit, sep = ", ") else
""),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("Alpha", "Betaa"),
lower.limit = .lower.limit ,
lAlpha = .lAlpha ,
lBetaa = .lBetaa ,
zero = .zero )
}, list( .zero = zero,
.lAlpha = lAlpha ,
.lBetaa = lBetaa ,
.lower.limit = lower.limit
))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
extra$lower.limit <- matrix( .lower.limit , n, ncoly, byrow = TRUE)
if (any(y < extra$lower.limit)) {
stop("some response values less than argument 'lower.limit'")
}
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
mynames1 <- param.names("Alpha", ncoly, skip1 = TRUE)
mynames2 <- param.names("Betaa", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lAlpha , earg = .eAlpha , tag = FALSE),
namesof(mynames2, .lBetaa , earg = .eBetaa , tag = FALSE))[
interleave.VGAM(M, M1 = M1)]
Alpha.init <- matrix(if (length( .iAlpha )) .iAlpha else 0 + NA,
n, ncoly, byrow = TRUE)
Betaa.init <- matrix(if (length( .iBetaa )) .iBetaa else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .iAlpha ) ||
!length( .iBetaa )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
aaa.init <- 1 / fit0$coef["X"]
bbb.init <- exp(fit0$coef["Intercept"])
if (!is.Numeric(Betaa.init[, ilocal]))
Betaa.init[, ilocal] <- aaa.init
if (!is.Numeric(Alpha.init[, ilocal]))
Alpha.init[, ilocal] <- (1 / bbb.init)^aaa.init
} # ilocal
} else {
Alpha.init <- rep_len( .iAlpha , n)
Betaa.init <- rep_len( .iBetaa , n)
}
etastart <-
cbind(theta2eta(Alpha.init, .lAlpha , earg = .eAlpha ),
theta2eta(Betaa.init, .lBetaa , earg = .eBetaa ))[,
interleave.VGAM(M, M1 = M1)]
}
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.iBetaa = iBetaa, .iAlpha = iAlpha,
.lower.limit = lower.limit,
.probs.y = probs.y,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
aTb <- Alpha * extra$lower.limit^Betaa
wingo3 <- pgamma.deriv.unscaled(q = aTb, shape = 1 + 1 / Betaa)
exp.aTb <- exp(aTb)
(gamma(1 + 1 / Betaa) - wingo3[, 1]) *
exp.aTb / Alpha^(1 / Betaa)
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit) )),
last = eval(substitute(expression({
aaa.hat <- Betaa
regnotok <- any(aaa.hat <= 2)
if (any(aaa.hat <= 1)) {
warning("MLE regularity conditions are violated",
"(Betaa <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < aaa.hat & aaa.hat < 2)) {
warning("MLE regularity conditions are violated",
"(1 < Betaa < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == aaa.hat)) {
warning("MLE regularity conditions are violated",
"(Betaa == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than when ",
"Betaa > 2")
}
M1 <- extra$M1
misc$link <-
c(rep_len( .lAlpha , ncoly),
rep_len( .lBetaa , ncoly))[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .eAlpha
misc$earg[[M1*ii ]] <- .eBetaa
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$nrfs <- .nrfs
misc$RegCondOK <- !regnotok # Save this for later
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.imethod = imethod,
.lower.limit = lower.limit,
.nrfs = nrfs ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
Shape <- Betaa
Scale <- 1 / Alpha^(1/Betaa)
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * (dweibull(x = y, shape = Shape,
scale = Scale, log = TRUE) -
pweibull(q = extra$lower.limit, shape = Shape,
scale = Scale, log.p = TRUE,
lower.tail = FALSE)))
}
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
vfamily = c("truncweibull"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
okay1 <- all(is.finite(Alpha)) && all(0 < Alpha) &&
all(is.finite(Betaa)) && all(0 < Betaa)
okay1
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
deriv = eval(substitute(expression({
M1 <- 2
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
Shape <- Betaa
Scale <- 1 / Alpha^(1/Betaa)
TTT <- extra$lower.limit
dl.dAlpha <- 1 / Alpha - y^Betaa + TTT^Betaa
dl.dBetaa <- (1 / Betaa) + log(y) -
Alpha * (y^Betaa * log(y) -
TTT^Betaa * log(TTT))
dAlpha.deta <- dtheta.deta(Alpha, .lAlpha, earg = .eAlpha )
dBetaa.deta <- dtheta.deta(Betaa, .lBetaa, earg = .eBetaa )
myderiv <- c(w) * cbind(dl.dAlpha * dAlpha.deta,
dl.dBetaa * dBetaa.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
weight = eval(substitute(expression({
aTb <- Alpha * TTT^Betaa
exp.aTb <- exp(aTb)
TblogT <- (TTT^Betaa) * log(TTT)
wingo3 <- pgamma.deriv.unscaled(q = aTb,
shape = 2) # 3-cols
Eyblogy <- (exp.aTb * (digamma(2) - wingo3[, 2]) -
(aTb + 1) * log(Alpha)) / (Alpha * Betaa)
Eyblog2y <- (exp.aTb * (digamma(2)^2 + trigamma(2) -
wingo3[, 3]) - 2 * log(Alpha) *
(digamma(2) - wingo3[, 2])) / (Alpha * Betaa^2) +
(log(Alpha)^2) * (aTb + 1) / (Alpha * Betaa^2)
ned2l.daa <- 1 / Alpha^2
ned2l.dab <- Eyblogy - TblogT
ned2l.dbb <- (1 / Betaa)^2 + Alpha * Eyblog2y -
aTb * (log(TTT))^2
wz <- array(c(c(w) * ned2l.daa * dAlpha.deta^2,
c(w) * ned2l.dbb * dBetaa.deta^2,
c(w) * ned2l.dab * dBetaa.deta * dAlpha.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .nrfs = nrfs ))))
} # truncweibull
Try the VGAM package in your browser
Any scripts or data that you put into this service are public.
VGAM documentation built on Sept. 18, 2024, 9:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions truncweibull pgamma.deriv.unscaled is.na.SurvS4 is.SurvS4 weibullR weibull.mean cens.rayleigh cens.normal ecens.exponential cens.poisson
Documented in cens.normal cens.poisson cens.poisson cens.rayleigh is.na.SurvS4 is.SurvS4 pgamma.deriv.unscaled pgamma.deriv.unscaled truncweibull weibull.mean weibullR
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
cens.poisson <-
function(link = "loglink", imu = NULL,
biglambda = 10,
smallno = 1e-10) {
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Censored Poisson distribution\n\n",
"Link: ", namesof("mu", link, earg = earg), "\n",
"Variance: mu"),
infos = eval(substitute(function(...) {
list(M1 = 1,
expected = FALSE,
multipleResponses = FALSE,
parameters.names = c("mu"),
link = .link ,
biglambda = .biglambda ,
smallno = .smallno )
}, list( .link = link,
.biglambda = biglambda,
.smallno = smallno ))),
initialize = eval(substitute(expression({
if (anyNA(y))
stop("NAs are not allowed in the response")
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 3,
Is.integer.y = TRUE)
centype <- attr(y, "type")
if (centype == "right") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- rep_len(FALSE, n)
extra$interval <- rep_len(FALSE, n)
init.mu <- pmax(y[, 1], 1/8)
} else
if (centype == "left") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- rep_len(FALSE, n)
extra$leftcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
init.mu <- pmax(y[, 1], 1/8)
} else
if (centype == "interval" ||
centype == "interval2") {
temp <- y[, 3]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- ifelse(temp == 2, TRUE, FALSE)
extra$intervalcensored <- ifelse(temp == 3, TRUE, FALSE)
init.mu <- pmax((y[, 1] + y[, 2])/2, 1/8) # intervalcensored
if (any(extra$uncensored))
init.mu[extra$uncensored] <-
pmax(y[extra$uncensored, 1], 1/8)
if (any(extra$rightcensored))
init.mu[extra$rightcensored] <-
pmax(y[extra$rightcensored, 1], 1/8)
if (any(extra$leftcensored))
init.mu[extra$leftcensored] <-
pmax(y[extra$leftcensored, 1], 1/8)
} else
if (centype == "counting") {
stop("type == 'counting' not compatible with cens.poisson()")
init.mu <- pmax(y[, 1], 1/8)
stop("currently not working")
} else
stop("response have to be in a class of SurvS4")
if (length( .imu )) init.mu <- 0 * y[, 1] + .imu
predictors.names <-
namesof("mu", .link, earg = .earg, short = TRUE)
if (!length(etastart))
etastart <- theta2eta(init.mu, link = .link , earg = .earg )
}), list( .link = link, .earg = earg, .imu = imu))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu <- eta2theta(eta, link = .link , earg = .earg )
as.vector(mu)
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c("mu" = .link )
misc$earg <- list("mu" = .earg )
}), list( .link = link, .earg = earg ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
loglikelihood = function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$intervalcensored
if (residuals){
stop("loglikelihood residuals not implemented yet")
} else {
sum(w[cen0] * dpois(y[cen0, 1], mu[cen0], log = TRUE)) +
sum(w[cenU] * log1p(-ppois(y[cenU, 1] - 1, mu[cenU]))) +
sum(w[cenL] * ppois(y[cenL, 1] - 1, mu[cenL], log.p = TRUE)) +
sum(w[cenI] * log(ppois(y[cenI, 2], mu[cenI]) -
ppois(y[cenI, 1], mu[cenI])))
}
},
vfamily = "cens.poisson",
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta, link = .link , earg = .earg )
okay1 <- all(is.finite(lambda)) && all(0 < lambda)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$intervalcensored
lambda <- eta2theta(eta, link = .link , earg = .earg )
dl.dlambda <- y[, 1] / lambda - 1 # uncensored
yllim <- yulim <- y[, 1] # uncensored
if (any(cenU)) {
yllim[cenU] <- y[cenU, 1]
densm1 <- dpois(yllim-1, lambda)
queue <- ppois(yllim-1, lambda, lower.tail = FALSE)
dl.dlambda[cenU] <- densm1[cenU] / queue[cenU]
if (any(fix.up <- (is.na(dl.dlambda) |
queue < .smallno ) &
cenU &
lambda > .biglambda )) {
zedd <- (yllim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
dl.dlambda[fix.up] <-
mills.ratio(-zedd) / sqrt(lambda[fix.up])
}
} # cenU
if (any(cenL)) {
yulim[cenL] <- y[cenL, 1] - 1
densm0 <- dpois(yulim, lambda)
Queue <- ppois(yulim, lambda) # Left tail probability
dl.dlambda[cenL] <- -densm0[cenL] / Queue[cenL]
if (any(fix.up <- (is.na(dl.dlambda) |
Queue < .smallno ) &
cenL &
lambda > .biglambda )) {
dl.dlambda[fix.up] <-
-mills.ratio(yulim[fix.up] / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])) / sqrt(lambda[fix.up])
}
} # cenL
if (any(cenI)) {
yllim[cenI] <- y[cenI, 1] + 1
yulim[cenI] <- y[cenI, 2]
Queue1 <- ppois(yllim-1, lambda)
Queue2 <- ppois(yulim , lambda)
densm02 <- dpois(yulim , lambda)
densm12 <- dpois(yllim-1, lambda)
dl.dlambda[cenI] <-
(-densm02[cenI]+densm12[cenI]) / (Queue2[cenI]-Queue1[cenI])
} # cenI
dlambda.deta <- dtheta.deta(theta = lambda,
link = .link , earg = .earg )
der1 <- c(w) * dl.dlambda * dlambda.deta
der1
}), list( .link = link, .earg = earg,
.biglambda = biglambda ,
.smallno = smallno ))),
weight = eval(substitute(expression({
d2lambda.deta2 <- d2theta.deta2(theta = lambda,
link = .link , earg = .earg )
d2l.dlambda2 <- 1 / lambda # uncensored; Fisher scoring
if (any(cenU)) {
densm2 <- dpois(yllim-2, lambda)
d2l.dlambda2[cenU] <- (dl.dlambda[cenU])^2 -
(densm2[cenU] - densm1[cenU]) / queue[cenU]
if (any(fix.up <- is.na(d2l.dlambda2) | fix.up)) {
warning("incorrect as one of the queue uses yulim, not yulim-1")
zeddm1 <- (yllim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
zeddm2 <- (yllim[fix.up] - 2) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])
d2l.dlambda2[fix.up] <- (dl.dlambda[fix.up])^2 -
(mills.ratio(-zeddm2) -
mills.ratio(-zeddm1)) / sqrt(lambda[fix.up])
}
} # cenU
if (any(cenL)) {
densm1 <- dpois(yulim-1, lambda)
d2l.dlambda2[cenL] <- (dl.dlambda[cenL])^2 -
(densm0[cenL] - densm1[cenL]) / Queue[cenL]
if (any(fix.up <- is.na(d2l.dlambda2) | fix.up)) {
d2l.dlambda2[fix.up] <- (dl.dlambda[fix.up])^2 -
(mills.ratio(yulim[fix.up] / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up])) -
mills.ratio((yulim[fix.up] - 1) / sqrt(lambda[fix.up]) -
sqrt(lambda[fix.up]))) / sqrt(lambda[fix.up])
}
} # cenL
if (any(cenI)) {
densm03 <- dpois(yulim-1, lambda)
densm13 <- dpois(yllim-2, lambda)
d2l.dlambda2[cenI] <- (dl.dlambda[cenI])^2 -
(densm13[cenI]-densm12[cenI]-densm03[cenI] +
densm02[cenI]) / (Queue2[cenI]-Queue1[cenI])
} # cenI
wz <- c(w) * (dlambda.deta^2) * d2l.dlambda2
wz
}), list( .link = link, .earg = earg,
.biglambda = biglambda ,
.smallno = smallno ))))
} # cens.poisson
if (FALSE)
cens.exponential <-
ecens.exponential <- function(link = "loglink", location = 0) {
if (!is.Numeric(location, length.arg = 1))
stop("bad input for 'location'")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Censored exponential distribution\n\n",
"Link: ", namesof("rate", link, tag = TRUE), "\n",
"Mean: ", "mu = ", location, " + 1 / ",
namesof("rate", link, tag = FALSE), "\n",
"Variance: ",
if (location == 0) "Exponential: mu^2" else
paste("(mu-", location, ")^2", sep = "")),
initialize = eval(substitute(expression({
extra$location <- .location
if (any(y[, 1] <= extra$location))
stop("all responses must be greater than ", extra$location)
predictors.names <- namesof("rate", .link , .earg , tag = FALSE)
type <- attr(y, "type")
if (type == "right" || type == "left"){
mu <- y[, 1] + (abs(y[, 1] - extra$location) < 0.001) / 8
} else if (type == "interval") {
temp <- y[, 3]
mu <- ifelse(temp == 3,
y[, 2] + (abs(y[, 2] - extra$location) < 0.001) / 8,
y[, 1] + (abs(y[, 1] - extra$location) < 0.001) / 8)
}
if (!length(etastart))
etastart <- theta2eta(1/(mu-extra$location), .link , .earg )
if (type == "right") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- rep_len(FALSE, n)
extra$interval <- rep_len(FALSE, n)
} else
if (type == "left") {
temp <- y[, 2]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- rep_len(FALSE, n)
extra$leftcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
} else
if (type == "counting") {
stop("type == 'counting' not recognized")
extra$uncensored <- rep(temp == 1, TRUE, FALSE)
extra$interval <- rep_len(FALSE, n)
extra$leftcensored <- rep_len(FALSE, n)
extra$rightcensored <- rep_len(FALSE, n)
extra$counting <- ifelse(temp == 0, TRUE, FALSE)
} else
if (type == "interval") {
temp <- y[, 3]
extra$uncensored <- ifelse(temp == 1, TRUE, FALSE)
extra$rightcensored <- ifelse(temp == 0, TRUE, FALSE)
extra$leftcensored <- ifelse(temp == 2, TRUE, FALSE)
extra$interval <- ifelse(temp == 3, TRUE, FALSE)
} else
stop("'type' not recognized")
}), list( .location = location, .link = link ))),
linkinv = eval(substitute(function(eta, extra = NULL)
extra$location + 1 / eta2theta(eta, .link , .earg ),
list( .link = link ) )),
last = eval(substitute(expression({
misc$location <- extra$location
misc$link <- c("rate" = .link)
misc$multipleResponses <- FALSE
}), list( .link = link ))),
link = eval(substitute(function(mu, extra = NULL)
theta2eta(1 / (mu - extra$location), .link , .earg ),
list( .link = link ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
rate <- 1 / (mu - extra$location)
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$interval
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cenL] * log1p(-exp(-rate[cenL] *
(y[cenL, 1] - extra$location)))) +
sum(w[cenU] * (-rate[cenU]*(y[cenU, 1]-extra$location))) +
sum(w[cen0] * (log(rate[cen0]) -
rate[cen0]*(y[cen0, 1]-extra$location))) +
sum(w[cenI] * log(-exp(-rate[cenI]*(y[cenI, 2]-extra$location))+
exp(-rate[cenI]*(y[cenI, 1]-extra$location))))
}, list( .link = link ))),
vfamily = c("ecens.exponential"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
rate <- 1 / (mu - extra$location)
okay1 <- all(is.finite(rate)) && all(0 < rate)
okay1
}, list( .link = link ))),
deriv = eval(substitute(expression({
rate <- 1 / (mu - extra$location)
cen0 <- extra$uncensored
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cenI <- extra$interval
dl.drate <- 1/rate - (y[, 1]-extra$location) # uncensored
tmp200 <- exp(-rate*(y[, 1]-extra$location))
tmp200b <- exp(-rate*(y[, 2]-extra$location)) # interval censored
if (any(cenL))
dl.drate[cenL] <- (y[cenL, 1]-extra$location) *
tmp200[cenL] / (1 - tmp200[cenL])
if (any(cenU))
dl.drate[cenU] <- -(y[cenU, 1]-extra$location)
if (any(cenI))
dl.drate[cenI] <- ((y[cenI, 2] - extra$location) *
tmp200b[cenI] - (y[cenI, 1] - extra$location) *
tmp200[cenI]) / (-tmp200b[cenI] + tmp200[cenI])
drate.deta <- dtheta.deta(rate, .link , .earg )
c(w) * dl.drate * drate.deta
}), list( .link = link ) )),
weight = eval(substitute(expression({
A123 <- ((mu-extra$location)^2) # uncensored d2l.drate2
Lowpt <- ifelse(cenL, y[, 1], extra$location)
Lowpt <- ifelse(cenI, y[, 1], Lowpt) #interval censored
Upppt <- ifelse(cenU, y[, 1], Inf)
Upppt <- ifelse(cenI, y[, 2], Upppt) #interval censored
tmp300 <- exp(-rate*(Lowpt - extra$location))
d2l.drate2 <- 0 * y[, 1]
ind50 <- Lowpt > extra$location
d2l.drate2[ind50] <- (Lowpt[ind50]-extra$location)^2 *
tmp300[ind50] / (1-tmp300[ind50])
d2l.drate2 <- d2l.drate2 + (exp(-rate*(Lowpt-extra$location)) -
exp(-rate*(Upppt-extra$location))) * A123
wz <- c(w) * (drate.deta^2) * d2l.drate2
wz
}), list( .link = link ))))
} # cens.exponential
cennormal <-
cens.normal <- function(lmu = "identitylink", lsd = "loglink",
imethod = 1, zero = "sd") {
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsd <- as.list(substitute(lsd))
esd <- link2list(lsd)
lsd <- attr(esd, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
new("vglmff",
blurb = c("Censored univariate normal\n\n",
"Links: ", namesof("mu", lmu, tag = TRUE), "; ",
namesof("sd", lsd, tag = TRUE), "\n",
"Conditional variance: sd^2"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
zero = .zero ,
multiple.responses = FALSE,
parameters.names = c("mu", "sd"),
expected = TRUE )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
if (!length(extra$leftcensored))
extra$leftcensored <- rep_len(FALSE, n)
if (!length(extra$rightcensored))
extra$rightcensored <- rep_len(FALSE, n)
if (any(extra$rightcensored & extra$leftcensored))
stop("some observations are both right and left censored!")
predictors.names <-
c(namesof("mu", .lmu , earg = .emu , tag = FALSE),
namesof("sd", .lsd , earg = .esd , tag = FALSE))
if (!length(etastart)) {
anyc <- extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # can be all data
junk <- lm.wfit(x = cbind(x[!i11, ]),
y = y[!i11], w = w[!i11])
sd.y.est <- sqrt(sum(w[!i11] * junk$resid^2) / junk$df.residual)
etastart <- cbind(mu = y,
rep_len(theta2eta(sd.y.est, .lsd), n))
if (any(anyc))
etastart[anyc, 1] <- x[anyc, , drop = FALSE] %*% junk$coeff
}
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd,
.imethod = imethod ))),
linkinv = eval(substitute( function(eta, extra = NULL) {
eta2theta(eta[, 1], .lmu , earg = .emu )
}, list( .lmu = lmu, .emu = emu ))),
last = eval(substitute(expression({
misc$link <- c("mu" = .lmu , "sd" = .lsd )
misc$earg <- list("mu" = .emu , "sd" = .esd )
misc$expected <- TRUE
misc$multipleResponses <- FALSE
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cen0 <- !cenL & !cenU # uncensored obsns
mum <- eta2theta(eta[, 1], .lmu , earg = .emu )
sdv <- eta2theta(eta[, 2], .lsd , earg = .esd )
Lower <- ifelse(cenL, y, -Inf)
Upper <- ifelse(cenU, y, Inf)
ell1 <- -log(sdv[cen0]) -
0.5 * ((y[cen0] - mum[cen0])/sdv[cen0])^2
ell2 <- log1p(-pnorm((mum[cenL] - Lower[cenL]) / sdv[cenL]))
ell3 <- log1p(-pnorm(( Upper[cenU] - mum[cenU]) / sdv[cenU]))
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cen0] * ell1) +
sum(w[cenL] * ell2) +
sum(w[cenU] * ell3)
}, list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
vfamily = c("cens.normal"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mum <- eta2theta(eta[, 1], .lmu )
sdv <- eta2theta(eta[, 2], .lsd )
okay1 <- all(is.finite(mum)) &&
all(is.finite(sdv)) && all(0 < sdv)
okay1
}, list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
deriv = eval(substitute(expression({
cenL <- extra$leftcensored
cenU <- extra$rightcensored
cen0 <- !cenL & !cenU # uncensored obsns
Lower <- ifelse(cenL, y, -Inf)
Upper <- ifelse(cenU, y, Inf)
mum <- eta2theta(eta[, 1], .lmu )
sdv <- eta2theta(eta[, 2], .lsd )
dl.dmu <- (y-mum) / sdv^2
dl.dsd <- (((y-mum)/sdv)^2 - 1) / sdv
dmu.deta <- dtheta.deta(mum, .lmu , earg = .emu )
dsd.deta <- dtheta.deta(sdv, .lsd , earg = .esd )
if (any(cenL)) {
mumL <- mum - Lower
temp21L <- mumL[cenL] / sdv[cenL]
PhiL <- pnorm(temp21L)
phiL <- dnorm(temp21L)
fred21 <- phiL / (1 - PhiL)
dl.dmu[cenL] <- -fred21 / sdv[cenL]
dl.dsd[cenL] <- mumL[cenL] * fred21 / sdv[cenL]^2
rm(fred21)
}
if (any(cenU)) {
mumU <- Upper - mum
temp21U <- mumU[cenU] / sdv[cenU]
PhiU <- pnorm(temp21U)
phiU <- dnorm(temp21U)
fred21 <- phiU / (1 - PhiU)
dl.dmu[cenU] <- fred21 / sdv[cenU] # Negated
dl.dsd[cenU] <- mumU[cenU] * fred21 / sdv[cenU]^2
rm(fred21)
}
c(w) * cbind(dl.dmu * dmu.deta,
dl.dsd * dsd.deta)
}), list( .lmu = lmu, .lsd = lsd,
.emu = emu, .esd = esd ))),
weight = eval(substitute(expression({
A1 <- 1 - pnorm((mum - Lower) / sdv) # Lower
A3 <- 1 - pnorm((Upper - mum) / sdv) # Upper
A2 <- 1 - A1 - A3 # Middle; uncensored
wz <- matrix(0, n, 3)
wz[, iam(1, 1,M)] <- A2 * 1 / sdv^2 # ed2l.dmu2
wz[, iam(2, 2,M)] <- A2 * 2 / sdv^2 # ed2l.dsd2
mumL <- mum - Lower
temp21L <- mumL / sdv
PhiL <- pnorm(temp21L)
phiL <- dnorm(temp21L)
temp31L <- ((1-PhiL) * sdv)^2
wz.cenL11 <- phiL * (phiL - (1-PhiL)*temp21L) / temp31L
wz.cenL22 <- mumL * phiL * ((1-PhiL) * (2 - temp21L^2) +
mumL * phiL / sdv) / (sdv * temp31L)
wz.cenL12 <- phiL * ((1-PhiL)*(temp21L^2 - 1) -
temp21L*phiL) / temp31L
wz.cenL11[!is.finite(wz.cenL11)] <- 0
wz.cenL22[!is.finite(wz.cenL22)] <- 0
wz.cenL12[!is.finite(wz.cenL12)] <- 0
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + A1 * wz.cenL11
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] + A1 * wz.cenL22
wz[, iam(1, 2, M)] <- A1 * wz.cenL12
mumU <- Upper - mum # often Inf
temp21U <- mumU / sdv # often Inf
PhiU <- pnorm(temp21U) # often 1
phiU <- dnorm(temp21U) # often 0
temp31U <- ((1-PhiU) * sdv)^2 # often 0
tmp8 <- (1-PhiU)*temp21U
wzcenU11 <- phiU * (phiU - tmp8) / temp31U
tmp9 <- (1-PhiU) * (2 - temp21U^2)
wzcenU22 <- mumU * phiU * (tmp9 +
mumU * phiU / sdv) / (sdv * temp31U)
wzcenU12 <- -phiU * ((1-PhiU)*(temp21U^2 - 1) -
temp21U*phiU) / temp31U
wzcenU11[!is.finite(wzcenU11)] <- 0 # Needed when Upper==Inf
wzcenU22[!is.finite(wzcenU22)] <- 0 # Needed when Upper==Inf
wzcenU12[!is.finite(wzcenU12)] <- 0 # Needed when Upper==Inf
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] + A3 * wzcenU11
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] + A3 * wzcenU22
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] + A3 * wzcenU12
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dmu.deta^2
wz[, iam(2, 2, M)] <- wz[, iam(2, 2, M)] * dsd.deta^2
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] * dmu.deta * dsd.deta
c(w) * wz
}), list( .lmu = lmu, .lsd = lsd ))))
} # cens.normal
cens.rayleigh <- function(lscale = "loglink",
oim = TRUE) {
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.logical(oim) || length(oim) != 1)
stop("bad input for argument 'oim'")
new("vglmff",
blurb = c("Censored Rayleigh distribution\n\n",
"f(y) = y*exp(-0.5*(y/scale)^2)/scale^2, y>0, scale>0\n",
"Link: ",
namesof("scale", lscale, earg = escale ), "\n", "\n",
"Mean: scale * sqrt(pi / 2)"),
initialize = eval(substitute(expression({
if (NCOL(y) != 1)
stop("response must be a vector or a one-column matrix")
if (length(extra$leftcensored))
stop("cannot handle left-censored data")
if (!length(extra$rightcensored))
extra$rightcensored <- rep_len(FALSE, n)
predictors.names <-
namesof("scale", .lscale , earg = .escale , tag = FALSE)
if (!length(etastart)) {
a.init <- (y+1/8) / sqrt(pi/2)
etastart <- theta2eta(a.init, .lscale , earg = .escale )
}
}), list( .lscale = lscale, .escale = escale ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
Scale * sqrt(pi/2)
}, list( .lscale = lscale, .escale = escale ))),
last = eval(substitute(expression({
misc$link <- c("scale" = .lscale )
misc$earg <- list("scale" = .escale )
misc$oim <- .oim
}), list( .lscale = lscale, .escale = escale,
.oim = oim ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
cen0 <- !extra$rightcensored # uncensored obsns
cenU <- extra$rightcensored
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else
sum(w[cen0] * (log(y[cen0]) - 2*log(Scale[cen0]) -
0.5*(y[cen0]/Scale[cen0])^2)) -
sum(w[cenU] * (y[cenU]/Scale[cenU])^2) * 0.5
}, list( .lscale = lscale, .escale = escale ))),
vfamily = c("cens.rayleigh"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Scale <- eta2theta(eta, .lscale , earg = .escale )
okay1 <- all(is.finite(Scale)) && all(0 < Scale)
okay1
}, list( .lscale = lscale, .escale = escale ))),
deriv = eval(substitute(expression({
cen0 <- !extra$rightcensored # uncensored obsns
cenU <- extra$rightcensored
Scale <- eta2theta(eta, .lscale , earg = .escale )
dl.dScale <- ((y/Scale)^2 - 2) / Scale
dScale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
dl.dScale[cenU] <- y[cenU]^2 / Scale[cenU]^3
c(w) * dl.dScale * dScale.deta
}), list( .lscale = lscale,
.escale = escale ))),
weight = eval(substitute(expression({
ned2l.dScale2 <- 4 / Scale^2
wz <- dScale.deta^2 * ned2l.dScale2
if ( .oim ) {
d2l.dScale2 <- 3 * (y[cenU])^2 / (Scale[cenU])^4
d2Scale.deta2 <- d2theta.deta2(Scale[cenU], .lscale ,
earg = .escale )
wz[cenU] <- (dScale.deta[cenU])^2 * d2l.dScale2 -
dl.dScale[cenU] * d2Scale.deta2
} else {
ned2l.dScale2[cenU] <- 6 / (Scale[cenU])^2
wz[cenU] <- (dScale.deta[cenU])^2 * ned2l.dScale2[cenU]
}
c(w) * wz
}), list( .lscale = lscale, .escale = escale,
.oim = oim ))))
} # cens.rayleigh
weibull.mean <-
function(lmean = "loglink", lshape = "loglink",
imean = NULL, ishape = NULL,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1, zero = "shape") {
imeann <- imean
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
lmeann <- as.list(substitute(lmean))
emeann <- link2list(lmeann)
lmeann <- attr(emeann, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (length(ishape))
if (!is.Numeric(ishape, positive = TRUE))
stop("argument 'ishape' values must be positive")
if (length(imeann))
if (!is.Numeric(imeann, positive = TRUE))
stop("argument 'imean' values must be positive")
blurb.vec <- c(namesof("mean", lmeann, earg = emeann),
namesof("shape", lshape, earg = eshape))
new("vglmff",
blurb = c("Weibull distribution (parameterized by the mean)\n\n",
"Links: ",
blurb.vec[1], ", ",
blurb.vec[2], "\n",
"Mean: mean\n",
"Variance: mean^2 * (gamma(1 + 2/shape) / ",
"gamma(1 + 1/shape)^2 - 1)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero,
.lmeann = lmeann ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("mean", "shape"),
lmean = .lmeann ,
lshape = .lshape ,
zero = .zero )
}, list( .zero = zero,
.lmeann = lmeann, .lshape = lshape ))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
mynames1 <- param.names("mean" , ncoly, skip1 = TRUE)
mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lmeann , earg = .emeann , tag = FALSE),
namesof(mynames2, .lshape , earg = .eshape , tag = FALSE))
predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)]
Meann.init <- matrix(if (length( .imeann ))
.imeann else 0.5 * colMeans(y),
n, ncoly, byrow = TRUE) + 0.5 * y
Shape.init <- matrix(if (length( .ishape )) .ishape else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .ishape ) ||
!length( .imeann )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
if (!is.Numeric(Shape.init[, ilocal]))
Shape.init[, ilocal] <- 1 / fit0$coef["X"]
} # ilocal
etastart <-
cbind(theta2eta(Meann.init, .lmeann , earg = .emeann ),
theta2eta(Shape.init, .lshape , earg = .eshape ))[,
interleave.VGAM(M, M1 = M1)]
}
}
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape,
.imeann = imeann, .ishape = ishape,
.probs.y = probs.y,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Meann
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
last = eval(substitute(expression({
regnotok <- any(Shape <= 2)
if (any(Shape <= 1)) {
warning("MLE regularity conditions are violated",
"(shape <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < Shape & Shape < 2)) {
warning("MLE regularity conditions are violated",
"(1 < shape < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == Shape)) {
warning("MLE regularity conditions are violated",
"(shape == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than ",
"when shape > 2")
}
M1 <- extra$M1
avector <- c(rep_len( .lmeann , ncoly),
rep_len( .lshape , ncoly))
misc$link <- avector[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .emeann
misc$earg[[M1*ii ]] <- .eshape
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$RegCondOK <- !regnotok # Save this for later
misc$expected <- TRUE # all(cen0)
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape,
.imethod = imethod ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * dweibull(x = y, shape = Shape,
scale = Meann / gamma(1 + 1/Shape),
log = TRUE))
}
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
vfamily = c("weibull.mean"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
okay1 <- all(is.finite(Meann)) && all(0 < Meann) &&
all(is.finite(Shape)) && all(0 < Shape)
okay1
}, list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
deriv = eval(substitute(expression({
M1 <- 2
Meann <- eta2theta(eta[, c(TRUE, FALSE)], .lmeann , earg = .emeann )
Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , earg = .eshape )
if (FALSE) {
} else {
EulerM <- -digamma(1.0)
AA <- (EulerM - 1)^2 + (pi^2) / 6
BB <- digamma(1 + 1/Shape)
CC <- y * gamma(1 + 1/Shape) / Meann
dl.dmeann <- (CC^Shape - 1) * Shape / Meann # Agrees
dl.dshape <- 1/Shape -
(log(y/Meann) + lgamma(1 + 1/Shape)) * (CC^Shape - 1) +
(BB / Shape) * (CC^Shape - 1)
}
dmeann.deta <- dtheta.deta(Meann, .lmeann , earg = .emeann )
dshape.deta <- dtheta.deta(Shape, .lshape , earg = .eshape )
myderiv <- c(w) * cbind(dl.dmeann * dmeann.deta,
dl.dshape * dshape.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lmeann = lmeann, .lshape = lshape,
.emeann = emeann, .eshape = eshape ) )),
weight = eval(substitute(expression({
if (FALSE) {
} else {
ned2l.dmeann <- (Shape / Meann)^2 #
ned2l.dshape <- AA / Shape^2 # Unchanged
ned2l.dshapemeann <- (EulerM - 1 + BB) / Meann
}
wz <- array(c(c(w) * ned2l.dmeann * dmeann.deta^2,
c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dshapemeann * dmeann.deta * dshape.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .eshape = eshape ))))
} # weibull.mean
weibullR <-
function(lscale = "loglink", lshape = "loglink",
iscale = NULL, ishape = NULL,
lss = TRUE,
nrfs = 1,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1, zero = "shape") {
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (!is.Numeric(nrfs, length.arg = 1) ||
nrfs < 0 ||
nrfs > 1)
stop("bad input for argument 'nrfs'")
if (length(ishape))
if (!is.Numeric(ishape, positive = TRUE))
stop("argument 'ishape' values must be positive")
if (length(iscale))
if (!is.Numeric(iscale, positive = TRUE))
stop("argument 'iscale' values must be positive")
scale.TF <- if (lss) c(TRUE, FALSE) else c(FALSE, TRUE)
scale.12 <- if (lss) 1:2 else 2:1
blurb.vec <- c(namesof("scale", lscale, earg = escale),
namesof("shape", lshape, earg = eshape))
blurb.vec <- blurb.vec[scale.12]
new("vglmff",
blurb = c("Weibull distribution\n\n",
"Links: ",
blurb.vec[1], ", ",
blurb.vec[2], "\n",
"Mean: scale * gamma(1 + 1/shape)\n",
"Variance: scale^2 * (gamma(1 + 2/shape) - ",
"gamma(1 + 1/shape)^2)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = if ( .lss )
c("scale", "shape") else
c("shape", "scale"),
lss = .lss ,
lscale = .lscale ,
lshape = .lshape ,
zero = .zero )
}, list( .zero = zero, .scale.12 = scale.12, .scale.TF = scale.TF,
.lscale = lscale ,
.lshape = lshape ,
.lss = lss
))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
if ( .lss ) {
mynames1 <- param.names("scale", ncoly, skip1 = TRUE)
mynames2 <- param.names("shape", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lscale , earg = .escale , tag = FALSE),
namesof(mynames2, .lshape , earg = .eshape , tag = FALSE))
} else {
mynames1 <- param.names("shape", ncoly, skip1 = TRUE)
mynames2 <- param.names("scale", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lshape , earg = .eshape , tag = FALSE),
namesof(mynames2, .lscale , earg = .escale , tag = FALSE))
}
predictors.names <- predictors.names[
interleave.VGAM(M, M1 = M1)]
Shape.init <- matrix(if (length( .ishape )) .ishape else 0 + NA,
n, ncoly, byrow = TRUE)
Scale.init <- matrix(if (length( .iscale )) .iscale else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .ishape ) ||
!length( .iscale )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
if (!is.Numeric(Shape.init[, ilocal]))
Shape.init[, ilocal] <- 1 / fit0$coef["X"]
if (!is.Numeric(Scale.init[, ilocal]))
Scale.init[, ilocal] <- exp(fit0$coef["Intercept"])
} # ilocal
etastart <- if ( .lss )
cbind(theta2eta(Scale.init, .lscale , earg = .escale ),
theta2eta(Shape.init, .lshape , earg = .eshape ))[,
interleave.VGAM(M, M1 = M1)] else
cbind(theta2eta(Shape.init, .lshape , earg = .eshape ),
theta2eta(Scale.init, .lscale , earg = .escale ))[,
interleave.VGAM(M, M1 = M1)]
}
}
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.iscale = iscale, .ishape = ishape,
.probs.y = probs.y,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
Scale * gamma(1 + 1 / Shape)
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
last = eval(substitute(expression({
regnotok <- any(Shape <= 2)
if (any(Shape <= 1)) {
warning("MLE regularity conditions are violated",
"(shape <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < Shape & Shape < 2)) {
warning("MLE regularity conditions are violated",
"(1 < shape < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == Shape)) {
warning("MLE regularity conditions are violated",
"(shape == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than when ",
"shape > 2")
}
M1 <- extra$M1
avector <- if ( .lss ) c(rep_len( .lscale , ncoly),
rep_len( .lshape , ncoly)) else
c(rep_len( .lshape , ncoly),
rep_len( .lscale , ncoly))
misc$link <- avector[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- if ( .lss ) .escale else .eshape
misc$earg[[M1*ii ]] <- if ( .lss ) .eshape else .escale
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$nrfs <- .nrfs
misc$RegCondOK <- !regnotok # Save this for later
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.imethod = imethod,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss,
.nrfs = nrfs ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * dweibull(y, shape = Shape, scale = Scale, log = TRUE))
}
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
vfamily = c("weibullR"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
okay1 <- all(is.finite(Scale)) && all(0 < Scale) &&
all(is.finite(Shape)) && all(0 < Shape)
okay1
}, list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
deriv = eval(substitute(expression({
M1 <- 2
Scale <- eta2theta(eta[, .scale.TF ], .lscale , earg = .escale )
Shape <- eta2theta(eta[, !( .scale.TF )], .lshape , earg = .eshape )
dl.dshape <- 1 / Shape + log(y / Scale) -
log(y / Scale) * (y / Scale)^Shape
dl.dscale <- (Shape / Scale) * (-1.0 + (y / Scale)^Shape)
dshape.deta <- dtheta.deta(Shape, .lshape, earg = .eshape )
dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
myderiv <- if ( .lss )
c(w) * cbind(dl.dscale * dscale.deta,
dl.dshape * dshape.deta) else
c(w) * cbind(dl.dshape * dshape.deta,
dl.dscale * dscale.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lscale = lscale, .lshape = lshape,
.escale = escale, .eshape = eshape,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ) )),
weight = eval(substitute(expression({
EulerM <- -digamma(1.0)
ned2l.dscale <- (Shape / Scale)^2
ned2l.dshape <- (6*(EulerM - 1)^2 + pi^2)/(6*Shape^2) # KK (2003)
ned2l.dshapescale <- (EulerM-1) / Scale
wz <- if ( .lss )
array(c(c(w) * ned2l.dscale * dscale.deta^2,
c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dshapescale * dscale.deta * dshape.deta),
dim = c(n, M / M1, 3)) else
array(c(c(w) * ned2l.dshape * dshape.deta^2,
c(w) * ned2l.dscale * dscale.deta^2,
c(w) * ned2l.dshapescale * dscale.deta * dshape.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .eshape = eshape, .nrfs = nrfs,
.scale.12 = scale.12, .scale.TF = scale.TF, .lss = lss ))))
} # weibullR
setOldClass(c("SurvS4", "Surv"))
SurvS4 <-
function (time, time2, event, type = c("right", "left", "interval",
"counting", "interval2"), origin = 0) {
nn <- length(time)
ng <- nargs()
if (missing(type)) {
if (ng == 1 || ng == 2)
type <- "right" else if (ng == 3)
type <- "counting" else stop("Invalid number of arguments")
} else {
type <- match.arg(type)
ng <- ng - 1
if (ng != 3 && (type == "interval" || type == "counting"))
stop("Wrong number of args for this type of survival data")
if (ng != 2 && (type == "right" || type == "left" ||
type == "interval2"))
stop("Wrong number of args for this type of survival data")
}
who <- !is.na(time)
if (ng == 1) {
if (!is.numeric(time))
stop("Time variable is not numeric")
ss <- cbind(time, 1)
dimnames(ss) <- list(NULL, c("time", "status"))
} else if (type == "right" || type == "left") {
if (!is.numeric(time))
stop("Time variable is not numeric")
if (length(time2) != nn)
stop("Time and status are different lengths")
if (is.logical(time2))
status <- 1 * time2 else if (is.numeric(time2)) {
who2 <- !is.na(time2)
if (max(time2[who2]) == 2)
status <- time2 - 1 else status <- time2
if (any(status[who2] != 0 & status[who2] != 1))
stop("Invalid status value")
} else stop("Invalid status value")
ss <- cbind(time, status)
dimnames(ss) <- list(NULL, c("time", "status"))
} else if (type == "counting") {
if (length(time2) != nn)
stop("Start and stop are different lengths")
if (length(event) != nn)
stop("Start and event are different lengths")
if (!is.numeric(time))
stop("Start time is not numeric")
if (!is.numeric(time2))
stop("Stop time is not numeric")
who3 <- who & !is.na(time2)
if (any(time[who3] >= time2[who3]))
stop("Stop time must be > start time")
if (is.logical(event))
status <- 1 * event else if (is.numeric(event)) {
who2 <- !is.na(event)
if (max(event[who2]) == 2)
status <- event - 1 else status <- event
if (any(status[who2] != 0 & status[who2] != 1))
stop("Invalid status value")
} else stop("Invalid status value")
ss <- cbind(time - origin, time2 - origin, status)
} else {
if (type == "interval2") {
event <- ifelse(is.na(time), 2, ifelse(is.na(time2),
0, ifelse(time == time2, 1, 3)))
if (any(time[event == 3] > time2[event == 3]))
stop("Invalid interval: start > stop")
time <- ifelse(event != 2, time, time2)
type <- "interval"
} else {
temp <- event[!is.na(event)]
if (!is.numeric(temp))
stop("Status indicator must be numeric")
if (length(temp) > 0 && any(temp != floor(temp) |
temp < 0 | temp > 3))
stop("Status indicator must be 0, 1, 2 or 3")
}
status <- event
ss <- cbind(time, ifelse(!is.na(event) & event == 3,
time2, 1), status)
}
attr(ss, "type") <- type
class(ss) <- "SurvS4"
ss
} # SurvS4
is.SurvS4 <- function(x) inherits(x, "SurvS4")
setIs(class1 = "SurvS4",
class2 = "matrix") # Forces vglm()@y to be a matrix
as.character.SurvS4 <- function (x, ...) {
class(x) <- NULL
type <- attr(x, "type")
if (type == "right") {
temp <- x[, 2]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "+", " "))
paste(format(x[, 1]), temp, sep = "")
} else if (type == "counting") {
temp <- x[, 3]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "+", " "))
paste("(", format(x[, 1]), ",", format(x[, 2]), temp, "]", sep = "")
} else if (type == "left") {
temp <- x[, 2]
temp <- ifelse(is.na(temp), "?", ifelse(temp == 0, "<", " "))
paste(temp, format(x[, 1]), sep = "")
} else {
stat <- x[, 3]
temp <- c("+", "", "-", "]")[stat + 1]
temp2 <- ifelse(stat == 3, paste("(", format(x[, 1]),
", ", format(x[, 2]), sep = ""), format(x[, 1]))
ifelse(is.na(stat), as.character(NA), paste(temp2, temp, sep = ""))
}
} # as.character.SurvS4
"[.SurvS4" <- function(x, i, j, drop = FALSE) {
if (missing(j)) {
temp <- class(x)
type <- attr(x, "type")
class(x) <- NULL
x <- x[i, , drop = FALSE]
class(x) <- temp
attr(x, "type") <- type
x
} else {
class(x) <- NULL
NextMethod("[")
}
} # "[.SurvS4"
is.na.SurvS4 <- function(x) {
as.vector( (1* is.na(unclass(x)))%*% rep(1, ncol(x)) >0)
}
show.SurvS4 <- function (object)
print.default(as.character.SurvS4(object), quote = FALSE)
setMethod("show", "SurvS4",
function(object)
show.SurvS4(object))
pgamma.deriv.unscaled <- function(q, shape) {
gam0 <- exp(lgamma(shape) + pgamma(q = q, shape = shape,
log.p = TRUE))
I.sq <- pgamma(q = q, shape = shape)
alld <- pgamma.deriv(q = q, shape = shape) # 6-coln matrix
tmp3 <- alld[, 3] / I.sq # RHS of eqn (4.5) of \cite{wing:1989}
G1s <- digamma(shape) + tmp3 # eqn (4.9)
gam1 <- gam0 * G1s
dG1s <- trigamma(shape) + alld[, 4] / I.sq - tmp3^2 # eqn (4.13)
G2s <- dG1s + G1s^2 # eqn (4.12)
gam2 <- gam0 * G2s
cbind("0" = gam0,
"1" = gam1,
"2" = gam2)
} # pgamma.deriv.unscaled
truncweibull <-
function(lower.limit = 1e-5,
lAlpha = "loglink", lBetaa = "loglink",
iAlpha = NULL, iBetaa = NULL,
nrfs = 1,
probs.y = c(0.2, 0.5, 0.8),
imethod = 1,
zero = "Betaa") {
lAlpha <- as.list(substitute(lAlpha))
eAlpha <- link2list(lAlpha)
lAlpha <- attr(eAlpha, "function.name")
lBetaa <- as.list(substitute(lBetaa))
eBetaa <- link2list(lBetaa)
lBetaa <- attr(eBetaa, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(probs.y, positive = TRUE) ||
length(probs.y) < 2 ||
max(probs.y) >= 1)
stop("bad input for argument 'probs.y'")
if (!is.Numeric(nrfs, length.arg = 1) ||
nrfs < 0 ||
nrfs > 1)
stop("bad input for argument 'nrfs'")
if (length(iAlpha))
if (!is.Numeric(iAlpha, positive = TRUE))
stop("argument 'iAlpha' values must be positive")
if (length(iBetaa))
if (!is.Numeric(iBetaa, positive = TRUE))
stop("argument 'iBetaa' values must be positive")
new("vglmff",
blurb = c("Truncated weibull distribution\n\n",
"Links: ",
namesof("Alpha", lAlpha, earg = eAlpha), ", ",
namesof("Betaa", lBetaa, earg = eBetaa), "\n",
if (length( lower.limit ) < 5)
paste("Truncation point(s): ",
lower.limit, sep = ", ") else
""),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("Alpha", "Betaa"),
lower.limit = .lower.limit ,
lAlpha = .lAlpha ,
lBetaa = .lBetaa ,
zero = .zero )
}, list( .zero = zero,
.lAlpha = lAlpha ,
.lBetaa = lBetaa ,
.lower.limit = lower.limit
))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 2
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
extra$lower.limit <- matrix( .lower.limit , n, ncoly, byrow = TRUE)
if (any(y < extra$lower.limit)) {
stop("some response values less than argument 'lower.limit'")
}
if (is.SurvS4(y))
stop("only uncensored observations are allowed; ",
"don't use SurvS4()")
mynames1 <- param.names("Alpha", ncoly, skip1 = TRUE)
mynames2 <- param.names("Betaa", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lAlpha , earg = .eAlpha , tag = FALSE),
namesof(mynames2, .lBetaa , earg = .eBetaa , tag = FALSE))[
interleave.VGAM(M, M1 = M1)]
Alpha.init <- matrix(if (length( .iAlpha )) .iAlpha else 0 + NA,
n, ncoly, byrow = TRUE)
Betaa.init <- matrix(if (length( .iBetaa )) .iBetaa else 0 + NA,
n, ncoly, byrow = TRUE)
if (!length(etastart)) {
if (!length( .iAlpha ) ||
!length( .iBetaa )) {
for (ilocal in 1:ncoly) {
anyc <- FALSE # extra$leftcensored | extra$rightcensored
i11 <- if ( .imethod == 1) anyc else FALSE # Can be all data
probs.y <- .probs.y
xvec <- log(-log1p(-probs.y))
fit0 <- lsfit(x = xvec,
y = log(quantile(y[!i11, ilocal],
probs = probs.y )))
aaa.init <- 1 / fit0$coef["X"]
bbb.init <- exp(fit0$coef["Intercept"])
if (!is.Numeric(Betaa.init[, ilocal]))
Betaa.init[, ilocal] <- aaa.init
if (!is.Numeric(Alpha.init[, ilocal]))
Alpha.init[, ilocal] <- (1 / bbb.init)^aaa.init
} # ilocal
} else {
Alpha.init <- rep_len( .iAlpha , n)
Betaa.init <- rep_len( .iBetaa , n)
}
etastart <-
cbind(theta2eta(Alpha.init, .lAlpha , earg = .eAlpha ),
theta2eta(Betaa.init, .lBetaa , earg = .eBetaa ))[,
interleave.VGAM(M, M1 = M1)]
}
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.iBetaa = iBetaa, .iAlpha = iAlpha,
.lower.limit = lower.limit,
.probs.y = probs.y,
.imethod = imethod ) )),
linkinv = eval(substitute(function(eta, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
aTb <- Alpha * extra$lower.limit^Betaa
wingo3 <- pgamma.deriv.unscaled(q = aTb, shape = 1 + 1 / Betaa)
exp.aTb <- exp(aTb)
(gamma(1 + 1 / Betaa) - wingo3[, 1]) *
exp.aTb / Alpha^(1 / Betaa)
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit) )),
last = eval(substitute(expression({
aaa.hat <- Betaa
regnotok <- any(aaa.hat <= 2)
if (any(aaa.hat <= 1)) {
warning("MLE regularity conditions are violated",
"(Betaa <= 1) at the final iteration: ",
"MLEs are not consistent")
} else if (any(1 < aaa.hat & aaa.hat < 2)) {
warning("MLE regularity conditions are violated",
"(1 < Betaa < 2) at the final iteration: ",
"MLEs exist but are not asymptotically normal")
} else if (any(2 == aaa.hat)) {
warning("MLE regularity conditions are violated",
"(Betaa == 2) at the final iteration: ",
"MLEs exist and are normal and asymptotically ",
"efficient but with a slower convergence rate than when ",
"Betaa > 2")
}
M1 <- extra$M1
misc$link <-
c(rep_len( .lAlpha , ncoly),
rep_len( .lBetaa , ncoly))[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .eAlpha
misc$earg[[M1*ii ]] <- .eBetaa
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
misc$nrfs <- .nrfs
misc$RegCondOK <- !regnotok # Save this for later
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.imethod = imethod,
.lower.limit = lower.limit,
.nrfs = nrfs ) )),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
Shape <- Betaa
Scale <- 1 / Alpha^(1/Betaa)
if (residuals) stop("loglikelihood residuals not ",
"implemented yet") else {
sum(c(w) * (dweibull(x = y, shape = Shape,
scale = Scale, log = TRUE) -
pweibull(q = extra$lower.limit, shape = Shape,
scale = Scale, log.p = TRUE,
lower.tail = FALSE)))
}
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
vfamily = c("truncweibull"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
okay1 <- all(is.finite(Alpha)) && all(0 < Alpha) &&
all(is.finite(Betaa)) && all(0 < Betaa)
okay1
}, list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
deriv = eval(substitute(expression({
M1 <- 2
Alpha <- eta2theta(eta[, c(TRUE, FALSE)], .lAlpha , earg = .eAlpha )
Betaa <- eta2theta(eta[, c(FALSE, TRUE)], .lBetaa , earg = .eBetaa )
Shape <- Betaa
Scale <- 1 / Alpha^(1/Betaa)
TTT <- extra$lower.limit
dl.dAlpha <- 1 / Alpha - y^Betaa + TTT^Betaa
dl.dBetaa <- (1 / Betaa) + log(y) -
Alpha * (y^Betaa * log(y) -
TTT^Betaa * log(TTT))
dAlpha.deta <- dtheta.deta(Alpha, .lAlpha, earg = .eAlpha )
dBetaa.deta <- dtheta.deta(Betaa, .lBetaa, earg = .eBetaa )
myderiv <- c(w) * cbind(dl.dAlpha * dAlpha.deta,
dl.dBetaa * dBetaa.deta)
myderiv[, interleave.VGAM(M, M1 = M1)]
}), list( .lBetaa = lBetaa, .lAlpha = lAlpha,
.eBetaa = eBetaa, .eAlpha = eAlpha,
.lower.limit = lower.limit ) )),
weight = eval(substitute(expression({
aTb <- Alpha * TTT^Betaa
exp.aTb <- exp(aTb)
TblogT <- (TTT^Betaa) * log(TTT)
wingo3 <- pgamma.deriv.unscaled(q = aTb,
shape = 2) # 3-cols
Eyblogy <- (exp.aTb * (digamma(2) - wingo3[, 2]) -
(aTb + 1) * log(Alpha)) / (Alpha * Betaa)
Eyblog2y <- (exp.aTb * (digamma(2)^2 + trigamma(2) -
wingo3[, 3]) - 2 * log(Alpha) *
(digamma(2) - wingo3[, 2])) / (Alpha * Betaa^2) +
(log(Alpha)^2) * (aTb + 1) / (Alpha * Betaa^2)
ned2l.daa <- 1 / Alpha^2
ned2l.dab <- Eyblogy - TblogT
ned2l.dbb <- (1 / Betaa)^2 + Alpha * Eyblog2y -
aTb * (log(TTT))^2
wz <- array(c(c(w) * ned2l.daa * dAlpha.deta^2,
c(w) * ned2l.dbb * dBetaa.deta^2,
c(w) * ned2l.dab * dBetaa.deta * dAlpha.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .nrfs = nrfs ))))
} # truncweibull
Try the VGAM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.