R Manual for Statistical Models

M.C.M. de Gunst and M.P. McAssey

Department of Mathematics, VU University Amsterdam

12/3/2013
Contents

1 Introduction 3

2 Analysis of Variance in R 5

3 Nonlinear regression in R 10

4 Generalized Linear Models in R 14

5 Time series in R 18
5.1 Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
5.2 ARIMA models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

6 Some literature on R 24

2
Chapter 1

Introduction

These notes provide a summary of the main R functions, data structures, and relevant ex-
amples, for each chapter of the lecture notes “Statistical Models” by M.C.M. de Gunst. For
more extended documentation see the R-project website [Link]
and the literature provided at the end of these notes.
In the sequel the symbol ‘>’ denotes the R command prompt. Each R command must
be entered in the R console after the command prompt. The prompt itself is not typed
by the user.
For the analysis of data based on statistical models using R, three new R-concepts
need to be introduced. The first is the object type factor. A factor is a vector which
is encoded with an additional feature designated as “Levels.” A factor is understood to
be a categorical vector, in which each component belongs to one of several categories.
For example, if an experiment involves measurements on 12 subjects who are randomized
among four treatment groups A, B, C or D, the data could be stored in a factor named
trtmnt as follows:

> trtmnt <- factor(c("D","A","C","B","A","C","D","B","B","C","D","A"))

To view the contents of trtmnt, we type, as usual, its name at the prompt:

> trtmnt
[1] D A C B A C D B B C D A
Levels: A B C D

The next new concept is that of a dataframe. The implementation of many R functions
is more convenient when, or requires that, the data are provided as a ‘dataframe’ object.
Also, several R functions have a dataframe as output. A dataframe is an R data structure
which is a rectangular array like a matrix, but each column need not be of the same data
type. For instance, one column can be a numeric vector, another column can be a character
vector, and a third column can be a factor, as long as each column has the same length.
The R function [Link] may be used to create a dataframe. For instance, suppose
we have four R vectors: y containing the observations, and x1, x2, and x3 containing the
values of three explanatory variables, created by

> y <- c(8,5,9)

> x1 <- c(2,3,4)

3
4 Chapter 1. Introduction

> x2 <- c("a","b","c")

> x3 <- factor(c(1,2,3))

Then a dataframe mydata can be created by the command

> mydata <- [Link](y, x1, x2, x3)

Its content looks like:

> mydata
y x1 x2 x3
1 8 2 a 1
2 5 3 b 2
3 9 4 c 3

Note that rows and columns are assigned names automatically. The components can be
extracted from a dataframe in a similar way as from a matrix, and can be used in the
same way as the elements of the vectors from which they originate:

> mydata[1,2] # extract component in first row, second column

[1] 2
> mydata[,1] # extract all components of first column
[1] 8 5 9
> mydata$y # extract column named "y"
[1] 8 5 9
> mydata[[2]] # extract second column
[1] 2 3 4

A matrix can be converted to a dataframe using the R command [Link]. Also,

if the command [Link] is used to read data into R from an external text file, the
resulting object will be a dataframe. Note that in that case the columns that contain
characters are objects of type factor, whereas the columns that contain levels represented
by integers are interpreted to be of type integer. To keep given column names or to assign
row and/or column names while using [Link], the arguments header, [Link] and
[Link], respectively, can be used. This is further explained on the R help page for
[Link].
The third new concept is a formula. For many functions the formula is a required
argument. It specifies the model, that is, the relationship between the observations and
the explanatory variables that is to be used in the analysis. The argument formula is
indeed a formula like a mathematical formula. Specific examples of formulas and their use
will be provided in the sequel.
Chapter 2

Analysis of Variance in R

The basic R expression for fitting an analysis of variance (ANOVA) model is

aov(formula, data)
where aov is short for ‘analysis of variance’, formula gives the relationship between the
response variable and the effects and their interactions, and data is a dataframe containing
the data to be analyzed. If the data argument is missing, R will look for the variables
used in the formula argument within the user’s local environment (that is, among the
variables currently existing in the R workspace). In general, it is much more convenient
to maintain all data pertaining to a set of subjects together in one data structure, the
dataframe. Hence we will consistently use a dataframe to pass data into R functions.
Consider the blood pressure example described in the lecture notes, involving 36 sub-
jects. In this example we investigate in what way the systolic blood pressure depends on
age and socio-economic class. The dataframe named bpdata contains the data pertaining
to the example. The first six and last six lines of bpdata are
bloodpressure age socio
1 116 1 1
2 108 1 1
3 160 1 1
4 116 1 1
5 136 1 2
6 124 1 2
. . . .
. . . .
. . . .
31 188 3 2
32 158 3 2
33 214 3 3
34 146 3 3
35 138 3 3
36 148 3 3
The first column above denotes the observation number, which is automatically supplied
by R. The second column provides the observed systolic blood pressure for each subject.

