Replacing negative values in a model (system of ODEs) with zero
My standard approach is to transform the state variables to an unconstrained scale. The most obvious/standard way to do this for positive variables is to write down equations for the dynamics of log(x) rather than of x.
For example, with the Susceptible-Infected-Recovered (SIR) model for infectious disease epidemics, where the equations are dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I we would naively write the gradient function as
gfun <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c(-beta*S*I,
beta*S*I-gamma*I,
gamma*I)
)
return(list(g))
}
If we make log(I) rather than I be the state variable (in principle we could do this with S as well, but in practice S is much less likely to approach the boundary), then we have d(log(I))/dt = (dI/dt)/I = beta*S-gamma; the rest of the equations need to use exp(logI) to refer to I. So:
gfun_log <- function(time, y, params) {
g <- with(as.list(c(y,params)),
c(-beta*S*exp(logI),
beta*S-gamma,
gamma*exp(logI))
)
return(list(g))
}
(it would be slightly more efficient to compute exp(logI) once and store/re-use it rather than computing it twice ...)
avoid negative values when resolving a ODE
I agree completely with @Lutz Lehmann, that the negative values are a result of the structure of the model.
The system of equations allows that derivatives still become negative, even if the states are already below zero, i.e. the states can further decrease. We don't have information about what the states are, so the following is only a technical demonstration. Here a dimensionless Monod-type feedback function fb is implemented as a safeguard. It is normally close to one. The km value should be small enough to act only for state values close to zero, and it should not be too small to avoid numerical errors. It can be formulated individually for each state. Other function types are also possible.
library(deSolve)
library(ggplot2)
library(reshape2)
values <- c(A = 1,
B = 1,
D = 1,
E = 20,
R = 1)
constants <- c(a = 1.2,
b = 0.5,
c = 1.2,
d = 1.5,
e = 0.3,
f = 0.5,
g = 1.5,
h = 0.9,
i = 1.3,
j = 1.3,
m = 0.8,
n = 0.6,
q = 1,
t = 0.0075,
u = 0.0009,
Pa = 100,
Pb = 0.05,
Pd = 0.1,
Pe = 10,
km = 0.001)
Dynamic_Model<-function(t, values, constants) {
with(as.list(c(values, constants)),{
fb <- function(x) x / (x+km) # feedback
dA <- (Pa + a*D - j*A - R) * fb(A)
dB <- (Pb + b*A + e*E - m*B) * fb(B)
dD <- (Pd + d*B + f*E - g*A - n*D) * fb(D)
dE <- (Pe - h*B + i*E - q*E) * fb(E)
dR <- (t*A*B - u*D*E) * fb(R)
list(c(dA, dB, dD, dE, dR))
})
}
times <- seq(0, 200, by = 0.1)
out <- ode(y = values, times = times, func = Dynamic_Model, parms = constants)
plot(out)
Additional hints:
- Removal of negative values afterwards (
out2 <- ifelse(out<0, 0, out)) is just wrong. - Removal of negative values in the model function, i.e.
use the ifelse in the main
would also be wrong as it can lead to a severe violation of mass balance.
- the time steps don't need to be very small. They are automatically adapted anyway by the solver. Too small time steps make your model slow and you get more outputs as needed.
- some of your parameters are quite large, so that the model becomes very stiff.
Replace negative values by zero
Thanks for the reproducible example. This is pretty basic R stuff. You can assign to selected elements of a vector (note an array has dimensions, and what you've given is a vector not an array):
> pred_precipitation[pred_precipitation<0] <- 0
> pred_precipitation
[1] 1.2091281 0.0000000 7.7665555 0.0000000 0.0000000 0.0000000 0.5151504 0.0000000 1.8281251
[10] 0.5098688 2.8370263 0.4895606 1.5152191 4.1740177 7.1527742 2.8992215 4.5322934 6.7180530
[19] 0.0000000 1.1914052 3.6152333 0.0000000 0.3778717 0.0000000 1.4940469
Benchmark wars!
@James has found an even faster method and left it in a comment. I upvoted him, if only because I know his victory will be short-lived.
First, I try compiling, but that doesn't seem to help anyone:
p <- rnorm(10000)
gsk3 <- function(x) { x[x<0] <- 0; x }
jmsigner <- function(x) ifelse(x<0, 0, x)
joshua <- function(x) pmin(x,0)
james <- function(x) (abs(x)+x)/2
library(compiler)
gsk3.c <- cmpfun(gsk3)
jmsigner.c <- cmpfun(jmsigner)
joshua.c <- cmpfun(joshua)
james.c <- cmpfun(james)
microbenchmark(joshua(p),joshua.c(p),gsk3(p),gsk3.c(p),jmsigner(p),james(p),jmsigner.c(p),james.c(p))
expr min lq median uq max
1 gsk3.c(p) 251.782 255.0515 266.8685 269.5205 457.998
2 gsk3(p) 256.262 261.6105 270.7340 281.3560 2940.486
3 james.c(p) 38.418 41.3770 43.3020 45.6160 132.342
4 james(p) 38.934 42.1965 43.5700 47.2085 4524.303
5 jmsigner.c(p) 2047.739 2145.9915 2198.6170 2291.8475 4879.418
6 jmsigner(p) 2047.502 2169.9555 2258.6225 2405.0730 5064.334
7 joshua.c(p) 237.008 244.3570 251.7375 265.2545 376.684
8 joshua(p) 237.545 244.8635 255.1690 271.9910 430.566
But wait! Dirk wrote this Rcpp thing. Can a complete C++ incompetent read his JSS paper, adapt his example, and write the fastest function of them all? Stay tuned, dear listeners.
library(inline)
cpp_if_src <- '
Rcpp::NumericVector xa(a);
int n_xa = xa.size();
for(int i=0; i < n_xa; i++) {
if(xa[i]<0) xa[i] = 0;
}
return xa;
'
cpp_if <- cxxfunction(signature(a="numeric"), cpp_if_src, plugin="Rcpp")
microbenchmark(joshua(p),joshua.c(p),gsk3(p),gsk3.c(p),jmsigner(p),james(p),jmsigner.c(p),james.c(p), cpp_if(p))
expr min lq median uq max
1 cpp_if(p) 8.233 10.4865 11.6000 12.4090 69.512
2 gsk3(p) 170.572 172.7975 175.0515 182.4035 2515.870
3 james(p) 37.074 39.6955 40.5720 42.1965 2396.758
4 jmsigner(p) 1110.313 1118.9445 1133.4725 1164.2305 65942.680
5 joshua(p) 237.135 240.1655 243.3990 250.3660 2597.429
That's affirmative, captain.
This modifies the input p even if you don't assign to it. If you want to avoid that behavior, you have to clone:
cpp_ifclone_src <- '
Rcpp::NumericVector xa(Rcpp::clone(a));
int n_xa = xa.size();
for(int i=0; i < n_xa; i++) {
if(xa[i]<0) xa[i] = 0;
}
return xa;
'
cpp_ifclone <- cxxfunction(signature(a="numeric"), cpp_ifclone_src, plugin="Rcpp")
Which unfortunately kills the speed advantage.
Very small populations/state variables to zero in differential equation solutions
Here an example. As said, very pragmatic.
library(package = "deSolve")
lv <- function(times, state, parms) {
with(as.list(c(state, parms)), {
dR <- 2*R - 0.5*R*C
dC <- 0.2*R*C - 0.6*C
return(list(c(dR, dC)))
})
}
time_vec <- seq(from = 0, to = 100, length.out = 1e4)
y_0 <- c(R = 50, C = 20)
out <- ode(y = y_0, times = time_vec, func = lv, parms = NULL, method = "lsoda")
min(out[,-1])
eps <- 1e-2
## event triggered if state variable <= eps
rootfun <- function (t, y, pars) {
return(y - eps)
}
## sets state variable = 0
eventfun <- function(t, y, pars) {
if (y[1] <= eps) y[1] <- 0
if (y[2] <= eps) y[2] <- 0
return(y)
}
out1 <- ode(y = y_0, times = time_vec, func = lv, parms = NULL, method = "lsoda",
rootfun = rootfun,
events = list(func = eventfun, root = TRUE))
plot(out, out1)
min(out1[,-1])
Ordinary differential equations (ODEs) - Is there any way to prevent negative values?
There are a few major issues here; to cut a long story short you have defined your model system incorrectly to use with method = "iteration", so no amount of light tinkering will get you sensible results. I'll get to this in the second part of my answer, but first, I'll answer your original question.
Using Events to Obtain Non-Negative State Values
You can force non-negative population sizes in deSolve using events. I suggest you read the deSolve documentation further, as there are many ways to trigger events, but we'll integrate one into your code that is triggered during each time step. Because it is time-based, it'll need some sort of maximum time to refer to, so I've created a maxtime value that you can also use in your for-loop. You should define this value somewhere accessible to your other functions; I simply put it before your mod1 function declaration.
# Here we declare the maximum time to which your system will evaluate
maxtime <- 2
# This is where we define your event function
# Add this directly above your call to ode()
posfun <- function(t, y, parms){
with(as.list(y), {
y[which(y<0)] <- 0
return(y)
})
}
# Here's your original call to ode(), with a small addition
# Notice that we added events; iteration is missing, more on that later
out <- ode(func=mod,
y=states,
times=seq(time_step, time_step + 1, by=1),
parms=input_parameters,
events=list(func = posfun, time = c(0:maxtime)))
out <- as.data.frame(out)
# Don't forget to add maxtime to your for-loop
for(time in 1:maxtime){
out_tab <- mod1(out_tab = out_tab, time_step = time)
## print(out_tab)
out_tab[time + 1, c("Nh")] <- Nh[time + 1]
out_tab[time + 1, c("Ih")] <- Ih[time + 1]
Looking at posfun, we see that it simply checks each of our state variables at each time step, and sets any negative values to zero. If we check our output we see it gives us non-negative population densities:
out_tab
Sv Ev Iv Nh Ih
1 10.0000 3.00000 2.000000 2 1
3 179.7493 16.67784 3.244288 1 0
4 416.2958 10.01499 6.133576 8 4
Great, right? Well, not quite. Unfortunately, there is currently no way to use events when method = iteration. Given what you've shared about your model so far, there are certainly arguments for modelling it in continuous time. Just because dispersal events are discretized does not necessarily mean that birth, death, and infection events must be discrete as well. It is important to conceptualize the different timescales at which these phenomena are occurring in real life, and ask yourself whether or not you are capturing them accurately with your model. However, this is already getting beyond the scope of StackOverflow, so on to the second part:
Coding discrete-time models into R with deSolve
Looking at the documentation for the iteration method in deSolve, we see that: "Method "iteration" is special in that here the function func should return the new value of the state variables rather than the rate of change."
You have clearly coded in a continuous-time model that returns the value of derivatives, and not a discrete-time model that returns the value of your state variables. Your birth component in dSv is a bare constant, so you are using an intrinsic birth rate; in a discrete-time model your birth component would be some constant number of offspring (almost always an integer) multiplied by the number of "parents" from the previous timestep.
Looking at your system of equations, we can see how this creates a problem: Instead of having each individual of Sv give birth to 100 individuals at each time step, you are setting Sv to 100. It is inevitable that Sv will plummet down past zero as you quickly accumulate a net loss.
An example discrete-time model would look something like:
# discrete-time host-parasite model
Parasite <- function(t, y, ks) {
P <- y[1]
H <- y[2]
f <- A * P / (ks + H)
Pnew <- H * (1 - exp(-f))
Hnew <- H * exp(rH * (1 - H) - f)
list (c(Pnew, Hnew))
}
Notice how we constantly refer to the "present" values of P and H in order to calculate their population dynamics for the next time step. I hope this has been helpful to you, and I wish you luck with your modelling adventures!
Replace dash with zero without affecting negative numbers
Just do
x <- c("-","-121512","123","-")
x[x == "-"] <- NA
x
#[1] NA "-121512" "123" NA
If you need a numeric vector instead of character wrap x in as.numeric().
If you want to replace all "-" in a dataframe we can use the same logic
df1 <- data.frame(x = c("-","-121512","123","-"),
y = c("-","-121512","123","-"),
z = c("-","A","B","-"), stringsAsFactors = FALSE)
df1[df1 == "-"] <- NA
If you want numeric columns if appropriate then you type.convert
df1[] <- lapply(df1, type.convert, as.is = TRUE)
str(df1)
'data.frame': 4 obs. of 3 variables:
$ x: int NA -121512 123 NA
$ y: int NA -121512 123 NA
$ z: chr NA "A" "B" NA
Solving log-transformed ODE system without overflow error
[…] is log-transform just not applicable in my case?
I don’t know where your transform went wrong, but it will certainly not achieve what you think it does. Log-transforming as a means to avoid negative values makes sense and works if and only if the following two conditions hold:
- If the value of a dynamical variable approaches zero (from above), its derivative also approaches zero (from above) in your model.
- Due to numerical noise, your derivative may turn negative though it actually isn’t.
Conversely, it is not necessary or doesn’t work in the following cases:
If Condition 1 fails because your derivative never approaches zero in your model, but is strictly positive, you have no problem to begin with, as your derivative should not become negative in any reasonable implementation of your model. (You might make it happen by implementing some spectacular numerical annihilation, but that’s quite a difficult feat to achieve and not what I would consider a reasonable implementation.)
If Condition 1 fails because your derivative becomes truly negative in your model, logarithms won’t save you, because the dynamics wants to push the derivative below zero and the logarithms cannot represent this. You usually get an overflow error due to the logarithms becoming extremely negative or the adaptive integration fails.
Even if Condition 1 applies, Condition 2 can usually be handled by avoiding numerical annihilations and similar when implementing your model.
Unless I am mistaken, your model falls into the first category. If M1 goes to zero, mdot_convert goes towards zero and thus M1dot = mdot_grow_now - mdot_convert is strictly positive, because mdot_grow_now is. M2dot is strictly positive anyway. Thus, you gain nothing from log-transforming. In fact, in the vast majority of cases, your dynamical variables will quickly increase.
With all that being said, some things you might want to look into are:
- Normalising your variables to be in the order of magnitude of 1.
- Stochastic differential equations.
Related Topics
Shiny Rcharts Multiple Chart Output
Match Two Columns with Two Other Columns
Connect R and Vertica Using Rodbc
Rstudio Calls Source() When Saving Script
How to Find the Package Name in R for a Specific Function
Ggplot Legend - Scale_Colour_Manual Not Working
Convert Factor to Date Class for Multiple Columns
Add Column to Data Frame Which Returns 1 If String Match a Certain Pattern
Annotate Ggplot2 Facets with Number of Observations Per Facet
R Ggplot Ordering Bars in "Barplot-Like " Plot
Remove Duplicates Column Combinations from a Dataframe in R
R: Matrix by Vector Multiplication
How to Retrieve the Client's Current Time and Time Zone When Using Shiny
Combine/Merge Columns While Avoiding Na