[R] Error survreg: Density function returned an an invalid matrix
Israel Ortiz
isra4884 at gmail.com
Mon Nov 16 15:56:43 CET 2015
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. <therneau at mayo.edu>:
> 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-request at 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]]
More information about the R-help
mailing list