Comments (9)
The workflow as of pomp version 2 is to construct an objective function that is then optimized using any general-purpose optimization tool. The optimization function created is stateful in that it remembers the point at which it was last evaluated.
For example, the following fits four parameters from an SIR model to simulated data (the model and data are those from the example included in the package).
library(pomp)
enames <- c("gamma","iota","rho","k")
sir() |>
traj_objfun(
dmeasure=Csnippet("
lik = dnbinom_mu(nearbyint(reports),1/k,rho*cases,give_log);
"),
rmeasure=NULL,
params=c(coef(sir()),k=1),
partrans=parameter_trans(
log=c("gamma","iota","k"),
logit="rho"
),
est=enames,
paramnames=enames,
statenames="cases"
) -> ofun
theta <- c(gamma=10,iota=1,rho=0.2,k=1)
library(subplex)
subplex(
fn=ofun,
par=partrans(as(ofun,"pomp"),theta,dir="toEst"),
control=list(trace=0)
) -> fit
fit
ofun(fit$par)
coef(ofun)
library(ggplot2)
ofun |>
trajectory(format="d") |>
ggplot(aes(x=time,y=coef(ofun,"rho")*cases))+
geom_line()+
geom_point(
data=as(ofun,"data.frame"),
mapping=aes(x=time,y=reports)
)+
theme_bw()
The chunk beginning with sir()
constructs the pomp object and pipes it to traj_objfun
, which constructs the objective function. Note that I have here modified the measurement model from the example that comes with the package. Parameter transformations are provided for each of the four parameters which are to be estimated. theta
is a (not very good) guess as to the unknown parameters, which will be used as a starting place for the search. The next chunk (beginning with library(subplex)
), uses subplex
to maximize the likelihood via trajectory-matching. The result is stored in fit
. I then examine fit
, call ofun
one last time at the estimated parameters (just in case), and then plot the data and the expected value of the data under the fitted model on the same axes.
from pomp.
I've tried to streamline your code down to focus on the essentials. Since you are doing trajectory matching, your rproc
is neither here nor there. I don't have dat3
, so I can't verify if this works properly. If the thread continues, please upload that (or something that acts like it; it isn't necessary to include the actual data) with your reply.
library(pomp)
model <- vectorfield(
Csnippet("
DVs = (beta-c*T-alpha)*Vs;
DVb = alpha*Vs+(beta-c*T)*Vb;
DT = (sigma*(Vs+Vb)-m)*T+p;"
)
)
init <- Csnippet("
Vs_r = 200;
Vb_r = 0;
background = 3000;
Vs = 3020;
Vb = 3000;
T = 20;
")
rmeas = Csnippet("
N_Vs = rpois(Vs*rho);
N_Vb = rpois(Vb*rho);")
dmeas <- Csnippet("
lik = dpois(N_Vs,Vs*rho,1)+
dpois(N_Vb,Vb*rho,1);
lik = (give_log) ? lik : exp(lik)")
Note that the dmeas
as you had written it was incorrect. In the case that give_log=1
it would return the correct answer (the sum of the log likelihoods), but if give_log=0
, it would return nonsense (the sum of the likelihoods).
dat3 |>
traj_objfun(
times="time", t0=0,
skeleton=model,
rmeasure=rmeas,
dmeasure=dmeas,
rinit = init,
partrans=parameter__trans(log=c("beta","c","alpha","sigma","m","p")),
paramnames=c("beta","c","alpha","sigma","m","p","rho","shape_fit","rate_fit","factor","area_sg","area_b"),
statenames=c("Vs_r","Vb_r","Vs","Vb","T","background"),
est=enames,
params = c(
beta=0.7,c = 0.01,alpha = 0.01,
sigma = 0.0002, m = 0.05, p = 1,rho = 0.9,
factor = factor,area_sg = area_sg, area_b = area_b,
shape_fit = shape_fit, rate_fit = rate_fit
)
) -> ofun
I am not sure where the warnings are coming from. You are postulating a Poisson measurement model: are the data all integers, as this model assumes? If not, this may be your problem. If they are, then perhaps you can include the data (dat3
or equivalent), and we can do some more detailed debugging.
from pomp.
One thing does occur to me: you define Vs_r
, Vb_r
, and background
as latent states, but you do not provide any derivatives. I suppose you are assuming that, left unspecified, these derivatives are effectively zero. This is not a safe assumption. In general, uninitialized memory locations contain undefined values. I would suggest either:
model <- vectorfield(
Csnippet("
DVs = (beta-c*T-alpha)*Vs;
DVb = alpha*Vs+(beta-c*T)*Vb;
DT = (sigma*(Vs+Vb)-m)*T+p;
DVs_r = DVb_r = Dbackground = 0;"
)
)
which ensures that these state variables do not change their values. More elegantly (and efficiently), you could make these three quantities into model parameters (or globals, using the globals
argument (see the documentation).
from pomp.
Your responses have been very helpful! I have the model working! Another issue has arisen though. So far I have been fitting to data capturing Vs and Vb over time (measured daily). However, I also have data capturing T, but this is not measured as often (every 4 days). Is there a way to include all the data when doing the fitting?
I tried to follow what was written in https://kingaa.github.io/pomp/FAQ.html under the heading "3.3 How do I deal with missing data?", but this doesn't seem to work. This is how I rewrote dmeasure:
dmeas <- Csnippet("
if (ISNA(N_T)) {
lik = dpois(N_Vs,Vs*rho1+1e-6,1)+
dpois(N_Vb,Vb*rho1+1e-6,1);
lik = (give_log) ? lik : exp(lik);
} else if(ISNA(N_Vb)){
lik = dpois(N_T,T*rho2+1e-6,1);
lik = (give_log) ? lik : exp(lik);
} else {
lik = dpois(N_Vs,Vs*rho1+1e-6,1)+
dpois(N_Vb,Vb*rho1+1e-6,1)+
dpois(N_T,T*rho2+1e-6,1);
lik = (give_log) ? lik : exp(lik);
}
")
It seems like since there are NAs in some of the N_T data, it ignores N_T entirely.
I've attached an example of what the data looks like.
from pomp.
I don't immediately see anything wrong with your code. Do you see error messages? On what basis do you conclude that N_T
is being entirely ignored?
from pomp.
I also note that your data variables N_Vs
and N_Vb
are not integers, as your model assumes. What are these variables? Why have you selected a Poisson measurement model?
from pomp.
Yes, you are right that N_Vs and N_Vb are not integers. We are infecting mice with a luciferase-tagged virus that emits light, and putting these mice in a bioimager. The values represent the amount of light emitted by the virus in units of photons/s. For simplicity, I'm assuming that each virus emits 1 photon/s. There is also some background signal in this measurement which I am also trying to account for. Perhaps a gamma distribution would be better. I've definitely been having some issues with writing an appropriate dmeasure and rmeasure.
N_T represents CD8 T cell amounts.
I think the issue I was having with fitting is not only due to perhaps a poorly written dmeasure expression, but also due to the fact that the values of N_Vs/N_Vb and N_T are on vastly different scales. Because the viral load measures are so large, it leads to large errors between the data and model-predicted values that outweigh any error related to T cells. As a result, the model fits the viral data well, and the T cell data poorly. Could this be what's happening?
from pomp.
The Poisson is a discrete distribution: It is defined only on the non-negative integers. I think it would be wise were you to give some more thought to the measurement model. For example: what are the errors you expect to see on these measurements? Do you think the error is better conceived of as additive or multiplicative?
I think the issue I was having with fitting is not only due to perhaps a poorly written dmeasure expression, but also due to the fact that the values of N_Vs/N_Vb and N_T are on vastly different scales. Because the viral load measures are so large, it leads to large errors between the data and model-predicted values that outweigh any error related to T cells. As a result, the model fits the viral data well, and the T cell data poorly.
The weight given to each distinct data stream is a function of the distribution of errors. If you are assuming that the error on the fluorescence values is small or large, then these will be given heavy or light weights, respectively. In other words, it is not the magnitudes per se but the scaling of the error relative to this magnitude (e.g., the coefficient of variation) that determines the relative weighting of the data streams.
Could this be what's happening?
It might, but I have no way of answering this question without more information from you.
from pomp.
I'm closing this issue now, but feel free to re-open if more discussion is warranted.
from pomp.
Related Issues (20)
- looking for more information on `bsmc2()` HOT 1
- Error I'm getting for pfilter function HOT 2
- Error: in ‘bsmc2’: ‘dmeasure’ with log=TRUE returns illegal value HOT 4
- Global Search failing bc of one particle HOT 2
- Measurement model for multi-observed variables/states HOT 12
- How to simulate the time-varying parameter model by pomp package? HOT 14
- encounter problems when conducting trajectory fitting HOT 15
- MIF unable to fit to the data HOT 7
- Can we fit the time-series with the weighted likelihood? HOT 3
- Implementing a time delay in pomp HOT 1
- Using R's internal sample function in CSnippet HOT 7
- feel confused when choosing the best fitted model with loglik HOT 1
- An error associated with the installation and C snippets HOT 4
- Can I access the files related to Csnippet and directly edit? HOT 4
- I can't resolve this error when trying to run pomp HOT 4
- Error in adapting old syntax (traj.match) to the new one(traj_objfun) HOT 8
- Pseudo-likelihood combining iterated filtering and probe matching HOT 4
- How to estimate parameters of SIR model HOT 2
- How to add the exp(cubic spline) function to the deterministic model? HOT 14
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from pomp.