5
6 Chapter 2. Analysis of Variance in R

The numbers 1, 2 and 3 in the column ‘age’ refer to the age levels of the subjects, while
the numbers in the column ‘socio’ refer to the three socio-economic levels. Note that the
numbers in the last two columns represent levels of the corresponding categorical variables.
The choice of the representation is arbitrary; the numbers could be replaced by the letters
“A”, “B” and “C” , for instance. Hence these two columns should be of the R class
“factor”, and if not they must be converted into factors.
How to create such a dataframe in R? As mentioned in the introduction, if the data are
stored in an external text file, the R function [Link] can be used to create bpdata.
In that case bpdata is a dataframe by default. However, the two columns containing
the effects consist of integers, and hence they are interpreted by [Link] as integer
vectors that are not of the factor format. We will discuss below how to convert the effects
into factors within the formula argument of the aov function so that it computes correct
results. However, one can also do the conversion in advance. For example,

> bpdata$age <- factor(bpdata$age)

> bpdata$socio <- factor(bpdata$socio)

will transform the two effects in bpdata into the required factor type.
But what if the data are not available in an external text file? As illustrated in the
introduction, we may first create the different components of the data frame and then use
the R function dataframe to create the data frame. Since systolic blood pressure is a
quantitative variable, an R vector bloodpressure containing the measured systolic blood
pressures can be constructed in the usual manner. The factor age can be made with the
command

> age <- factor(c(rep(1,12), rep(2,12), rep(3,12)))

where rep is the replication function. Since ‘age’ is categorical rather than quantitative,
we must use the factor function to transform the numerical vector into a factor. The
factor socio can be created in a similar way. Finally, to create the dataframe bpdata we
use the command

> bpdata <- [Link](bloodpressure, age, socio)

To see just the first few rows of a dataframe, the head function is useful. For this example,
type head(bpdata) at the R prompt.
Once the dataframe containing the data is available in the R workspace, we are ready
to analyze the relationship between the response variable and the the effects with the R
function aov. The first argument of aov is a formula which establishes that relationship.
In our example the formula can be given by

formula = bloodpressure ~ 1 + age + socio

where 1 indicates that an intercept (or general mean) µ is present in the model. Since we
usually want an intercept, R includes it by default. Thus this formula is equivalent to

formula = bloodpressure ~ age + socio

 7

Excluding the intercept can be done by using -1 in the formula

formula = bloodpressure ~ age + socio - 1

If the effects in the formula are not already in the factor format within the dataframe,
they must be transformed within the formula, such as

formula = bloodpressure ~ -1 + factor(age) + factor(socio)

We will assume below that the effects are already in the factor format. If you are not
sure, it will never hurt to use the factor function on an effect that is already in the factor
format.
To perform an ANOVA on the dataset bpdata and store the output in an object named
bpfit we could use the command

> bpfit <- aov(formula = bloodpressure ~ age + socio, data = bpdata)

or the shorter version

> bpfit <- aov(bloodpressure ~ age + socio, bpdata)

If the model is to contain all the factors in the dataframe, this command can be shortened
further by

> bpfit <- aov(bloodpressure ~ . , bpdata)

The output of the R function aov is an object of class aov. This means that several
generic functions, like summary or plot, can extract specific “aov-related” information
from it. There are also R functions that require aov objects. If a generic and a specific
“aov” function do the same thing, often the latter one does it more efficiently or gives
more information. To illustrate this, the resulting ANOVA table (Table 1.3 in the lecture
notes) can be found by using the summary function as well as the function anova. The
latter would give

> anova(bpfit)
Analysis of Variance Table

Response: bloodpressure
Df Sum Sq Mean Sq F Value Pr(>F)
age 2 6890.4 3445.2 6.5438 0.004258 **
socio 2 747.7 373.9 0.7101 0.499409
Residuals 31 16320.9 526.5
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

