aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2019-04-08 18:00:25 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2019-04-08 18:00:25 +0200
commit5814be02f286ce96d6cff8d698aea6844e4025f1 (patch)
treebdc45984af83e9629abc968e5dd5193ace2c3a95 /R
parent5e8190f26611c094a1a5d877a314cbca53e3e530 (diff)
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
Diffstat (limited to 'R')
-rw-r--r--R/logLik.mkinfit.R28
-rw-r--r--R/mkinfit.R7
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

Contact - Imprint