From 5814be02f286ce96d6cff8d698aea6844e4025f1 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 8 Apr 2019 18:00:25 +0200
Subject: 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
---
 R/logLik.mkinfit.R | 28 ++++------------------------
 R/mkinfit.R        |  7 +++++++
 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
-- 
cgit v1.2.1