In the column ‘F Value’ the values of the test statistics FA and FB , as defined in the
lecture notes, are given. The column ‘Pr(>F)’ gives the corresponding p-values for the
significance of each factor. The symbol ** indicates that this p-value is between 0.001 and
8 Chapter 2. Analysis of Variance in R

0.01, according to the “Signif. codes” legend below the ANOVA table. If instead of anova
the function summary would have been used, the numbers would have been the same, but
the two first lines would have been omitted. Check this yourself.
Until now we have only considered a model without interactions. Suppose we want to
include an interaction term between ‘age’ and ‘socio-economic class’. The command to fit
the model with this interaction term is
bpfit <- aov(bloodpressure ~ age + socio + age:socio, bpdata)
or, more efficiently,
> bpfit <- aov(bloodpressure ~ age*socio, bpdata)
Here age:socio stands for the interaction between ‘age’ and ‘socio-economic class’, while
age*socio represents factor crossing, i.e., each factor individually and their interaction.
If the model is to contain all the factors in the dataframe, and all interaction terms are to
be included in the model, the command can be expressed in the shorter form
> bpfit <- aov(bloodpressure ~ .^2 , bpdata)
The ANOVA table is obtained in the same way, using either the summary or anova com-
mands:
> anova(bpfit)
Analysis of Variance Table

Response: bloodpressure
Df Sum Sq Mean Sq F Value Pr(>F)
age 2 6890.4 3445.2 6.3579 0.005463 **
socio 2 747.7 373.9 0.6899 0.510236
age:socio 4 1690.1 422.5 0.7797 0.548171
Residuals 27 14630.8 541.9
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1
In general, if we wish to include all kth-order interaction terms, and the model is to contain
all the factors in the dataframe, we use the command
> bpfit <- aov(Y ~ .^k , mydata)
for observations Y and a data frame mydata; Y must be the first column of the data frame
mydata.
Note that the functions anova and summary do not provide an estimate for the coeffi-
cient parameter vector β, which consists of the general mean and the model effects. We
can compute the estimate easily in R, since the expressions for each coefficient in terms
of the observations contain only the sample means (see (1.19) and (1.24) in the lecture
notes). However, with the R function [Link] we can directly extract this infor-
mation from the output of aov. For the blood pressure example, to obtain the estimated
effects of the levels for both factors and their interactions, use
 9

> [Link](bpfit)

For estimates of the standard errors of the effects, add the optional setting se=T. To
obtain an estimate for the general mean µ, include the argument type="means". This will
produce a table of sample means; the number entitled Grand mean is the estimate µ̂ as
presented in the lecture notes. In this case, if the setting se=T is also present, the output
will include the standard errors for the differences between the sample means.
Chapter 3

Nonlinear regression in R

The process of fitting nonlinear regression models in R is similar to that for fitting linear
regression models. However, there are two important differences:

– expressions for the parameter estimates cannot be derived analytically, so an iterative

numerical procedure is used, for which initial values will be needed.

– the regression function is more general than in the linear case, so a more flexible
protocol will be needed.

The R function for fitting a nonlinear regression model is nls. This function requires the
argument

– formula: a formula which specifies the nonlinear regression model.

It also has several optional arguments; the most important ones are

– data: a dataframe containing the observations on the response variable and the
values of explanatory variables. If this argument is missing, R will look for the
variables used in the formula argument within the R workspace.

– start: a list or numerical vector specifying initial parameter values. Although it is

optional, use of start is recommended for unambiguous specification of the param-
eters. If start is omitted, the assumption is that any names occurring in formula
that are not variables in the dataframe are parameters that must be estimated.

The output of the function nls is an object of class nls containing the estimated parameter
values, among others. See the R help facility for nls for further details.
Let us illustrate the use of nls by means of the Puromycin example. The dataframe
Puromycin is a built-in R dataframe:

> Puromycin
conc rate state
1 0.02 76 treated
2 0.02 47 treated
. . . .
. . . .

10
 11

