# this script uses a modification of the function "ssm" from the package "sspir"
# to compare differences in estimates of states in models with similar form, but
# different calls with respect to the glm-esque interface

library(sspir)
library(boot)
source("ssm2")

# generate some data
# 2 cols are "successes" and "failures"
y <- ts(cbind(rpois(20,20), rpois(20,100)), start=1, freq=1)

#------------------------------
# use ssm with binomial errors
#------------------------------

# set the ts model via ssm
ssm.ex1 <- ssm2(y ~ tvar(1), family=binomial(link="logit"), m0=c(y[1,1]/(y[1,1]+y[1,2])), fit=TRUE)
# get the model fit
ssm.ex1.fit <- getFit(ssm.ex1)

# this bit uses another approach, but returns same results
# set the ts model via ssm
#ssm.ex1 <- ssm2(y ~ tvar(1), family=binomial(link="logit"), m0=c(y[1,1]/(y[1,1]+y[1,2])), fit=FALSE)
# get the model fit
#ssm.ex1.fit <- extended(ssm.ex1$ss)
#ssm.ex1.fit <- getFit(ssm.ex1.fit)

# back-transform the estimates
ssm.ex1.pre <- inv.logit(ssm.ex1.fit$m)
ssm.ex1.up <- inv.logit(ssm.ex1.fit$m-unlist(ssm.ex1.fit$C))
ssm.ex1.lo <- inv.logit(ssm.ex1.fit$m+unlist(ssm.ex1.fit$C))


#------------------------------
# use ssm with gaussian errors
#------------------------------

# ts of the logit(p)
y2 <- ts(logit(y[,1]/(y[,1]+y[,2])), start=1, freq=1)

# set the ts model via ssm
ssm.ex2 <- ssm2(y2 ~ tvar(1), family=gaussian, m0=y2[1], fit=TRUE)
# get the model fit
ssm.ex2.fit <- getFit(ssm.ex2)

# back-transform the estimates
ssm.ex2.pre <- inv.logit(ssm.ex2.fit$m)
ssm.ex2.up <- inv.logit(ssm.ex2.fit$m-unlist(ssm.ex2.fit$C))
ssm.ex2.lo <- inv.logit(ssm.ex2.fit$m+unlist(ssm.ex2.fit$C))


#-----------------
# draw some plots
#-----------------

# plot the ts of logit(p) for the data
dev.new()
plot(y2, type="b", pch=16, ylab="logit(p)")
# plot the ts of estimated logit(p)
points(ssm.ex1.fit$m, type="b", pch=16, col="red")
points(ssm.ex2.fit$m, type="b", pch=16, col="blue")

# plot the ts of p for data (black) & model (blue)
dev.new()
par(mfrow = c(2, 1))
plot(y[,1]/(y[,1]+y[,2]), type="b", pch=16, ylim=c(0.1, 0.3), ylab="p")
# plot the estimated RW
points(ssm.ex1.pre, type="b", pch=16, col="blue")
points(ssm.ex1.up, type="l", pch=16, col="red")
points(ssm.ex1.lo, type="l", pch=16, col="red")

plot(y[,1]/(y[,1]+y[,2]), type="b", pch=16, ylim=c(0.1, 0.3), ylab="p")
# plot the estimated RW
points(ssm.ex2.pre, type="b", pch=16, col="blue")
points(ssm.ex2.up, type="l", pch=16, col="red")
points(ssm.ex2.lo, type="l", pch=16, col="red")
