[R] r question

[謝孟珂]

Hi ,I have some question about simulate, I don't know how to paste question
to this website,so I paste below.
I use genoud to find the maximum likelihood value, but when I use numcut=3
,it will get error message,like this " coxph.wtest(fit$var[nabeta, nabeta],
temp, control$toler.chol) : NA/NaN/Inf"
and this is my code
library(rgenoud)
library(survival)
N=500
ei=rexp(N)
censor=rep(1,N)
x1=runif(N)
x2=runif(N)
x3=runif(N)
truecut=c(0.3,0.6)
dum1=1*(x1>truecut[1] & x1<truecut[2])
dum2=1*(x1>truecut[2])
x_true=cbind(dum1,dum2,x2,x3)
relativerisk<- matrix(log(c(1,1,2,2),ncol=4,byrow = TRUE)
beta_true<-relativerisk
Ti=exp(x_true%*%beta_true)*ei
confound=cbind(x2,x3)
initial2<-c(0.09,0.299,0.597,-0.17,-1.3,-3.1,-1.4,-1.12)
numcut=2
 domain2<<-matrix(c(rep(c(0.05,0.95),numcut),rep(c(0,5),numcut+dim(confound)[2])),ncol=2,byrow
= TRUE)

loglikfun <- function(beta, formula) {
  beta1 <- coxph(formula, init = beta,
control=list('iter.max'=0))#iteration is zero
  return(beta1$loglik[2])
}
obj <- function(xx){
  cutoff <- xx[1:numcut_global] #cutpoint
  cut_design <-
cut(target_global,breaks=c(0,sort(cutoff)+seq(0,gap_global*(length(cutoff)-1),by=gap_global),target_max),quantile=FALSE,labels=c(0:numcut_global))
  beta <- -xx[(numcut+1):nvars]  #coefficients of parameters
  logliks <-
loglikfun(beta,Surv(time_global,censor_global)~cut_design+confound_global)
  return(logliks)
}
maxloglik<-function(target,numcut,time,censor,confound,domain2,initial2,gap){
  time_global<<-time
  censor_global<<-censor
  target_global<<-target
  nvars<<-2*numcut+dim(confound)[2]
  confound_global<<-confound
  numcut_global<<-numcut
  target_max<<-max(target)
  gap_global<<-gap
  ccc<-genoud(obj, nvars, max=TRUE, pop.size=100, max.generations=6,
wait.generations=10,
              hard.generation.limit=TRUE, starting.values=initial2,
MemoryMatrix=TRUE,
              Domains=domain2, solution.tolerance=0.001,
              gr=NULL, boundary.enforcement=2, lexical=FALSE,
gradient.check=TRUE)
  ccc$par_corr<-ccc$par #the coefficients of genoud

ccc$par_corr[1:numcut]<-sort(ccc$par[1:numcut])+seq(0,gap_global*(numcut-1),by=gap_global)
#sort cutpoint
  return(ccc)
}


maxloglik(x1,3,Ti,censor,confound,domain2,initial2,0.02)$par_corr

I have no idea about the error ,maybe is my initial is wrong.
thank you
 From Meng-Ke

	[[alternative HTML version deleted]]