Measles Example

[stefanpetrevski]

Hello Dr. King,

I am trying to run the measles example (https://kingaa.github.io/clim-dis/parest/parest.html) and have gone through all the installations and FAQs to ensure that everything is in place.

I have replaced the "initializer" with "rinit" but then obtain an error: "Error: in ‘trajectory’: in ‘rinit’: variable 'gamma' not found among the parameters." It appears after the following block:

coef(niameyA,"Beta") <- beta.hat
x <- trajectory(niameyA,as.data.frame=TRUE)
dats <- join(as.data.frame(niameyA),x,by='time')

ggplot(dats,aes(x=time))+
  geom_line(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')

I am a bit perplexed as to what went wrong given all the code is the same as your examples.

[kingaa]

Have you furnished a value for the parameter "gamma"?

[stefanpetrevski]

Yes. In the definition of pomp, gamma is inputed as a parameter name. Then, when building the function which minimizes the square error for beta, gamma is set to 1. The complete code reads (taken from your website without the approximate R0 estimation at the beginning):

library(ggplot2) 
library(pomp) 
library(plyr)
library(reshape2)

niamey <- read.csv("http://kingaa.github.io/clim-dis/parest/niamey.csv")
ggplot(niamey,mapping=aes(x=biweek,y=measles,color=community))+
  geom_line()+geom_point()

pomp(
  data=subset(niamey,community=="A",select=-community),
  times="biweek",t0=0,
  skeleton=vectorfield(
    Csnippet("
      DS = -Beta*S*I/N;
      DI = Beta*S*I/N-gamma*I;
      DR = gamma*I;")),
  rinit=Csnippet("
      S = S_0;
      I = I_0;
      R = N-S_0-I_0;"),
  statenames=c("S","I","R"),
  paramnames=c("Beta","gamma","N","S_0","I_0")) -> niameyA

f1 <- function (beta) {
  params <- c(Beta=beta,gamma=1,N=50000,S_0=10000,I_0=10)
  sse(params)
}

beta <- seq(from=0,to=40,by=0.5)
SSE <- sapply(beta,f1)

plot(beta,SSE,type='l')
beta.hat <- beta[which.min(SSE)]
abline(v=beta.hat,lty=2)

coef(niameyA,"Beta") <- beta.hat
x <- trajectory(niameyA,as.data.frame=TRUE)
dats <- join(as.data.frame(niameyA),x,by='time')

ggplot(dats,aes(x=time))+
  geom_line(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')

Thank you very much.

[kingaa]

In f1, gamma is provided as a parameter. However, in niameyA, is there a parameter named "gamma"?

coef(niameyA)

[stefanpetrevski]

Yes it seems so:

paramnames=c("Beta","gamma","N","S_0","I_0")) -> niameyA

[kingaa]

No. coef(niameyA).

[stefanpetrevski]

Running

coef(niameyA)

gives an output of:

 Beta gamma 
    8     8 

So "gamma" is indeed the parameter. Following https://kingaa.github.io/clim-dis/parest/parest.html, if instead I specify:

coef(niameyA) <- c(Beta=beta.hat,gamma=1,N=50000,S_0=10000,I_0=10)
x <- trajectory(niameyA,as.data.frame=TRUE)
ggplot(data=join(as.data.frame(niameyA),x,by='time'),
       mapping=aes(x=time))+
  geom_line(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')

then the error is instead "Error in [.data.frame(x, by) : undefined columns selected".

[kingaa]

Good. What are the variable names in as.data.frame(niameyA)? Is there one called "time"?

[stefanpetrevski]

The output of the following is NULL:

coef(as.data.frame(niameyA))

That explains the error, but I must provide the arguments when I use data=join(..). What should I use instead, if no other variables appear?

[kingaa]

No. You can't apply coef to a data frame. Try: names(as.data.frame(niameyA)).

[stefanpetrevski]

I see, thank you for the correction.

names(as.data.frame(niameyA))

gives

[1] "biweek"  "measles"

But then if I change "time" to "biweek", so that:

coef(niameyA) <- c(Beta=beta.hat,gamma=1,N=50000,S_0=10000,I_0=10)
x <- trajectory(niameyA,as.data.frame=TRUE)
ggplot(data=join(as.data.frame(niameyA),x,by='biweek'),
       mapping=aes(x=time))+
  geom_line(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')

the error becomes "Error: All inputs to rbind.fill must be data.frames".

[kingaa]

Is x a data frame?

[stefanpetrevski]

No, typeof says it is a double. It shouldn't be, however, as it is storing a trajectory which we want to plot alongside the existing data for the Niamey outbreak. What should be corrected?

[kingaa]

With version 2, trajectory syntax changed. To get a data frame result, use format="data.frame" instead of as.data.frame=TRUE.

[kingaa]

There is a document detailing all the changes that came with version 2.

[stefanpetrevski]

I see. I spotted the "initializer" to "rinit" change in the document, but I didn't see the data formatting change, thank you very much. So now I have changed the definition of x to:

x <- trajectory(niameyA, format="data.frame")

and it follows that x is now a list. To reflect this change, I should also modify the new syntax for:

ggplot(data=join(as.data.frame(niameyA),x,by='biweek'),
       mapping=aes(x=time))+
  geom_line(aes(y=measles),color='black')+
  geom_line(aes(y=I),color='red')

So finally, how to join niameyA as a data.frame with this new format command?

[kingaa]

x is a list, but is also a data frame (all data frames are lists). Does the ggplot code above not work?

[stefanpetrevski]

Yes, I have just modified the aesthetics from x=time to x=biweek. The new reformatting was sufficient, ggplot does not need to be modified and the code now works. Thank you very much, your comments were extremely helpful!

[kingaa]

Glad to hear it.