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 ([EMAIL PROTECTED]) or Janez Stare ([EMAIL PROTECTED]) 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(...)
"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 }
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html