[R] possible bug in "sspir" package?

[Mark Scheuerell]

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>