[R] Error survreg: Density function returned an an invalid matrix

[Israel Ortiz]

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]]