R-Lum/Luminescence

calc_MinDose: display gamma, sigma, p0 and mu in their real, absolute units

Closed this issue · 20 comments

SebastienHuot commented

Expected behaviour

gamma = 34.32

Observed behaviour

gamma = 3.54

Running mini example

library(Luminescence)

in the final table of values, but also with the various graphs

It is very hard to relate to the logarithm transportation

RLumSK commented

@tzerk Please have a look and decide. Thanks.

tzerk commented

I would be fine with that, though I do not see why this would be necessary. For the console output reporting the non-logged values would be straightforward, changing the likelihood profile plot will be a major problem. These are automatically generated by the bbmle package when calculating the profiles. I can only think of some really "hacky" ways to overplot the x-axis labels, which is likely prone to visual bugs and very likely also not worth the effort.

RLumSK commented

@tzerk I do agree, would you mind changing the console output accordingly? We do not touch the graphical output.

tzerk commented

Will do (when I find the time).

RLumSK commented

Sure, this is not at all urgent. Thank you!

SebastienHuot commented
RLumSK commented

While I do agree in general, we should keep in mind that what the function reflects is a translation of the model by Rex which expects gamma beeing the minimum dose on the log scale. The output should reflect this also.

@tzerk I double-checked how 'bbmle' produces the plot, the following example from the manual, slightly modified, should help:

library(bbmle)

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
d <- data.frame(x,y)
## we have a choice here: (1) don't impose boundaries on the parameters,
##  put up with warning messages about NaN values: 
fit1 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)),
             start=list(ymax=1,xhalf=1),
             data=d)
p1 <- suppressWarnings(profile(fit1))

##modify S4-object, kick out 2nd object 
p1@profile$xhalf <- NULL

plot(p1,main=c("first","second"),
     xlab=c(~y[max],~x[1/2]),ylab="Signed square root deviance",
     show.points=TRUE, xaxt = "n")

Now you have only one plot and you can set a double x-axis (the log and the unlogged) values.
It is a little hack, but not particularly complicated.

SebastienHuot commented
RLumSK commented

Don't worry, your suggestions are much appreciated, therefore I had given @tzerk the code snippets at hand to allow a straight forward integration.

tzerk commented

The calc_MinDose() function now always prints the parameter estimates and the confidence interval values in their absolute values, independent of whether the logged or the un-logged model was applied. I made the same changes to the the profile log-likelilhood graphs by directly overwriting the raw data in the bbmle objects with their exponential values (if applicable) before plotting them. I found this to be the most straightforward way of doing so, but I only did some fairly limited testing due to time constraints. @SebastienHuot please let me know whether it is now working as you'd expect or if there are any bugs.

  • 👍1
RLumSK commented

@tzerk Thank you very much! I double-checked the new code, however, do you see any possibility to print both, the log and the absolute values? If it takes too much of your time, we stick with the current implementation.

tzerk commented

I could certainly do that for the console output. For the plot however this is probably a bit more intricate, as this would require to forgo the automatic 'single plot' feature of the bbmle::plot() function and to do it manually with layout() or par() instead. This way it should be possible to add a second set of axis-labels to each individual plot before moving on to the next plot. I will see what I can do.

SebastienHuot commented

thanks Christoph!

displaying the results in their original (input) dimensional scale makes far better sense.

regarding Sebastian's suggestion, I will not re-iterate my arguments that triggered this issue. Only to suggest that presenting < log output > should be left as a optional argument, when calling calc_MinDose(). The default should be to present results in the same input dimensional scale.

Now, in regarding to testing. I have successfully installed the dev_0.9.x version

the various textual outputs appears ok.

likelihood plot for p0 is off-centered! (it is also the case with the 0.9.3 built)

for the sigma output result, it is properly returned, as seen in min@data[["summary"]][["sig"]]. However, in the console output (along with the likelihood profile plot for sigma, the value is exp(sigma)

Also, from inspecting the revised code, it seems you still require init.values in their log() transformation (for gamma and mu). It would be great to update this part.

thanks a lot for your time

Sebastien

RLumSK commented

@SebastienHuot Regarding p0, I am not a 100 % sure, but for me, it looks like a 'problem' in the plot method of 'bbmle'. If this is the case, I am not sure whether we can tackle this easily.

tzerk commented

The default should be to present results in the same input dimensional scale.

Either I am not understanding correctly, or this contradicts the initial request to have the output always in the absolute units for the output!? This is getting more and more confusing...

Only to suggest that presenting < log output > should be left as a optional argument, when calling calc_MinDose()

I've added the new argument log.output (FALSE by default). If TRUE, the logged parameter estimate and confidence level values are presented alongside the values in their absolute unit.

tzerk commented

Also, from inspecting the revised code, it seems you still require init.values in their log() transformation (for gamma and mu). It would be great to update this part.

I've change this in 7e73ed3. init.values must now always be in their absolute units. If the logged model is applied (log = TRUE) the values must still be given in their absolute units, but are then automatically log transformed before any calculations.

tzerk commented

Output after 65ddedc:

data(ExampleData.DeValues, envir = environment())
calc_MinDose(
  data = ExampleData.DeValues$CA1, 
  sigmab = 0.1, 
  par = 4,
  log = TRUE,
  log.output = TRUE # new argument, set to FALSE to hide all output with logged values
)

RLumSK commented

@tzerk Very nice!

tzerk commented

for the sigma output result, it is properly returned, as seen in min@data[["summary"]][["sig"]]. However, in the console output (along with the likelihood profile plot for sigma, the value is exp(sigma)

Isn't it the other way around? The min@data$summary data frame has all its output on the absolute scale, independent of whether the logged or non-logged model was applied. It is only for sigma that the value was either logged or non-logged, which is of course inconsistent. With 2a0dba1 it is now so that sigma is now also always given in its absolute value.

tzerk commented

@SebastienHuot Regarding p0, I am not a 100 % sure, but for me, it looks like a 'problem' in the plot method of 'bbmle'. If this is the case, I am not sure whether we can tackle this easily.

Yes, this is a problem in the bbmle package itself. I haven't had a look at their code for their custom plot method, but its logic for finding some proper xlim values sometimes fails. In 8c3fb24 I now override the xlim values with a super simple approach, but at least for the one tested data set it gets the job done. Do we want to stick with this approach? If yes, someone please test the code with some more data sets to see whether the plots are reasonable. If not, feel free to revert this commit.

Before

After