Dear R-users,
I apply to your kind attention to know if someone have used the Splus software
koq.q (Kent & O'Quigley's measure of dependence for censored data)
in R and
kindly can help me.
I have tried several times to contact the authors Andrej Blejec
(andrej.blejec at uni-lj.si) or Janez Stare (janez.stare at mf.uni-lj.si) but
unfortunately no one answered me.
Following you'll see the function nlminb that i have changed with R-function
optim() and the error that came out when i try to run this software.
As attached file there is the file koq.q
Thanking in advance for what you can do about it
Yours Faithfully
Dr. Antonello Romani
Dipartimento di Medicina Sperimentale
Sezione di Patologia Molecolare e Immunologia
via Volturno 39
43100 Parma
Italy
#-----------------------------------------------------------------------------------------------
find.mu.alfa <- function(theta, theta1, x, ell, which,...)
{
#
# find theta=c(beta,mu,alfa) which maximize
# Expected Log-Likelihood given by ell
#
# Set lower and upper bounds for mu (-Inf,Inf)
# and alfa (0,Inf)
#
lower <- c(which, T, 0) * theta
lower[c(which, T, 0) == T] <- - Inf
upper <- c(which, T, T) * theta
upper[c(which, T, T) == T] <- Inf #
optim(par=theta,fn=ell,method="L-BFGS-B",lower = lower, upper =
upper,
x = x,theta1 =theta1)
#nlminb(start = theta, objective = ell, x = x,lower = lower, upper = upper,
theta1 = theta1, ...)
}
#
#----------------------------------------------------------------------------------------------------------
fit<-cph(Surv(futime,fustat)~age,data=ovarian,x=T,y=T,surv=T,method="breslow",type="kaplan-meier")> fit
Cox Proportional Hazards Model
cph(formula = Surv(futime, fustat) ~ age, data = ovarian, method =
"breslow",
x = T, y = T, surv = T, type = "kaplan-meier")
Obs Events Model L.R. d.f. P Score Score P
26 12 14.29 1 2e-04 12.26 5e-04
R2
0.454
coef se(coef) z p
age 0.162 0.0497 3.25 0.00116
>koq(fit)
Beta = 0
Beta1 = 0.1616199
Error in ELL(theta = theta0, theta1 = theta1, x = x) :
only 0's may mix with negative subscripts
In addition: Warning message:
Replacement length not a multiple of the elements to replace in matrix(...)
-------------- next part --------------
"koq"<-
function(beta1, x = NULL, p = NULL, verbose = T)
{
#
#--------------------------------------------------------------------
check.parameters <- function(beta1, x, p)
{
# Check parameters and
# determine variables to be tested
#
if(data.class(beta1) == "cph") x <-
coxph.detail(beta1)$x
if(is.null(x))
stop("Independent variables (x) not set")
if(data.class(beta1) == "coxreg" || data.class(
beta1) == "cph")
beta1 <- beta1$coef
if(!(data.class(beta1) == "numeric"))
stop("Illegal parameter beta1 \n A vector of numeric coefficients or an
object of class \"coxreg\" expected"
)
m <- length(beta1)
if(is.null(p)) p <- 1:m
# p is a list of vars to be tested
if(max(p) > m | (min(p) < 1 & max(p) > 1))
stop(paste("Illegal variable selected p=",
list(p)))
if(max(p) == 1 & (min(p) == 0 | max(p) == 1) &
length(p) == length(beta1)) {
which <- !p # p is indicator variable
}
else {
which <- rep(T, m)
which[p] <- F
}
beta <- beta1 * which
list(beta = beta, beta1 = beta1, x = x, p = p,
which = which)
}
#
#----------------------------------------------------------------
find.mu.alfa <- function(theta, theta1, x, ell, which,
...)
{
#
# find theta=c(beta,mu,alfa) which maximize
# Expected Log-Likelihood given by ell
#
# Set lower and upper bounds for mu (-Inf,Inf)
# and alfa (0,Inf)
#
lower <- c(which, T, 0) * theta
lower[c(which, T, 0) == T] <- - Inf
upper <- c(which, T, T) * theta
upper[c(which, T, T) == T] <- Inf #
optim(par=theta,ell,x = x,theta1 =theta1)
#nlminb(start = theta, objective = ell, x = x,
# lower = lower, upper = upper, theta1 # theta1, ...)
}
#
#----------------------------------------------------------------
ELL <- function(theta, theta1, x)
{
#
# Expected Log-Likelihood function for the Weibull regression model
# (see reference 1 in help file)
#
# Note: negative Log-likelihood value is returned
# to facilitate finding of extreme value of ELL in
# find.mu.alfa
#
np <- length(theta) - 2
alfa <- theta[np + 2]
mu <- theta[np + 1]
alfa1 <- theta1[np + 2]
mu1 <- theta1[np + 1]
beta <- theta[1:np]
beta1 <- theta1[1:np]
a <- alfa/alfa1
b.beta <- t(as.matrix(beta - a * beta1)) %*% t(x)
b <- (mu - a * mu1) + b.beta
ga1 <- gamma(a + 1) #
# negative value of ELL is returned !
#
- (log(alfa) - 0.57721566 * a + mean(b - exp(b) * ga1))
}
#
#
#-- koq ---------------------------------------------------------
#
params <- check.parameters(beta1 = beta1, x = x, p = p)
# are parameters OK?
beta1 <- params$beta1
beta <- params$beta
which <- params$which # which - coefficients not tested
np <- length(beta1)
x <- as.matrix(params$x)
if(length(beta1) != ncol(x))
stop("Number of coefficients in beta1 should be equal to the number of
variables (columns) in x"
)
if(verbose) {
cat("Beta = ", beta, "\n")
cat("Beta1 = ", beta1, "\n")
}
theta1 <- c(beta1, 0, 1)
theta <- c(beta, 0, 1)
b0 <- find.mu.alfa(theta = theta, theta1 = theta1, x = x,
ell = ELL, which = which)
theta0 <- b0$parameters #
GAMMA <- 2 * (ELL(theta = theta0, theta1 = theta1, x = x) -
ELL(theta = theta1, theta1 = theta1, x = x))
koq <- 1 - exp( - GAMMA)
koq
}
