Hi, thanks for your help.I run again my code
this is my cold.
run"maxloglik(x1,numcut,Ti,censor,confound,domain2,initial2,0.02)$par_corr"
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,
1.25,1.5,2,2,
1.5,2.5,2,2,
1.75,2.5,2,2,
2,3,2,2,
3,5,2,2,
4,7,2,2,
5,9,2,2,
6,11,2,2)),ncol=4,byrow = TRUE)
beta_true<-relativerisk[1,] #I wanna see the effect of different
relativerisk,this way I choose relativerisk[1,]
Ti=exp(x_true%*%beta_true)*ei
confound=cbind(x2,x3)
numcut=3
initial2<-c(0.1303100 ,0.3259150 ,0.6264413,0.47507966 ,-0.03886098
,-0.21062905 ,-0.31606071 ,-0.17620576)
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)
}
but I get error message like this
> maxloglik(x1,numcut,Ti,censor,confound,domain2,initial2,0.02)$par_corr
Thu Mar 23 10:38:27 2017
Domains:
5.000000e-02 <= X1 <= 9.500000e-01
5.000000e-02 <= X2 <= 9.500000e-01
5.000000e-02 <= X3 <= 9.500000e-01
0.000000e+00 <= X4 <= 5.000000e+00
0.000000e+00 <= X5 <= 5.000000e+00
0.000000e+00 <= X6 <= 5.000000e+00
0.000000e+00 <= X7 <= 5.000000e+00
0.000000e+00 <= X8 <= 5.000000e+00
Data Type: Floating Point
Operators (code number, name, population)
(1) Cloning........................... 15
(2) Uniform Mutation.................. 12
(3) Boundary Mutation................. 12
(4) Non-Uniform Mutation.............. 12
(5) Polytope Crossover................ 12
(6) Simple Crossover.................. 12
(7) Whole Non-Uniform Mutation........ 12
(8) Heuristic Crossover............... 12
(9) Local-Minimum Crossover........... 0
HARD Maximum Number of Generations: 6
Maximum Nonchanging Generations: 10
Population size : 100
Convergence Tolerance: 1.000000e-03
Using the BFGS Derivative Based Optimizer on the Best Individual Each
Generation.
Checking Gradients before Stopping.
Not Using Out of Bounds Individuals and Not Allowing Trespassing.
Maximization Problem.
GENERATION: 0 (initializing the population)
Fitness value... -2.579176e+03
mean............ -3.480611e+03
variance........ 1.187916e+05
#unique......... 100, #Total UniqueCount: 100
var 1:
best............ 1.303100e-01
mean............ 5.286724e-01
variance........ 7.963177e-02
var 2:
best............ 3.259150e-01
mean............ 4.996113e-01
variance........ 5.772924e-02
var 3:
best............ 6.264413e-01
mean............ 5.465332e-01
variance........ 6.837739e-02
var 4:
best............ 4.750797e-01
mean............ 2.524373e+00
variance........ 2.104958e+00
var 5:
best............ -3.886098e-02
mean............ 2.500733e+00
variance........ 1.892662e+00
var 6:
best............ -2.106291e-01
mean............ 2.578064e+00
variance........ 2.202828e+00
var 7:
best............ -3.160607e-01
mean............ 2.470152e+00
variance........ 2.159971e+00
var 8:
best............ -1.762058e-01
mean............ 2.550891e+00
variance........ 1.880531e+00
Show Traceback
Rerun with Debug
Error in coxph.wtest(fit$var[nabeta, nabeta], temp, control$toler.chol) :
NA/NaN/Inf in foreign function call (arg 3)
I don't know what happened.Thanks for your help.
Meng-Ke
2017-03-23 5:53 GMT+08:00 Rui Barradas <ruipbarra...@sapo.pt
<mailto:ruipbarra...@sapo.pt>>:
Hello,
There's a paenthesis missing in
> relativerisk<- matrix(log(c(1,1,2,2),ncol=4,byrow = TRUE)
+ beta_true<-relativerisk
Error: unexpected symbol in:
"relativerisk<- matrix(log(c(1,1,2,2),ncol=4,byrow = TRUE)
beta_true"
The correct instruction would be
relativerisk<- matrix(log(c(1,1,2,2)),ncol=4,byrow = TRUE)
And there's an error in your matrix multiply.
> Ti=exp(x_true%*%beta_true)*ei
Error in x_true %*% beta_true : non-conformable arguments
It can be corrected if you transpose beta_true.
Ti=exp(x_true %*% t(beta_true))*ei
At this point I've stoped running your code, I believe you must
revise it and try to see what's wrong instruction by instruction. Do
that and post again.
ALso, get rid of the <<-, use <-
If you do this, I'll explain the difference in the next answer to
your doubts.
Hope this helps,
Rui Barradas
Em 22-03-2017 08 <tel:22-03-2017%2008>:11, 謝孟珂 escreveu:
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 <mailto:R-help@r-project.org> mailing list
-- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
<https://stat.ethz.ch/mailman/listinfo/r-help>
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
<http://www.R-project.org/posting-guide.html>
and provide commented, minimal, self-contained, reproducible code.
