[R] replicated time series - lme?
Spencer Graves
spencer.graves at pdf.com
Tue Mar 7 03:56:35 CET 2006
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
>
More information about the R-help
mailing list