summary method does not show fixed parameters
Closed this issue · 9 comments
The summary method does not show fixed parameters, as the show and coef methods do:
> mag5 <- c(103, 115, 13, 2, 67, 36, 51, 8, 6, 61, 10, 21,
+ 7, 65, 4, 49, 92, 37, 16, 6, 23, 9, 2, 6, 5, 4,
+ 1, 3, 1, 9, 2)
> summary(fitbs(mag5))
Maximum likelihood estimation
Call:
mle2(minuslogl = function (N, S)
-sum(dbs(x = x, N = N, S = S, log = TRUE)), fixed = list(N = 834,
S = 31L), data = list(x = c(103, 115, 13, 2, 67, 36, 51,
8, 6, 61, 10, 21, 7, 65, 4, 49, 92, 37, 16, 6, 23, 9, 2, 6, 5,
4, 1, 3, 1, 9, 2)), eval.only = TRUE)
Coefficients:
Estimate Std. Error z value Pr(z)
-2 log L: 266.5852
> summary(fitls(mag5))
Maximum likelihood estimation
Call:
mle2(minuslogl = function (N, alpha)
-sum(dls(x, N, alpha, log = TRUE)), start = list(alpha = 6.3443128102516),
method = "Brent", fixed = list(N = 834), data = list(x = c(103,
115, 13, 2, 67, 36, 51, 8, 6, 61, 10, 21, 7, 65, 4, 49, 92,
37, 16, 6, 23, 9, 2, 6, 5, 4, 1, 3, 1, 9, 2)), lower = 0,
upper = 31L)
Coefficients:
Estimate Std. Error z value Pr(z)
alpha 6.3443 2.8323 2.24 0.02509 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-2 log L: 264.1686
> summary(fitvolkov(mag5))
Maximum likelihood estimation
Call:
mle2(minuslogl = function (theta, m, J)
-sum(dvolkov(x, theta = theta, m = m, J = J, log = TRUE)), start = list(
theta = 17.2803698475857, m = 0.0271718232953198), method = "L-BFGS-B",
fixed = list(J = 834), data = list(x = c(103, 115, 13, 2,
67, 36, 51, 8, 6, 61, 10, 21, 7, 65, 4, 49, 92, 37, 16, 6,
23, 9, 2, 6, 5, 4, 1, 3, 1, 9, 2)), lower = c(3.45607396951713,
1e-04), upper = c(86.4018492379283, 0.9999))
Coefficients:
Estimate Std. Error z value Pr(z)
theta 17.280383 13.033451 1.3258 0.1849
m 0.028630 0.029031 0.9862 0.3240
-2 log L: 259.1805 Raised by Giorgios Kokkoris (gkok@aegean.gr)
The problem here is that the function displaying these coefficients is `stats::printCoefmat``, which expects the coefficients to have a Std.error. I'm not sure where in the call chain we should do this editing: summary method, printCoefmat, or coef? Maybe we should just issue a warning in the summary?
More than that, I'm getting a feeling that there's something wrong in the way that we are dealing with the coef method:
> f <- fitlnorm(moths)
> coef(f)
meanlog sdlog
2.579059 1.778636
> printCoefmat(coef(f))
Error in printCoefmat(coef(f)) :
'x' must be coefficient matrix/data frameRelevant code from bbmle is here: https://github.com/cran/bbmle/blob/1e1ebb408e65d4e1f0d2466bc99cbee1111729a8/R/mle2-methods.R
I go for including in the summary a line with the fixed parameters.
My view is that summary methods are expected to return the relevant info about estimated coefficients and stats::printCoefmat is an utility function to ease this task. But all that coef is expected to do is to extract the coefficients from a model object. So I think it is ok that coef method for sads extracts both fixed and free parameters, as the user may need also fixed parameters (e.g. to use [dpq] probablity functions). Moreover, summary methods vary a lot among different models, in order to provide all relevant information. Given that the lack of fixed parameters in `summary' can be confusing for the user I think it does not hurt to include fixed parameters in the summary, but apart from the table with estimates and std errors of the free parameters.
Yes, that is certainly the right way of doing this.
Also, please ignore the error on 'x' must be coefficient matrix/data frame, I have figured out this is the right call:
> printCoefmat(coef(summary(f)))Ok, the problem here (as it already happened on a lot of our issues, such as the coef issues, the x@convergence issues, etc) is actually a problem with the bbmle package. In particular, this is a problem with the summary.mle2 method, as it does not record the fullcoef of the model:
> fitls(moths)->f
> summary(f)->sf
> class(sf)
[1] "summary.mle2"
attr(,"package")
[1] "bbmle"
> str(sf)
Formal class 'summary.mle2' [package "bbmle"] with 3 slots
..@ call : language mle2(minuslogl = function (N, alpha) -sum(dls(x, N, alpha, log = TRUE)), start = structure(list(alpha = 40.247281791951), .Names = "alpha"), method = "Brent", ...
..@ coef : num [1, 1:4] 4.02e+01 6.96 5.78 7.39e-09
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr "alpha"
.. .. ..$ : chr [1:4] "Estimate" "Std. Error" "z value" "Pr(z)"
..@ m2logL: num 2175We could reimplement the summary.mle2 class and generating method, or we can implement new summary.fitsad and summary.fitrad classes to deal with this. However, I feel that we are creating a large amount of "puxadinhos" to deal with problems that should be fixed in the bbmle package.
Implemented:
> summary(fitls(moths))
Maximum likelihood estimation
Call:
mle2(minuslogl = function (N, alpha)
-sum(dls(x, N, alpha, log = TRUE)), start = list(alpha = 40.247281791951),
method = "Brent", fixed = list(N = 15609L), data = list(x = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L,
8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 12L, 12L,
13L, 13L, 13L, 13L, 13L, 14L, 14L, 15L, 15L, 15L, 15L, 16L,
16L, 16L, 17L, 17L, 17L, 18L, 18L, 18L, 19L, 19L, 19L, 20L,
20L, 20L, 20L, 21L, 22L, 22L, 22L, 23L, 23L, 23L, 24L, 25L,
25L, 25L, 26L, 27L, 28L, 28L, 28L, 29L, 29L, 32L, 34L, 34L,
36L, 36L, 36L, 37L, 37L, 43L, 43L, 44L, 44L, 45L, 49L, 49L,
49L, 51L, 51L, 51L, 51L, 52L, 53L, 54L, 54L, 57L, 58L, 58L,
60L, 60L, 60L, 61L, 64L, 67L, 73L, 76L, 76L, 78L, 84L, 89L,
96L, 99L, 109L, 112L, 120L, 122L, 129L, 135L, 141L, 148L,
149L, 151L, 154L, 177L, 181L, 187L, 190L, 199L, 211L, 221L,
226L, 235L, 239L, 244L, 246L, 282L, 305L, 306L, 333L, 464L,
560L, 572L, 589L, 604L, 743L, 823L, 2349L)), lower = 0, upper = 240L)
Coefficients:
Estimate Std. Error z value Pr(z)
alpha 40.247 6.961 5.7818 7.391e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Fixed parameters:
N
15609
-2 log L: 2175.425 Looks nice. Where is the code?
Concerning your previous comment, in general I agree. But in this specific case the "puxadinho" is mostly harmless, or am I missing something?
Writing the docs for the summary.sads class before pushing to github ;)