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]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