. . . .
12 1.10 200 treated
13 0.02 67 untreated
14 0.02 51 untreated
. . . .
. . . .
. . . .
23 1.10 160 untreated
For now we extract the data only for the 12 treated cases:
> pur.t <- Puromycin[Puromycin$state == "treated",]
The result pur.t is also a data frame, consisting of the data for the treated cases only:
> pur.t
conc rate state
1 0.02 76 treated
. . . .
. . . .
. . . .
12 1.10 200 treated
For the analysis we give the command
> [Link] <- nls(rate ~ theta1*conc/(theta2+conc), pur.t,
list(theta1=195.8027, theta2=0.04840654))
which computes optimal parameter values and stores the results in the object [Link].
Here rate ~ theta1*conc/(theta2+conc) is the argument formula, pur.t is the ar-
gument data, and list(theta1=195.8027, theta2=0.0480654) is the argument start
giving the initial parameter values for the iterative parameter estimation procedure. The
argument formula always consists of the vector of observations (rate in this case) con-
nected to the regression function (here theta1*conc/(theta2+conc)) by means of the
symbol ~ . Obviously, the regression function is expressed as the mathematical formula
which defines it. The results of the analysis, as reported in the lecture notes, can be ob-
tained by typing the name of the nls object (in this example, [Link]) at the command
prompt:
> [Link]
Nonlinear regression model
model: rate ~ theta1 * conc/(theta2 + conc)
data: pur.t
theta1 theta2
212.68356 0.06412
residual sum-of-squares: 1195

Number of iterations to convergence: 5

Achieved convergence tolerance: 9.852e-06
12 Chapter 3. Nonlinear regression in R

As in the analysis of variance, more information—but now “nls-related”—is given by the

function summary. We can apply this function to the object in which the output of nls is
stored:

> summary([Link])

Formula: rate ~ theta1 * conc/(theta2 + conc)

Parameters:
Estimate Std. Error t value Pr(>|t|)
theta1 2.127e+02 6.947e+00 30.615 3.24e-11 ***
theta2 6.412e-02 8.281e-03 7.743 1.57e-05 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 10.93 on 10 degrees of freedom

Number of iterations to convergence: 5

Achieved convergence tolerance: 9.852e-06

The t value is the t-ratio defined by the ratio of the point estimate (Estimate) to the
estimated standard deviation of the point estimate (Std. Error). The other results are
self-evident.
To obtain the residuals and store them in a vector names resid, use the command

> resid <- summary([Link])$residuals

The more complicated analysis of the full model which involves both the treated and
untreated groups, requires an indicator variable (called x2 in the lecture notes) as an
additional explanatory variable. It can easily be obtained by transforming the column
Puromycin$state into the appropriate 0–1-variable, for instance by means of the com-
mands

> pur <- Puromycin

> pur$state <- ifelse(Puromycin$state == "treated", 1, 0)

The column pur$state in pur now represents the desired ”x2-variable”. Applying nls
yields

> [Link] <- nls(rate ~ ((theta1+fi1*state)*conc)/((theta2+fi2*state)+conc),

pur, list(theta1=160, theta2=0.04, fi1=52, fi2=0))
> [Link]
Nonlinear regression model
model: rate ~ ((theta1 + fi1*state)*conc)/((theta2 + fi2*state) + conc)
data: pur
theta1 theta2 fi1 fi2
 13

160.27996 0.04771 52.40371 0.01641

residual sum-of-squares: 2055

Number of iterations to convergence: 6

Achieved convergence tolerance: 4.303e-06

and

> summary([Link])

Formula: rate ~ ((theta1 + fi1*state)*conc)/((theta2 + fi2*state) + conc)

Parameters:
Estimate Std. Error t value Pr(>|t|)
theta1 1.603e+02 6.896e+00 23.242 2.04e-15 ***
theta2 4.771e-02 8.281e-03 5.761 1.50e-05 ***
fi1 5.240e+01 9.551e+00 5.487 2.71e-05 ***
fi2 1.641e-02 1.143e-02 1.436 0.167
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Residual standard error: 10.4 on 19 degrees of freedom

Number of iterations to convergence: 6

Achieved convergence tolerance: 4.303e-06

From the above we see that an analysis with more than one explanatory variable is
performed in exactly the same way as an analysis with only one explanatory variable.
Moreover, R provides us with everything in which we are currently interested: all other
results mentioned in the lecture notes can be obtained based on the output of the function
nls.
Chapter 4

Generalized Linear Models in R

