Remove zero observations, adapt logLik

Zero observations at time zero made fitting the two-component error
model fail. A concentration of exactly zero does not make sense anyways,
as we generally have a limit of detection

2 files changed, 11 insertions, 24 deletions

diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
index d3f4232d..d812f177 100644
--- a/R/logLik.mkinfit.R
+++ b/R/logLik.mkinfit.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2018 Johannes Ranke
+# Copyright (C) 2018,2019 Johannes Ranke
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -16,29 +16,9 @@
 # You should have received a copy of the GNU General Public License along with
 # this program. If not, see <http://www.gnu.org/licenses/>
 logLik.mkinfit <- function(object, ...) {
-  y_ij <- object$data$observed
-  yhat_ij <- object$data$predicted
-  if (is.null(object$data$err)) {
-    # For unweighted fits we estimate a single value for sigma from the residuals
-    err <- sd(object$data$residual)
-    n_var_comp <- 1 # Number of variance components estimated
-  } else {
-    err <- object$data$err
-    # For weighted fits we check for variance models used in IRLS
-    # If the variance values (err) were given and were not
-    # reweighted, the number of variance components estimated is zero
-    if (is.null(object$reweight.method)) {
-      n_var_comp <- 0
-    } else {
-      n_var_comp <- switch(object$reweight.method,
-        obs = length(object$var_ms_unweighted),
-        tc = 2)
-    }
-  }
-  prob_dens <- dnorm(y_ij, yhat_ij, err)
-  val <- log(prod(prob_dens))
-  class(val) <- "logLik"
-  attr(val, "df") <- length(coef(object)) + n_var_comp
+  val <- object$logLik
+  # Number of estimated parameters
+  attr(val, "df") <- length(object$bparms.optim) + length(object$errparms)
   return(val)
 }
 # vim: set ts=2 sw=2 expandtab:
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 73c2f485..6c12d027 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -64,6 +64,12 @@ mkinfit <- function(mkinmod, observed,
   observed <- subset(observed, name %in% obs_vars)

   observed <- subset(observed, !is.na(value))

 

+  # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model)

+  if (any(observed$value == 0)) {

+    warning("Observations with value of zero were removed from the data")

+    observed <- subset(observed, value != 0)

+  }

+

   # Obtain data for decline from maximum mean value if requested

   if (from_max_mean) {

     # This is only used for simple decline models

@@ -405,6 +411,7 @@ mkinfit <- function(mkinmod, observed,
   } else {

     if(!quiet) cat("Optimisation successfully terminated.\n")

   }

+  fit$logLik <- - nlogLik.current

 

   # We need to return some more data for summary and plotting

   fit$solution_type <- solution_type