aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--R/logLik.mkinfit.R35
-rw-r--r--R/mkinfit.R2
-rw-r--r--man/logLik.mkinfit.Rd40
-rw-r--r--man/mkinfit.Rd5
-rw-r--r--tests/testthat/test_irls.R4
5 files changed, 82 insertions, 4 deletions
diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
new file mode 100644
index 00000000..c30cc099
--- /dev/null
+++ b/R/logLik.mkinfit.R
@@ -0,0 +1,35 @@
+# Copyright (C) 2018 Johannes Ranke
+# Contact: jranke@uni-bremen.de
+
+# This file is part of the R package mkin
+
+# mkin is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+
+# 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)) {
+ err <- sd(object$data$residual)
+ n_var_comp <- 1 # Number of variance components estimated
+ } else {
+ err <- object$data$err
+ if (object$reweight.method == "obs") n_var_comp <- length(object$var_ms_unweighted)
+ else n_var_comp <- 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
+ return(val)
+}
+# vim: set ts=2 sw=2 expandtab:
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 8c7549ad..b27f67b4 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -859,7 +859,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
invisible(x)
}
-# Fit the mean absolute deviance against the observed values,
+# Fit the median absolute deviation against the observed values,
# using the current error model for weighting
.fit_error_model_mad_obs <- function(tmp_res, tc, iteration) {
mad_agg <- aggregate(tmp_res$res.unweighted,
diff --git a/man/logLik.mkinfit.Rd b/man/logLik.mkinfit.Rd
new file mode 100644
index 00000000..2f3e58a1
--- /dev/null
+++ b/man/logLik.mkinfit.Rd
@@ -0,0 +1,40 @@
+\name{logLik.mkinfit}
+\alias{logLik.mkinfit}
+\title{
+ Calculated the log-likelihood of a fitted mkinfit object
+}
+\description{
+ This function simply calculates the product of the likelihood densities
+ calc
+}
+\usage{
+\method{logLik}{mkinfit}(object, ...)
+}
+\arguments{
+ \item{object}{
+ An object of class \code{\link{mkinfit}}.
+ }
+ \item{\dots}{
+ For compatibility with the generic method
+ }
+}
+\value{
+ An object of class \code{\link{logLik}} with the number of
+ estimated parameters (degradation model parameters plus variance
+ model parameters) as attribute.
+}
+\examples{
+ sfo_sfo <- mkinmod(
+ parent = mkinsub("SFO", to = "m1"),
+ m1 = mkinsub("SFO")
+ )
+ d_t <- FOCUS_2006_D
+ d_t[23:24, "value"] <- c(NA, NA) # can't cope with zero values at the moment
+ f_nw <- mkinfit(sfo_sfo, d_t, quiet = TRUE)
+ f_obs <- mkinfit(sfo_sfo, d_t, reweight.method = "obs", quiet = TRUE)
+ f_tc <- mkinfit(sfo_sfo, d_t, reweight.method = "tc", quiet = TRUE)
+ AIC(f_nw, f_obs, f_tc)
+}
+\author{
+ Johannes Ranke
+}
diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd
index 00d7eb47..228bab24 100644
--- a/man/mkinfit.Rd
+++ b/man/mkinfit.Rd
@@ -179,7 +179,7 @@ mkinfit(mkinmod, observed,
\item{weight}{
only if \code{err}=\code{NULL}: how to weight the residuals, one of "none",
"std", "mean", see details of \code{\link{modCost}}, or "tc" for the
- two component error model. The option "manual" is available for
+ two component error model. The option "manual" is available for
the case that \code{err}!=\code{NULL}, but it is not necessary to specify it.
}
\item{tc}{The two components of the error model as used for (initial)
@@ -245,6 +245,9 @@ mkinfit(mkinmod, observed,
Plotting methods \code{\link{plot.mkinfit}} and
\code{\link{mkinparplot}}.
+ Comparisons of models fitted to the same data can be made using \code{\link{AIC}}
+ by virtue of the method \code{\link{logLik.mkinfit}}.
+
Fitting of several models to several datasets in a single call to
\code{\link{mmkin}}.
}
diff --git a/tests/testthat/test_irls.R b/tests/testthat/test_irls.R
index 07b74078..c2bfc373 100644
--- a/tests/testthat/test_irls.R
+++ b/tests/testthat/test_irls.R
@@ -65,7 +65,7 @@ test_that("Reweighting method 'tc' works", {
cores = if (Sys.getenv("TRAVIS") != "") 1 else 15)
parms_2_10 <- apply(sapply(f_2_10, function(x) x$bparms.optim), 1, mean)
parm_errors_2_10 <- (parms_2_10 - parms_DFOP_optim) / parms_DFOP_optim
- expect_true(all(abs(parm_errors_2_10) < 0.2))
+ expect_true(all(abs(parm_errors_2_10) < 0.3))
f_2_10_tc <- mmkin("DFOP", d_2_10, reweight.method = "tc", quiet = TRUE,
cores = if (Sys.getenv("TRAVIS") != "") 1 else 15)
@@ -76,7 +76,7 @@ test_that("Reweighting method 'tc' works", {
tcf_2_10_tc <- apply(sapply(f_2_10_tc, function(x) x$tc_fitted), 1, mean, na.rm = TRUE)
tcf_2_10_error_model_errors <- (tcf_2_10_tc - c(0.5, 0.07)) / c(0.5, 0.07)
- expect_true(all(abs(tcf_2_10_error_model_errors) < 0.2))
+ expect_true(all(abs(tcf_2_10_error_model_errors) < 0.3))
f_tc_100_1 <- suppressWarnings(mkinfit(DFOP, d_100_1[[1]], reweight.method = "tc", quiet = TRUE))
parm_errors_100_1 <- (f_tc_100_1$bparms.optim - parms_DFOP_optim) / parms_DFOP_optim

Contact - Imprint