Fitting generalized linear models is in many ways similar to fitting linear models. This
is evident from the types of R functions that are available for the analysis of data under
GLMs. For fitting generalized linear models in R the most important R function is glm,
a function similar to the function lm for fitting linear regression models to data stored
in dataframes. The function glm has the required argument formula, and, among other
ones, the optional argument family and the optional argument data:

– formula: a formula expression with the response on the left of the ~ operator, and
the explanatory variables, separated by + operators, on the right.

– family: an expression that defines the distribution and the link function. The
default is the normal distribution (gaussian), with the identity link. The pos-
sible distributions are gaussian, binomial, poisson, Gamma, [Link],
and quasi. Family distributions can also take arguments, as in binomial(link =
probit). Each family has a default link, so it is not necessary to specify the link
function unless a link different from the default is preferred. For details, see the help
page for family.

– data: a dataframe containing the observations on the response variable and the
values of the explanatory variables.

The output of glm is an object of class glm.

The format of the data must be correct before using the glm function. For logistic
regression, the response is usually a binary variable, e.g., a 1 for the presence of a condition
or a 0 for its absence. In this case, the response should be of the factor class. However,
the response can also be given in two columns, as when each row of data pertains not to
a single subject but to a group of subjects. In this case, the first column is the number of
subjects in each group who meet a certain condition, and the second column is the number
of subjects in each group who do not meet the condition. The response should then be of
the matrix class, with two columns of numeric elements.
Likewise, some or all of the predictors can be factors, and should be set as such. For
example, a predictor might be the sex of the subject, or the ethnicity of the subject. In this
situation, the glm function will automatically choose one level of the factor as the baseline
and each additional level as a binary predictor, so that the model will be identifiable.

14
 15

For the analysis of the Kyphosis data from the example in the lecture notes we chose
the binomial family, corresponding to logistic regression. Its default link function is the
logit function, so we do not need to specify this. The response is in factor format, and
the three explanatory variables are numeric. For the analysis under a logistic regression
model with the three explanatory variables, we give the command

> [Link] <- glm(Kyphosis ~ Age + Number + Start, family = binomial,

data=kyphosis)

obtaining the results

> [Link]
Call:
glm(formula = Kyphosis ~ Age + Number + Start, family = binomial,
data = kyphosis)

Coefficients:
(Intercept) Age Number Start
-2.03693 0.01093 0.41060 -0.20651

Degrees of Freedom: 80 Total (i.e. Null); 77 Residual

Null Deviance: 83.23
Residual Deviance: 61.38 AIC: 69.38

We see that an intercept is included in the model by default. Deleting it from the model
can be done by adding the “explanatory variable” -1

> glm(Kyphosis ~ -1 + Age + Number + Start, family = binomial,

data = kyphosis)

Interactions between variables can be included in the same way as in the analysis of
variance. For instance,

> [Link] <- glm(Kyphosis ~ Age + Number + Start + Age:Number,

family = binomial, data = kyphosis)

or, including all main effects and all pairwise interactions,

> [Link] <- glm(Kyphosis ~ . ^2, family = binomial, data = kyphosis)

As before, we can get more details using the summary function, but now this information
is “glm-related”:

> summary([Link])

Call:
glm(formula = Kyphosis ~ Age + Number + Start, family = binomial,
data = kyphosis)
16 Chapter 4. Generalized Linear Models in R

Deviance Residuals:
Min 1Q Median 3Q Max
-2.3124 -0.5484 -0.3632 -0.1659 2.1613

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.036934 1.449575 -1.405 0.15996
Age 0.010930 0.006446 1.696 0.08996 .
Number 0.410601 0.224861 1.826 0.06785 .
Start -0.206510 0.067699 -3.050 0.00229 **
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 83.234 on 80 degrees of freedom

Residual deviance: 61.380 on 77 degrees of freedom
AIC: 69.38

Number of Fisher Scoring iterations: 5

