[R] Gradient problem in nlm
Hello everyone!
I am having some trouble supplying the gradient function to nlm in R for
windows version 2.2.1.
What follows are the R-code I use:
fredcs39<-function(a1,b1,b2,x){return(a1+exp(b1+b2*x))}
loglikcs39<-function(theta,len){
value<-sum(mcs39[1:len]*fredcs39(theta[1],theta[2],theta[3],c(8:(7+len))) -
pcs39[1:len] * log(fredcs39(theta[1],theta[2],theta[3],c(8:(7+len)
a1<-theta[1]; b1<-theta[2]; b2<-theta[3]
df.a1<-sum(-mcs39[1:len] + pcs39[1:len]/(a1+exp(b1+b2*c(8:(7+len)
df.b1<-sum( -mcs39[1:len] * exp(b1+b2*c(8:(7+len))) + (pcs39[1:len] *
exp(b1+b2*c(8:(7+len))) ) /(a1+exp(b1+b2*c(8:(7+len)
df.b2<- sum(-mcs39[1:len] * exp(b1+b2*c(8:(7+len))) * c(8:(7+len)) +
(pcs39[1:len] * exp(b1+b2*c(8:(7+len))) * c(8:(7+len)) )
/(a1+exp(b1+b2*c(8:(7+len )
attr(value,"gradient")<-c(df.a1,df.b1,df.b2)
return(value)
}
theta.start<-c(0.01 ,-1.20, -0.0005)
outcs39<-nlm(loglikcs39,theta.start,len=50)
Error in nlm(loglikcs39, theta.start, len = 50) :
probable coding error in analytic gradient
Any light that can be shed on this would be highly appreciated.
Many thanks
Singyee Ling
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
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.
[R] warning message in nlm
Dear R-users,
I am trying to find the MLEs for a loglikelihood function (loglikcs39) and
tried using both optim and nlm.
fredcs39<-function(b1,b2,x){return(exp(b1+b2*x))}
loglikcs39<-function(theta,len){
sum(mcs39[1:len]*fredcs39(theta[1],theta[2],c(8:(7+len))) - pcs39[1:len] *
log(fredcs39(theta[1],theta[2],c(8:(7+len)
}
theta.start<-c(0.1,0.1)
1. The output from using optim is as follow
--
optcs39<-optim(theta.start,loglikcs39,len=120,method="BFGS")
> optcs39
$par
[1] -1.27795226 -0.03626846
$value
[1] 7470.551
$counts
function gradient
133 23
$convergence
[1] 0
$message
NULL
2. The output from using nlm is as follow
---
> outcs39<-nlm(loglikcs39,theta.start,len=120)
Warning messages:
1: NA/Inf replaced by maximum positive value
2: NA/Inf replaced by maximum positive value
3: NA/Inf replaced by maximum positive value
4: NA/Inf replaced by maximum positive value
5: NA/Inf replaced by maximum positive value
6: NA/Inf replaced by maximum positive value
7: NA/Inf replaced by maximum positive value
> outcs39
$minimum
[1] 7470.551
$estimate
[1] -1.27817854 -0.03626027
$gradient
[1] -8.933577e-06 -1.460512e-04
$code
[1] 1
$iterations
[1] 40
As you can see, the values obtained from using both functions are very
similar. But, what puzzled is the warning message that i got from using
nlm. Could anyone please shed some light on how this warning message come
about and whether it is a cause for concern?
Many thanks in advance for any advice!
singyee
[[alternative HTML version deleted]]
__
R-help@stat.math.ethz.ch mailing list
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.
[R] question about survSplit
Dear R-users, I use the survsplit function in the survival package to change my data into counting-process format and the transformed format is as follow: (a) start stop event DP age 0 5 01 20 5 1001 20 1025 11 20 looking at the above three entries that belong to the same person, if an event happen at time 5, won't the person actually enter the risk set twice since there is another entry that start at time 5 and Cox proportional hazard model "won't know" that it actually belong to the same person. Shouldn't it be like this? (b) start stop event DP age 0 5 01 20 610 01 20 11 2511 20 or the R-function coxph has already take this into account when calculating its risk set although I believe (a) is actually the correct one. Your advice is greatly appreciated! regards, sing yee [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] left-truncation in survreg
Dear R-users, I know that a few people have asked whether survreg handles left-truncation data and the reply that i have seen so far is that it does. However, when I try to use survreg on left-truncated data, I got the following error message. > survcs3<-survreg(Surv(start,end,status)~AG, data=DPONEcs3, dist="exponential") Error in survreg(Surv(start, end, status) ~ AG, data = DPONEcs3, dist = "exponential") : Invalid survival type Then i tried > cc<-survreg(Surv(end,status)~AG, data=DPONEcs3, dist="exponential") and it works perfectly fine. May I know why survreg does not accept my left-truncated data? I used similar set of data on coxph and it works fine Many many thanks for your time and any comments given kind regards, sing yee [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] andersen plot vs score process or scaled Schoenfeld residuals to test for proporti0nal hazards
Dear all, I use the Andersen plot to check for proportional hazards assumption for a factor (say x) in the Cox regression model and obtained a straight line that pass through the origin. However, the formal test done by the R-function cox.zph, which is based on the plot of Schonefeld residuals against time, indicates that proportional hazards assumption is violated. Further, a plot of the score process (cumulative sums of schoenfeld residuals) against time again give the same conclusion as the cox.zph function and i am really stumped by this. Klein et al (1997, pp 354) mentioned that graphical checks for proportional hazards assumption is often preferred to formal test as formal test based on a large enough sample (in my case is 4000 data entries), will often reject the null hypothesis of proportionality. Is that what is happening in my case? any suggestion? Thanks! kind regards, singyee ling [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] breslow estimator for cumulative hazard function
Dear R-users, I am checking the proportional hazard assumption of a cox model for a given covariate, let say Z1, after adjusting for other relavent covariates in the model. To this end, I fitted cox model stratified on the discrete values of Z1 and try to get beslow estimator for the baseline cumulative hazard function (H(t)) in each stratum. As far as i know, if the proportionality assumption holds, the plot of ln[H(t)] of each stratum versus time should be approximately parallel. i.e fit<-coxph(Surv(start,end,status)~sx+rated+AGLEVEL+strata(Z1),data=ALLDPinfectionandbronchitis) ss<-survfit(fit) plot(ss,fun="cumhaz") My question is on whether the cumulative hazard given by the above command is actually a breslow estimator for baseline cumulative hazard ,i.e, estimator=sum( number of death/ (sum(risk score in risk set)) or a nelson-Aalen estimator. if the above command does not give me breslow estimator, please advise on how I can get it. Thanks for any help given. kind regards, sing yee [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] question about plotting survfit object on log-log scale
Dear all, I tried to plot cumulative hazard H(t) against t on a log-log sclale by putting the optional function "fun ='cloglog'" in the plot function. However, I came across the following error message > survcs20fit2<-survfit(Surv(start,end,status)~sx,data=ALLDPcs20) >plot(survcs20fit2,fun="cloglog") Error in rep.default (2, n2 - 1) : invalid number of copies in rep ( ) In addition : Warning message: 2 x values <= 0 omitted from logarithmic plot in : xy.coords ( x, y, xlabel, ylabel, log) What does the error message :" Error in rep.default (2, n2 - 1) : invalid number of copies in rep ( )" actually means and what can I do to rectify it. Thanks for any help given! kind regards, sing yee [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] cox.zph
Dear R-users, I am sorry if this is obvious. I am testing the proportional hazard assumptions using cox.zph. If i am not wrong, a g(t) function must be assumed. Four possibilities available in R are "km","identity" and "rank". may i know what functions of time are these transformation assuming? Thanks a lot in advance for your wisdom. kind regards, sing yee ling [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] questions about coxph.detail
Dear R- users, I obtained a coxph object using the coxph function and use the coxph.detail function to understand what coxph function actually does. what follow are the R output fit<-coxph(Surv(B3TODEATH,STATUS)~GROUP+AGE+SEX+SEX:AGE) detail<- coxph.detail(fit) <-detail $time .(rows of events time) $means ..(a matrix consisting of mean of each variable at event time) $nevent ..( number of event at unique event time) $nrisk ..(number of ppl at risk at each event time) $hazard ...(hazard at each event time) etc...etc okay, my question is on how the hazard at each event time is calculated. I suspect that it is h(t) = (nevent(t)/nrisk(t))* exp( sum(coeff* means(t)). Am I right? if not, how is the hazard actually calculated? [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] rescale x-axis
Dear all, I am trying to draw a survival curve with probability of surviving as the y-axis and days (0- 500 days )as the x-axis. however, i do not want the days to be equally spaced on the x-axis as i am more interested in looking at the behaviour of the curve in the first 50 days. I am reluctant to use xlim=c(0,1000) as i want to see the whole picture. Hence, what I am interested in is a scale in which the days are not equally spaced. By that , I mean the length of the interval between the days get smaller and smaller, which gives greater emphasis to the intial period. (i.e the length of the interval betwen 0-1 days is longer then the interval between 1-2 days and so on) .Hope what i say above make sense. any advise? thanks! sing yee __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
