[R] koq.q ---- Kent O' Quigley R2

[Antonello Romani]

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
}