Robustifying the generalized gamma fit - Project notebook for the Invasion Heterogeneity project

I’m still a bit concerned about the generalized gamma, for two reasons: am I getting good enough starting values; and why do I get errors in the function evaluation.

For the first, my previous work on moments was not helpful. So I think that the only robust approach is to start from lots of random start values. However, that brings up the second issue of rather frequent failures. I think the issue is that flexsurv is not quite standard in how it deals with edge conditions. So let’s look at this a bit more.

First, generate the error:

library(flexsurv)
temp <- filter(disperseLer, ID == "73_0")
cens_data_tble <- cens_dispersal_data(temp, 7)
startgg <- start_params(cens_data_tble, "gengamma")

     dist      AIC       p1        p2
1   lnorm 1555.858 1.135456 0.8428757
2 weibull 1530.989 1.342346 4.6042573
3   gamma 1533.055 1.630594 0.3855333

startgg

$mu
[1] 1.526981

$sigma
[1] 0.7449643

$Q
[1] 1

fit <- try(fitdistcens(cens_data_tble, "gengamma", start = startgg))
summary(fit)

Fitting of the distribution ' gengamma ' By maximum likelihood on censored data 
Parameters
       estimate Std. Error
mu    1.5702093 0.09299152
sigma 0.7276694 0.04804765
Q     1.1178364 0.22615356
Loglikelihood:  -763.355   AIC:  1532.71   BIC:  1544.015 
Correlation matrix:
              mu      sigma          Q
mu     1.0000000 -0.7392010  0.8825485
sigma -0.7392010  1.0000000 -0.7009284
Q      0.8825485 -0.7009284  1.0000000

That actually works now – what I had to fix was some mistakes in translating from the limiting distribugtions in start_gengamma. It seems like it may work now; next step is to test this on all the data.

Old results

That threw an error, as expected. Now let’s go back to my experiments on vectorization and see what happens when I pass in the values to the d and p functions; we’ll ontrast that with a well-functioning distribution.

attach(startgg)
dgengamma(numeric(0), mu,sigma, Q)

numeric(0)

dlnorm(numeric(0))

numeric(0)

dgengamma(c(0, 1, Inf, NaN, -1), mu,sigma, Q)

[1]      NaN 1.238682      NaN      NaN 0.000000

dlnorm(c(0, 1, Inf, NaN, -1))

[1] 0.0000000 0.3989423 0.0000000       NaN 0.0000000

dgengamma(c(0, 1, NA), mu,sigma, Q)

[1]      NaN 1.238682       NA

dlnorm(c(0, 1, NA))

[1] 0.0000000 0.3989423        NA

pgengamma(numeric(0), mu,sigma, Q)

numeric(0)

plnorm(numeric(0))

numeric(0)

pgengamma(c(0, 1, Inf, NaN, -1), mu,sigma, Q)

[1] 0.0000000 0.3014678 1.0000000       NaN 0.0000000

plnorm(c(0, 1, Inf, NaN, -1))

[1] 0.0 0.5 1.0 NaN 0.0

pgengamma(c(0, 1, NA), mu,sigma, Q)

[1] 0.0000000 0.3014678        NA

plnorm(c(0, 1, NA))

[1] 0.0 0.5  NA

detach(startgg)

Things that need fixing in dgengamma: - Should return 0 instead of NaN when x = 0 - Should return 0 instead of NaN when x = Inf

pgengamma looks ok.

So let’s try a wrapper:

dmygengamma <- function(x, mu, sigma, Q, ...) {
  print((x))
  print(c(mu, sigma, Q))
  res <- flexsurv::dgengamma(x, mu, sigma, Q, ...)
  res[x == 0] <- 0
  res[!is.finite(x)] <- 0
  print((res))
  res
}
pmygengamma <- function(q, mu, sigma, Q, ...) {
  print((q))
  print(c(mu, sigma, Q))
  res <- flexsurv::pgengamma(q, mu, sigma, Q, ...)
  print((res))
  res
}
fit <- try(fitdistcens(cens_data_tble, "mygengamma", start = startgg))

