Greetings,
I sent the message below to the developer of the contributed R package
"sspir", but have yet to receive any response. I would be very
grateful
for any advice people have on the matter.
Thanks,
Mark
-------- Original Message --------
Subject: possible bug in sspir?
Date: Tue, 19 May 2009 16:08:41 -0700
From: Mark Scheuerell <mark.scheuerell at noaa.gov>
To: aas.claus.dethlefsen at nja.dk, dethlef at math.aau.dk, cld at rn.dk
Hi Claus,
I have been using the "sspir" package that you developed for use with
R
and have found it very useful. Previously, all of the time series
models I had been fitting had gaussian errors, and everything seemed to
work quite well.
Recently, however, I became interested in some models for time series of
survival probabilities. Given the response variable is a proportion that
lies on [0,1], it seemed reasonable to fit a model with a binomial error
structure and logit link function. It is my understanding that if I were
to fit a similar model with "glm", the response variable can take one
of
3 forms: (1) a vector with 2 levels (eg, 0 and 1) if the data are from
individuals; (2) a vector of the proportions in which case the
additional function call "weights" would be needed with a vector of
the
number of trials; or (3) a 2xN matrix where the columns are successes
and failures, respectively.
When I try to fit a basic random-walk model using the "ssm" function,
however, I can only get option (3) above to work. That in itself is OK,
but it seems as though ssm and subsequent functions (eg, "kfs") want
to
work on a Nx2 matrix instead. For example, the following lines should
produce an estimated state vector of length=20, but the result is length=2:
>y <- ts(cbind(rpois(20,20), rpois(20,100)), start=1, freq=1)
>ssm.fit <- ssm(t(y) ~ tvar(1), family=binomial(link="logit"),
fit=FALSE)
>ssm.ex1.fit <- extended(ssm.ex1$ss)
>ssm.ex1.fit$m
Time Series:
Start = 1
End = 2
Frequency = 1
Series 1
[1,] -0.0498745
[2,] -0.1357787
If I try to trick the code and pass in the transpose of y, I get this error:
>y <- t(y)
>ssm.fit <- ssm(y ~ tvar(1), family=binomial(link="logit"),
fit=FALSE)
>ssm.ex1.fit <- extended(ssm.ex1$ss)
Error in ss$Fmat(tt, ss$x, ss$phi) : subscript out of bounds
When I went in and looked at the details of the ssm code, I saw that on
lines 146-151, the code is trying to establish the vectors for the
number of trials ("ntotal") and the response ("y") based on
the rows of
the input matrix y rather than based on the columns. Thus, if I change
those lines and rerun the code above, I do in fact get an estimated
state vector of length=20 (and the estimates seem reasonable). My
modified version of your code is attached to this msg (ssm2.r).
I also have an additional question that may or may not be related to
this. If I fit a random-walk model using the approach outlined above and
my ssm2 function (assuming it is correct), the point estimates of the
state vector are nearly identical to the "observed" data. However, if
I
first do a logit transform of the proportions and then use ssm (or my
version) with gaussian errors, the point estimates of the state vector
are not nearly as close as those obtained with the binomial error
structure (see my attached code "ssm_bin_ex.r"). I can understand why
the variance of the estimates would differ, but why the means? If I fit
some dummy models for non-time series models with glm, the estimates of
the means using the 2 approaches just outlined are nearly identical (see
my attached code "glm_bin_ex.r").
I would be sincerely grateful for any help or insights you could offer
with respect to my problems.
Best wishes,
Mark
_____________________________
Mark D. Scheuerell, Ph.D.
National Marine Fisheries Service
Northwest Fisheries Science Center
2725 Montlake Blvd E
Seattle, WA 98112
USA
206.302.2437 tel
206.860.3400 fax
http://www.nwfsc.noaa.gov/mark.scheuerell
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: ssm_bin_ex.r
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20090529/13a248ed/attachment-0006.pl>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: glm_bin_ex.r
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20090529/13a248ed/attachment-0007.pl>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: ssm2.r
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20090529/13a248ed/attachment-0008.pl>
