[R] replicated time series - lme?

[Spencer Graves]

	  What are your primary research objectives?  Are you primarily 
interested in a simple model of how the system evolves over time?  Or do 
your 17 different series include variations in simulation conditions?

	  If I were largely interested in understanding better how the system 
evolves over time, I might approach the analysis as follows:

	  1.  Make separate normal probability plots of all series.  If they 
aren't normal, can I transform to normality (e.g., taking logarithms -- 
or using the famous Box-Cox transformation, making lambda plots as 
described in VR = Venables and Ripley (2004) Modern Applied Statistics 
with S, 4th ed. (Springer) and the companion MASS package for R).

	  2.  If I could get each of the 17 normal plots to look relatively 
normal (possibly after a transformation), I'd then try to fit a separate 
ARIMA model to each of the 17 series using acf, pacf, arima, etc., as 
described in VR (see step # 1).

	  NOTE:  You do NOT want to ask lme to estimate an autoregression of 
order 100.  That's nonsense.  The genius of the Box-Jenkins / ARIMA 
system is that it provides a very flexible system for modeling series 
with potentially very long lags with parsimoneous models expressed as 
rational functions (i.e., ratios of polynomials) in the back-shift 
operator.  To understand this, see any good reference on ARIMA modeling, 
e.g., the original Box and Jenkins, Time Series Analysis, Forecasting 
and Control.

	  3.  After you have fit parsimoneous models to the several of the 17 
series, you should begin to see a pattern in what kinds of models should 
fit reasonably well.  Then you can return to "lme" to look for further 
patterns in how the parameters of the parsimoneous time series models 
are related.

	  Does this help?
	  spencer graves

Katrin Meyer wrote:

> Dear R-helpers,
> 
> I have a time series analysis problem in R:
> 
> I want to analyse the output of my simulation model which is proportional
> cover of shrubs in a savanna plot for each of 500 successive years. I have
> run the model (which includes stochasticity, especially in the initial
> conditions) 17 times generating 17 time series of shrub cover. 
> 
> I am interested in a possible periodicity of shrub cover (the time lag(s) of
> significant autocorrelation) and would like to include the information from
> the 17 replicate simulations. 
> 
>>From the acf and pacf plots, I can assume that time lags of up to 150 years
> are relevant. I tried  a mixed effects model as following:
> 
> #short example version of input file with 2 runs and 5 time steps (instead
> of 17 runs and 500 time steps)
> run	t	cover
> 1	1	0.234306
> 1	2	0.188896
> 1	3	0.198193
> 1	4	0.213959
> 1	5	0.184952
> 2	1	0.189316
> 2	2	0.185631
> 2	3	0.20211
> 2	4	0.216064
> 2	5	0.216064
> 
> #calculate the correlation of lag 1 over 17 replicates
> a<-0
> for (i in 1:17)
> {
> c<-ts( cover[run==i] )
> d<-acf( c, lag=1, plot=F)
> a<-a+d$acf[2]
> }
> a<-a/17
> a
> #[1] 0.9021463
> 
> #mixed effects model
> model1<-lme(cover~t,random=~t|run, method="ML")
> model2<-update(model1,correlation=corCAR1(0.902,form=~t|run))
> anova(model1,model2)
> 
> But this just gives significance for a lag of 1, so I tried to find out the
> correlation at greater lags with arima to be able to use corARMA() as
> correlation structure:
> 
> arima(cover[run==1],order=c(100,0,0)) 
> #does not work: “error in polyroot(z): polynomial degree too high”
> 
> Any ideas to solve this? Maybe, I don’t even need a mixed effects model?
> 
> I would be very grateful for any help
> Katrin
>