numeric(0)
[1] 0.2944190 0.2171903 1.0000000
numeric(0)
[1]   0   1 Inf NaN  -1
[1] 0.2944190 0.2171903 1.0000000
[1] 0.0000000 0.9172299 0.0000000 0.0000000 0.0000000
[1]  0  1 NA
[1] 0.2944190 0.2171903 1.0000000
[1] 0.0000000 0.9172299 0.0000000
[1] 0 1
[1] -0.2944190 -0.2171903 -1.0000000
[1]  0 NA
[1] 0 1
numeric(0)
[1] 0.2944190 0.2171903 1.0000000
numeric(0)
[1]   0   1 Inf NaN  -1
[1] 0.2944190 0.2171903 1.0000000
[1] 0.0000000 0.2272483 1.0000000       NaN 0.0000000
[1]  0  1 NA
[1] 0.2944190 0.2171903 1.0000000
[1] 0.0000000 0.2272483        NA
[1] 0 1
[1] -0.2944190 -0.2171903 -1.0000000
[1] NA NA
[1] 0 1
numeric(0)
[1] 0.2944190 0.2171903 1.0000000
numeric(0)
numeric(0)
[1] 0.2944190 0.2171903 1.0000000
numeric(0)
numeric(0)
[1] 0.2944190 0.2171903 1.0000000
numeric(0)
  [1]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
 [24]  1  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2
 [47]  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
 [70]  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
 [93]  2  2  2  2  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
[116]  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
[139]  3  3  3  3  3  3  3  3  3  3  4  4  4  4  4  4  4  4  4  4  4  4  4
[162]  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
[185]  4  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5
[208]  5  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6
[231]  6  6  6  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
[254]  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
[277]  8  8  9  9  9  9  9  9  9  9  9  9  9  9  9  9 10 10 10 10 10 10 10
[300] 10 10 10 10 11 11 11 11 12 12 12 12 12 12 12 12 12 14 15 16 16
[1] 0.2944190 0.2171903 1.0000000
  [1] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
  [8] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [15] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [22] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [29] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [36] 0.2272483 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [43] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [50] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [57] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [64] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [71] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [78] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [85] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
 [92] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 1.0000000 1.0000000
 [99] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[106] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[113] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[120] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[127] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[134] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[141] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[148] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[155] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[162] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[169] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[176] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[183] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[190] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[197] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[204] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[211] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[218] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[225] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[232] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[239] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[246] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[253] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[260] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[267] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[274] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[281] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[288] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[295] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[302] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[309] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[316] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
  [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 [24]  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1
 [47]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
 [70]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
 [93]  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
[116]  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
[139]  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3  3  3  3  3  3
[162]  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3
[185]  3  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
[208]  4  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5
[231]  5  5  5  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6
[254]  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
[277]  7  7  8  8  8  8  8  8  8  8  8  8  8  8  8  8  9  9  9  9  9  9  9
[300]  9  9  9  9 10 10 10 10 11 11 11 11 11 11 11 11 11 13 14 15 15
[1] 0.2944190 0.2171903 1.0000000
  [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
  [8] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [15] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [22] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [29] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [36] 0.0000000 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [43] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [50] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [57] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [64] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [71] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [78] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [85] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483
 [92] 0.2272483 0.2272483 0.2272483 0.2272483 0.2272483 0.9981086 0.9981086
 [99] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[106] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[113] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[120] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[127] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[134] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[141] 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086 0.9981086
[148] 0.9981086 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[155] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[162] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[169] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[176] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[183] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[190] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[197] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[204] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[211] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[218] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[225] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[232] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[239] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[246] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[253] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[260] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[267] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[274] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[281] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[288] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[295] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[302] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[309] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
[316] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
<simpleError in optim(par = vstart, fn = fnobjcens, fix.arg = fix.arg, gr = gradient,     rcens = rcens, lcens = lcens, icens = icens, ncens = ncens,     ddistnam = ddistname, pdistnam = pdistname, hessian = TRUE,     method = meth, lower = lower, upper = upper, ...): function cannot be evaluated at initial parameters>

So there really doens’t seem to be a way to stop this. It does seem that, with the start values I’m using, the CDF seems to be plateauing after just a few centimeters. So maybe I’m getting the wrong translation from the limiting distributions?