Re: [R] Error survreg: Density function returned an an invalid matrix
Thanks Terry, I use the following formula for density:
[image: f_X(x)= \begin{cases} \frac{\alpha
x_\mathrm{m}^\alpha}{x^{\alpha+1}} & x \ge x_\mathrm{m}, \\ 0 & x <
x_\mathrm{m}. \end{cases}]
Where *x*m is the minimum value for x. I get this fórmula in
https://en.wikipedia.org/wiki/Pareto_distribution but there are a lot of
books and sites that use the same fórmula. This part of the code use that
formula:
distribution <- function(x, alpha) ifelse(x > min(x) ,
alpha*min(x)**alpha/(x**(alpha+1)), 0)
Also, I support my sintax in the following post:
http://stats.stackexchange.com/questions/78168/how-to-know-if-my-data-fits-pareto-distribution
Another option is transform my variable for time from pareto to exponential
(but this solution it's not very elegant):
If X is pareto distributed then
[image: Y = \log\left(\frac{X}{x_\mathrm{m}}\right)]
it's exponential distributed.
The syntax:
library(foreign)
library(survival)
library(VGAM)
set.seed(3)
X <- rpareto(n=100, scale = 5,shape = 1)
Y <- log(X/min(X))
hist(X,breaks=100)
hist(Y,breaks=100)
b <- rnorm(100,5,1)
c <- rep(1,100)
base <- cbind.data.frame(X,Y,b,c)
mod1<-survreg(Surv(Y+1, c) ~ b, base, dist = "exponential")# +1 it's
because time should be > 1
summary(mod1)
This solution works but I don´t like it.
Thanks.
2015-11-16 7:40 GMT-06:00 Therneau, Terry M., Ph.D. :
> The error message states that there is an invalid value for the density.
> A long stretch of code is not very helpful in understanding this. What we
> need are the definition of your density -- as it would be written in a
> textbook. This formula needs to give a valid response for the range
> -infinity to +infinity. Or more precisely, for any value that the
> maximizer might guess at some point during the iteration.
>
> Terry T.
>
>
> On 11/14/2015 05:00 AM, r-help-requ...@r-project.org wrote:
>
>> Thanks Terry but the error persists. See:
>>
>> >library(foreign)> library(survival)> library(VGAM) > mypareto <-
>>> list(name='Pareto',+ init=
>>>
>>
> remainder of message trucated
>
[[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.
Re: [R] Error survreg: Density function returned an an invalid matrix
The error message states that there is an invalid value for the density. A long stretch of code is not very helpful in understanding this. What we need are the definition of your density -- as it would be written in a textbook. This formula needs to give a valid response for the range -infinity to +infinity. Or more precisely, for any value that the maximizer might guess at some point during the iteration. Terry T. On 11/14/2015 05:00 AM, r-help-requ...@r-project.org wrote: Thanks Terry but the error persists. See: >library(foreign)> library(survival)> library(VGAM) > mypareto <- list(name='Pareto',+ init= remainder of message trucated __ 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.
Re: [R] Error survreg: Density function returned an an invalid matrix
Thanks Terry but the error persists. See:
> library(foreign)> library(survival)> library(VGAM) > mypareto <-
> list(name='Pareto',+ init= function(x, weights,parms){+
> alpha <- length(x)/(sum(log(x)))#this is a MLE for alpha+
> c(media <-(alpha*1/(alpha-1)),varianza <-
> ((1/alpha)^2)*(alpha/(alpha-2)))},+ density=
> function(x,weights) {+alpha <- length(x)/(sum(log(x)))+
> cdf1 <- function(x, alpha) ifelse(x > 1 , 1 - (1/x)**alpha,
> 0 )+cdf2 <- function(x, alpha) ifelse(x > 1, (1/x)**alpha
> ,0)+distribution <- function(x, alpha) ifelse(x > 1 ,
> alpha/(x**(alpha+1)), 0)+firstdev <- function(x, alpha)
> ifelse(x > 1, -(alpha+x)/x, 0)+seconddev <- function(x,
> alpha) ifelse(x > 1, (alpha+1)*(alpha+2)/x^2,0)+
> cbind(cdf1(x,alpha),cdf2(x, alpha),
> distribution(x,alpha),firstdev(x,alpha),seconddev(x,alpha))},+
> devian!
ce=function(x) {stop('deviance residuals not defined')},+
quantile= function(p, alpha) ifelse(p < 0 | p > 1, NaN, 1*(1-p)**(-1/alpha)))>
> survregDtest(mypareto, TRUE)[1] TRUE> set.seed(1)> a <- rpareto(100, 1, 1) >
b <- rnorm(100,5,1)> c <- rep(1,100)> base <- cbind.data.frame(a,b,c)>
mod1<-survreg(Surv(a, c) ~ b, base, dist = mypareto)Error in survreg.fit(X, Y,
weights, offset, init = init, controlvals = control, :
Density function returned an invalid matrix
2015-11-04 7:52 GMT-06:00 Therneau, Terry M., Ph.D. :
> Hi, I want to perform a survival analysis using survreg procedure from
>> survival library in R for a pareto distribution for a time variable, so I
>> set the new distribution using the following sintax:
>>
>> library(foreign)
>> library(survival)
>> library(VGAM)
>>
>> mypareto <- list(name='Pareto',
>> init= function(x, weights,parms){
>>
> etc.
>
> The survreg routine fits location-scale distributions such that (t(y) -
> Xb)/s ~ F, where t is an optional transformation, F is some fixed
> distribution and X is a matrix of covariates. For any distribution the
> questions to ask before trying to add the distribution to survreg are
> - can it be written in a location-scale form?
> - if so, how do the parameters of the distribution map to the location
> (Xb) and scale (s).
>
> In fitting data we normally have per-subject location (X b) but an
> intercept-only model is of course possible.
>
> If y is Weibull then log(y) fits into the framework, which is how survreg
> fits it. The transformation of parameters location and scale parameters
> for log(y) back to the usual Weibull parameterization for y often trips
> people up (see comments in the Examples section of ?survreg).
>
> The log of a Pareto can be written in this form (I think?). The two
> parameters are the scale a and lower limit b, with survival function of
> S(x)= (b/x)^a, for x >= b. If y = log(x) the survival function for y is
> S(y) = (b/exp(y))^a = exp[-(y - log(b))/(1/a)], which has location log(b)
> and scale 1/a. But even if I am correct the discontinuity at b will cause
> the underlying Newton-Raphson method to fail.
>
> Terry Therneau
>
[[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.
Re: [R] Error survreg: Density function returned an an invalid matrix
Hi, I want to perform a survival analysis using survreg procedure from
survival library in R for a pareto distribution for a time variable, so I
set the new distribution using the following sintax:
library(foreign)
library(survival)
library(VGAM)
mypareto <- list(name='Pareto',
init= function(x, weights,parms){
etc.
The survreg routine fits location-scale distributions such that (t(y) - Xb)/s ~ F, where t
is an optional transformation, F is some fixed distribution and X is a matrix of
covariates. For any distribution the questions to ask before trying to add the
distribution to survreg are
- can it be written in a location-scale form?
- if so, how do the parameters of the distribution map to the location (Xb)
and scale (s).
In fitting data we normally have per-subject location (X b) but an intercept-only model is
of course possible.
If y is Weibull then log(y) fits into the framework, which is how survreg fits it. The
transformation of parameters location and scale parameters for log(y) back to the usual
Weibull parameterization for y often trips people up (see comments in the Examples section
of ?survreg).
The log of a Pareto can be written in this form (I think?). The two parameters are the
scale a and lower limit b, with survival function of S(x)= (b/x)^a, for x >= b. If y =
log(x) the survival function for y is
S(y) = (b/exp(y))^a = exp[-(y - log(b))/(1/a)], which has location log(b) and scale 1/a.
But even if I am correct the discontinuity at b will cause the underlying Newton-Raphson
method to fail.
Terry Therneau
__
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.
