1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
#
# The CAKE R modules are based on mkin
# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
# Authors: Rob Nelson, Richard Smith
# Tessella Project Reference: 6245
# This program 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/>.”
#
CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$parms)),
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0,
upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
plot = FALSE, quiet = FALSE, err = NULL, weight = "none",
scaleVar = FALSE, atol=1e-6, control=list(),...)
{
### This is a modification based on the "mkinfit" function.
### version 0.1 July 20
###
# This version has been modified to expect SFO parameterised as k and flow fractions
# Based on code in IRLSkinfit
NAind <-which(is.na(observed$value))
mod_vars <- names(mkinmod$diffs)
observed <- subset(observed, name %in% names(mkinmod$map))
ERR <- rep(1,nrow(observed))
observed <- cbind(observed,err=ERR)
obs_vars = unique(as.character(observed$name))
if (is.null(names(parms.ini)))
names(parms.ini) <- mkinmod$parms
mkindiff <- function(t, state, parms) {
time <- t
diffs <- vector()
for (box in mod_vars) {
diffname <- paste("d", box, sep = "_")
diffs[diffname] <- with(as.list(c(time, state, parms)),
eval(parse(text = mkinmod$diffs[[box]])))
}
return(list(c(diffs)))
}
if (is.null(names(state.ini)))
names(state.ini) <- mod_vars
parms.fixed <- parms.ini[fixed_parms]
optim_parms <- setdiff(names(parms.ini), fixed_parms)
parms.optim <- parms.ini[optim_parms]
state.ini.fixed <- state.ini[fixed_initials]
optim_initials <- setdiff(names(state.ini), fixed_initials)
state.ini.optim <- state.ini[optim_initials]
state.ini.optim.boxnames <- names(state.ini.optim)
if (length(state.ini.optim) > 0) {
names(state.ini.optim) <- paste(names(state.ini.optim),
"0", sep = "_")
}
costFunctions <- CakeInternalCostFunctions(mkinmod, state.ini.optim, state.ini.optim.boxnames,
state.ini.fixed, parms.fixed, observed, mkindiff, scaleVar, quiet, atol=atol)
############### Iteratively Reweighted Least Squares#############
## Start with no weighting
fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower,
upper = upper,control=control,...)
if(length(control)==0)
{
irls.control <- list(maxIter=10,tol=1e-05)
control <- list(irls.control=irls.control)
}else{
if(is.null(control$irls.control))
{
irls.control <- list(maxIter=10,tol=1e-05)
control <- list(irls.control=irls.control)
}
}
irls.control <- control$irls.control
maxIter <- irls.control$maxIter
tol <- irls.control$tol
####
niter <- 1
## insure one IRLS iteration
diffsigma <- 100
olderr <- rep(1,length(mod_vars))
while(diffsigma>tol & niter<=maxIter)
{
err <- sqrt(fit$var_ms_unweighted)
ERR <- err[as.character(observed$name)]
costFunctions$set.error(ERR)
diffsigma <- sum((err-olderr)^2)
cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
olderr <- err
fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
niter <- niter+1
### If not converged, reweight and fit
}
###########################################
fit$mkindiff <- mkindiff
fit$map <- mkinmod$map
fit$diffs <- mkinmod$diffs
out_predicted <- costFunctions$get.predicted()
# mkin_long_to_wide does not handle ragged data
fit$observed <- reshape(observed, direction="wide", timevar="name", idvar="time")
names(fit$observed) <- c("time", as.vector(unique(observed$name)))
predicted_long <- mkin_wide_to_long(out_predicted, time = "time")
fit$predicted <- out_predicted
fit$start <- data.frame(initial = c(state.ini.optim, parms.optim))
fit$start$type = c(rep("state", length(state.ini.optim)),
rep("deparm", length(parms.optim)))
fit$start$lower <- lower
fit$start$upper <- upper
fit$fixed <- data.frame(value = c(state.ini.fixed, parms.fixed))
fit$fixed$type = c(rep("state", length(state.ini.fixed)),
rep("deparm", length(parms.fixed)))
fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
parms.all = c(fit$par, parms.fixed)
fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed)
fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
for (obs_var in obs_vars) {
type = names(mkinmod$map[[obs_var]])[1]
fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
}
data <- merge(observed, predicted_long, by = c("time", "name"))
names(data) <- c("time", "variable", "observed","err-var", "predicted")
data$residual <- data$observed - data$predicted
data$variable <- ordered(data$variable, levels = obs_vars)
fit$data <- data[order(data$variable, data$time), ]
class(fit) <- c("CakeFit", "mkinfit", "modFit")
return(fit)
}
|
