kingaa
commented
11 months ago First of all, I point you to FAQ 1.4. In particular, it appears that the code you are running is from a very old version of pomp and much has changed. Many of these were in the transition to version 2. These are described in the pomp version 2 upgrade guide linked there.
However, not all of the problems were there. Some are due to simple typos (for example you refer to a parameter with the invalid name 21 in your call to parameter_trans. (Incidentally, this is the source of the specific error you mentioned). There are other mistakes as well. Rather than explain these all, I simply ask you to compare your code with the following.
library(tidyverse)
library(pomp)
stopifnot(packageVersion("pomp")>="5.4")
theme_set(theme_bw())
set.seed(1173439184)Download the data from here: LV_data.csv
read_csv("LV_data.csv") |>
rename(time=Time,species=Species) |>
group_by(species,time) |>
mutate(rep=seq_along(time)) |>
ungroup() |>
arrange(time,rep,species) -> dat
dat |>
ggplot(aes(x=time,y=counts,color=species,
fill=species,group=species))+
geom_point(size=2)+
geom_smooth(se=FALSE,linetype=2)+
labs(y="counts")dat |>
unite(col="name",species,rep) |>
pivot_wider(values_from="counts") |>
arrange(time) |>
traj_objfun(
times="time",
t0=0,
skeleton=vectorfield(
Csnippet("
DN1 = r1*N1*(K1-N1-a12*N2)/K1;
DN2 = r2*N2*(K2-N2-a21*N1)/K2;"
)
),
rinit=Csnippet("
N1 = N1_0;
N2 = N2_0;"
),
dmeasure=Csnippet("
lik = dnbinom_mu(N1_1,1/theta,p*N1,1)+
dnbinom_mu(N1_2,1/theta,p*N1,1)+
dnbinom_mu(N1_3,1/theta,p*N1,1)+
dnbinom_mu(N2_1,1/theta,p*N2,1)+
dnbinom_mu(N2_2,1/theta,p*N2,1)+
dnbinom_mu(N2_3,1/theta,p*N2,1);
lik = (give_log) ? lik : exp(lik);"
),
rmeasure=Csnippet("
N1_1 = rnbinom_mu(1/theta,p*N1);
N1_2 = rnbinom_mu(1/theta,p*N1);
N1_3 = rnbinom_mu(1/theta,p*N1);
N2_1 = rnbinom_mu(1/theta,p*N2);
N2_2 = rnbinom_mu(1/theta,p*N2);
N2_3 = rnbinom_mu(1/theta,p*N2);"
),
partrans=parameter_trans(
log=c("theta","r1","r2", "a12", "a21", "K1", "K2"),
logit=c("p")
),
params=c(
N1_0=10,N2_0=10,
r1=0.5,r2=0.5,
K1=180,K2=100,
a12=0.3,a21=0.4,
theta=0.9,p=0.8
),
est=c("r1","r2","a12","a21","theta","p","K1","K2"),
statenames=c("N1","N2"),
paramnames=c("p","theta","r1","r2","a12","a21","K1","K2","N1_0","N2_0")
) -> fn
library(subplex)
subplex(
fn=fn,
par=coef(fn,c("r1","r2","a12","a21","theta","p","K1","K2"),transform=TRUE),
control=list(maxit=50000)
) -> fit
fit
fn(fit$par)
coef(fn)
fn |>
trajectory(
times=seq(0,200,by=1),
format="data.frame"
) |>
pivot_longer(
c(N1,N2),
names_to="species",values_to="counts"
) |>
mutate(
counts=coef(fn,"p")*counts
) -> fitted
dat |>
ggplot(aes(x=time,y=counts,color=species,
fill=species,group=species))+
geom_point(size=2)+
geom_smooth(se=FALSE,linetype=2)+
geom_line(data=fitted)+
labs(y="counts")Also, note that no optimization is performed because you did not ask for any to be done. traj_objfun merely constructs an objective function. You must run an optimizer of some sort to actually find its minima. See the call to subplex in the above.
zz100chan
Hello, I am trying to test out one of the example shown earlier here by @kingaa from the post (here) But since the example was from an older version, the traj.match does not work anymore and I am trying to adapt it to the newer syntax of traj_objfun with the following code:
But errors come up firstly from the partrans function with the following error about a parameter of "T_21" which never appears. Removing the partrans, the code did run but then fit gives the same parameter as params, i.e. no estimation is done. I am wondering if you have insights in this issue. Thank you in advance!