Inclusion of the argument correlation=T outputs in addition the estimated correlations
between the estimators for the coefficients. We note that what is called the Residual
Deviance here is equal to what is called the deviance D in the lecture notes; it is the de-
viance D corresponding to the model that was fitted. The Null Deviance is the deviance
D corresponding to the model with no predictors (but with an intercept if the fitted model
has an intercept). To obtain the deviance residuals, use the residuals function with the
name of the fitted model. For example, residuals([Link]). The sum of the squares
of these residuals should equal D.
To obtain an Analysis of Deviance table, use the anova function. Although this func-
tion usually pertains to analysis of variance, it will recognize a glm object and provide
an analysis of deviance instead. To include a test for the significance of each predictor,
add the optional setting test="Chisq". One can use test="F" or test="Cp" instead, but
these tests are not appropriate for logistic or Poisson regression, and you will receive a
message stating so. For example:
> anova([Link], test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: Kyphosis

Terms added sequentially (first to last)

 17

Df Deviance Resid. Df Resid. Dev P(>|Chi|)

NULL 80 83.234
Age 1 1.302 79 81.932 0.253851
Number 1 10.306 78 71.627 0.001326 **
Start 1 10.247 77 61.380 0.001369 **
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

The Pearson residuals may be obtained using the residuals function with the setting
type="pearson". To obtain the Pearson chi-squared statistic P , compute the sum of the
squares of the Pearson residuals. For example:

> P <- sum( residuals([Link], type="pearson")^2 )

produces a value of P = 70.3152.

By applying the R function plot to the output of glm, e.g., plot([Link]), several
plots can be produced for the fitted model, for example,

– residuals against fitted values;

– response against fitted values;

– normal QQ-plot of standardized residuals (useful for normal models).

The use of the optional argument ask=TRUE in plot is recommended, since otherwise all
plots are made consecutively, one after the other, and one can only examine the last one on
the screen. (An alternative, of course is to specify the appropriate par(mfrow=c(.,.)).)
If we want to make different plots we need more than the result of printing the output
of glm or its summary. Since the glm function produces an object of class glm, we can
extract much more information from it than by just printing it. It contains also fitted
values, residuals, deviances, and other useful information. From the glm object [Link]
we can, for instance, obtain the vector of fitted expected responses by

> [Link]$fitted

while

> [Link]$residuals

yields the vector of unscaled residuals. The diagonal elements of the working weight matrix
W
c at the last iteration are obtained using

> [Link]$weights

Check the R help page for glm to see what else can be extracted from a glm object.
Chapter 5

Time series in R

5.1 Basics
The ts data structure in R is meant for storing times series observations. A ts object can be
seen as a matrix with several attributes: start, end, deltat, frequency, .... The columns of
the matrix correspond to different variables, i.e., to different time series that together form
one multivariate time series; each row represents observations at a given time, the rows
being ordered consecutively in time. If a univariate time series is considered, the matrix
reduces to a vector which is, confusingly, sometimes printed by R as if it were a matrix.
The attribute start specifies the time of the first observation, and it corresponds to the
observation time of the first row; the attribute end denotes the time of the last observation
and corresponds to the observation time of the last row. The attribute frequency is the
sampling rate (number of observations per time unit), and deltat denotes the interval be-
tween observation times, which is the reciprocal of frequency. Only one of the attributes
frequency and deltat should be provided. The attribute frequency is used when the
series are sampled an integral number of times in each unit time interval. For instance, one
could use a value of 7 for frequency when the data are sampled daily and the natural time
period is a week, or 12 when the data are sampled monthly and the natural time period is
a year. Values of 4 and 12 are automatically assumed to imply a quarterly and monthly
series, respectively. The current version of R only deals with regularly spaced time series.
A ts object can be created with the R function ts: its argument data contains the vector
or matrix of observations, and the appropriate attributes form the other arguments:

> mydata=cbind(rnorm(26), rnorm(26)) #matrix of bivariate time series

> myts=ts(mydata,start=c(2009,6), freq=12) #ts object of bivariate
time series

Names can be given to the different series via the argument names. Within one ts object
multiple series can be stored only if they have the same number of observations taken at
the same time points. The functions [Link] and plot applied to an object of class ts
provide an appropriate plot of a time series. To see this, compare the plots that result
from the two commands

> [Link](mydata)
> [Link](myts)

18
5.1. Basics 19

The R function acf can be used to compute and plot the sample autocorrelation, au-
tocovariance, and partial autocorrelation for a time series over a range of lags starting at
lag 0. The autocorrelation is the default setting, while the autocovariance and partial au-
tocorrelation are obtained using the argument type. To obtain the partial autocorrelation
the function pacf can also be used. In a plot of the (partial) autocorrelation function,
horizontal bands are placed at heights 1.96/length(x)1/2 , which are the approximate 95%
confidence limits around zero. By default, the acf function produces a plot; this can be
suppressed using the optional argument plot=F. One can also specify the maximum lag
at which to compute the function value, using the optional argument [Link]. Otherwise
the maximum lag is computed based on an internal formula. Some examples:

> acf(myts) # plot of the autocorrelation

> pacf(myts) # plot of the partial autocorrelation
> acf(myts, type="covariance", [Link]=25) # plot of the autocovariance
> acf(myts, type="covariance", plot=F) # outputs numerical values

For multivariate series the function acf not only produces autocorrelations, but also the
cross correlations of the constituting univariate series, i.e. if {Xt } and {Yt } are two such
univariate series, acf also gives the sample version of

γXY (h) = E[(Xt − EX)(Yt+h − EY )].

Note that in general γXY (h) 6= γY X (h). The R function ccf produces the cross correlation
function of two univariate time series.
The R function diff applied to a time series produces differences of specified order
and lag:

> diff(y,lag=2, differences=1) # ∇2 y

> diff(y,d=3) # ∇3 y
> diff(diff(x,lag=12),2) # ∇2 ∇12 y

The function filter produces a new time series that is a filtered version of the time
series to which it is applied. It requires two arguments: a univariate or multivariate
time series (or a vector, or a matrix with each column representing a univariate compo-
nent of a multivariate series) x, and an argument filter, which defines the vector of
filter coefficients—the first coefficient of the filter applies to the last time point in the
filter window, and the last coefficient applies to the earliest time point in the filter win-
dow. To use a moving average filter, use the option method="convolution". Otherwise,
use method="recursive" to use autoregression. For the convolution method, the option
sides=1 results in a filtered time series depending only on the past and present values of
the original time series. The default sides=2 centers the filter coefficients at the present,
so that the resulting filtered time series depends on the past, present and future values
of the original time series in a symmetric way, and hence an odd number of coefficients
should be provided. The function filter produces a time series of the same length as the
20 Chapter 5. Time series in R

original but with missing values for coordinates that cannot be computed at the begin-
ning and end of the time series; the output is a ts object. As an illustration compare the
following outputs for x = (1, 2, 3, 4):

> filter(x, filter=c(1,2,3)) # x[t+1] + 2x[t] + 3x[t-1]

Time Series:
Start = 1
End = 4
Frequency = 1
[1] NA 10 16 NA
> filter(x, filter=c(1,2,3), sides=1) # x[t] + 2x[t-1] + 3x[t-2]
Time Series:
Start = 1
End = 4
Frequency = 1
[1] NA NA 10 16
> filter(x, filter=c(2,3), method="recursive")
Time Series:
Start = 1
End = 4
Frequency = 1
[1] 1 4 14 44

Note that the last command also computes

xflt = x[t] + 2x[t-1] + 3x[t-2]

but now this is done recursively, i.e., only the present value is taken from the time series
x, whereas the past values are taken from the already computed part of the filtered time
series, and initial values of zero are used: x[-1]=0, x[0]=0. If we want to use initial
values different from zero, we include the optional argument init.

> # now use x[-1]=10, x[0]=0

> filter(x, filter=c(2,3), method= "recursive", init=c(0,10))

5.2 ARIMA models

Several R functions deal with simulation and estimation of ARMA or, more generally,
ARIMA models. In these functions the argument model is a list that specifies the type of
model. Some simple examples for models are

> mod1 <- list(ar=c(0.5, 0.2)) # AR(2)

> mod2 <- list(ma=c(0.5)) # MA(1)
> mod3 <- list(ar=0.2, ma=c(0.5,0.2)) # ARMA(1,2)
5.2. ARIMA models 21

The vectors ar and ma specify the coefficients of the polynomials in the AR and MA part
of the model. Note that R automatically adds the leading constant β0 = 1.
A simulation of an AR(I)MA process can be performed by means of the R function
[Link]. Its required arguments are the model and the length of the series n. The
argument [Link] is optional (default: rnorm). It specifies the distribution of the white
noise process {Zt } in the definition of the ARIMA process. Here are three examples:

> ts1 <- [Link](model=mod1, n=500, sd=2) # AR(2) with N(0,4) white noise
> ts2 <- [Link](model=mod2, n=500, [Link]=rt, df=5) # MA(1) with
# t(5) white noise
> ts3 <- [Link](model=mod3, n=100) # ARMA(1,2) with N(0,1) white noise

The output of the [Link] function is a ts object.

The parameters of ARIMA processes can be estimated with the R function arima.
The first argument is a ts object, while the second argument is the order (p, d, q) of the
ARIMA model to be estimated. The command

> [Link] <- arima(ts3, order = c(1,0,2))

will fit an ARMA(1,2) model to the time series ts3. The default estimation method uses
conditional sum-of-squares to obtain the initial estimate, then maximum likelihood to find
an optimal estimate. To force the function to use maximum-likelihood only, include the
optional argument method="ML". The output of arima includes the estimated AR and
MA coefficients with the intercept and the standard error on each estimate, an estimate of
the variance of the innovations (the Zt in the lecture notes), the maximized log-likelihood,
and the AIC for the fit:

> [Link]

Call:
arima(x = ts3, order = c(1, 0, 2))

Coefficients:
ar1 ma1 ma2 intercept ## compare with coefficients in
0.1174 0.5143 0.3833 -0.0568 ## mod3 used for simulating ts3
s.e. 0.3073 0.2931 0.2018 0.2205

sigma^2 estimated as 1.069: log likelihood = -145.51, aic = 301.02

The intercept is the mean of the time series. If you want to force the mean estimate to
equal zero, insert the setting [Link]=F within the arima function. We see that σ 2 ,
which was equal to 1 in the simulation of ts3, is estimated by 1.069.
To view a plot of the residuals Rt , use

> [Link](residuals([Link]))
22 Chapter 5. Time series in R

Whereas arima uses the exact likelihood, the R function ar can be used to fit an
autoregressive (AR) model by solving the Yule–Walker equations. If the optional setting
aic = FALSE is included, one should also include the setting [Link] = p, where p
is the highest order that should be considered. The fitted model includes the estimated
coefficients and an estimate of the (innovations) variance σ 2 . For example,

> ts4 <- [Link](n=500, model=list(ar=c(0.1,-0.2,0.3)))

> [Link] <- ar(ts4, aic=FALSE, [Link]=3)
> [Link]

Call:
ar(x = ts4, aic = FALSE, [Link] = 3)

Coefficients:
1 2 3
0.0551 -0.2335 0.2306 ## compare these with the true coefficients above

Order selected 3 sigma^2 estimated as 0.9584

The fitting method can be changed to maximum-likelihood or ordinary least-squares using

the optional arguments method="mle" or method="ols", respectively.
To perform a diagnostic check on the goodness of the fit, we may examine a plot of
the residuals or a correlogram of the sample autocorrelation of the residuals to determine
whether they resemble white noise:

> plot([Link]$resid) # plot the residuals;

> acf([Link]$resid) # sample autocorrelation of the residuals

We may also perform a Portmanteau test on the residuals, using either the Box-Pierce or
Ljung-Box test statistic, provided the fit was obtained using arima (The function ar does
not provide residuals).

> [Link]([Link]$resid, fitdf=3, lag=10) # Box-Pierce test

> [Link]([Link]$resid, type="Ljung", fitdf=3, lag=10) # Ljung-Box test

The setting fitdf should be set to equal p + q, while the lag setting should be larger than
p + q, but no larger than n.
To predict the next n values of the time series based on the fitted AR model, use
the generic R function predict with the fitted AR model as the first argument and the
parameter [Link]:

> predict([Link], [Link]=10)

This produces a ts object containing predictions for the next 10 values of the original time
series, plus a time series object containing the standard errors on those 10 predictions.
The function predict will also predict the next n values based on a fitted ARMA or
ARIMA model, obtained with arima. The R command is the same:
5.2. ARIMA models 23

> predict([Link], [Link]=20)

This produces the next 20 predicted values of the original time series, with the standard
errors.
Chapter 6

Some literature on R

P. Dalgaard, Introductory Statistics with R, Springer, 2008.

J.J. Faraway, Linear Models with R, Chapman & Hall/CRC, 2009.

J.J. Faraway, Extending the Linear Model with R: Generalized Linear, Mixed Effects and
Nonparametric Regression Models, Chapman & Hall/CRC, 2005.

C. Ritz, J.C. Streibig, Nonlinear Regression with R, Springer, 2008.

P.S.P. Cowpertwait, A.V. Metcalfe, Introductory Time Series with R, Springer, 2009.

24