[R] r question
Rui Barradas
ruipbarradas at sapo.pt
Wed Mar 22 22:53:52 CET 2017
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: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 at 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.
>
More information about the R-help
mailing list