[R] User-specified correlation structure (e.g., 2-banded Toeplitz)
Reid Landes
rdlandes at gmail.com
Fri Feb 8 16:39:42 CET 2008
Dear All:
I am trying to fit a special case of a 2-banded Toeplitz correlation
structure. A 2-banded Toeplitz has ones on the diagonal, a
correlation, RHO1, on the first off-diagonal, and a correlation, RHO2,
on the second off-diagonal, with zeros on all subsequent
off-diagonals. After reading relevant sections in Mixed-Effects
Models in S and S-PLUS (Pinheiro & Bates, 2000) and searching on the
R-help archives, I've figured out how to get the 2-banded Toeplitz,
but not the desired special case.
In the example below, the initial value RHO1 = 0 and RHO2= -0.3. The
output matrix below is an example of the special case I'd like to fit
-- a 2-banded Toeplitz constraining RHO1=0.
Fitting the 2-banded Toeplitz structure to the ``Orthodont'' example
dataset provided in R-Help, we estimate RHO1 also (since the example
matrix below contains INITIAL values).
---------------------Start R-code & output -------------------------------
#This intilizes a 2-banded Toeplitz structure
cs1ARMA <- corARMA(value = c(0,-.3), form = ~ 1 | Subject, p = 2, q = 0)
cs1ARMA <- Initialize(cs1ARMA, data = Orthodont)
corMatrix(cs1ARMA)$M01
[,1] [,2] [,3] [,4]
[1,] 1.0 0.0 -0.3 0.0
[2,] 0.0 1.0 0.0 -0.3
[3,] -0.3 0.0 1.0 0.0
[4,] 0.0 -0.3 0.0 1.0
> TOEP2 <- gls(distance ~ Sex * I(age - 11), Orthodont,
+ correlation = corARMA(value = c(0,-.3), form = ~
1 | Subject, p = 2, q = 0),
+ weights = varIdent(form = ~ 1 | age))
#-- Selected output follows-----
Correlation Structure: ARMA(2,0)
Formula: ~1 | Subject
Parameter estimate(s):
Phi1 Phi2
0.3269544 0.4897645
--------------------- End R-code & output --------------------------------
I cannot figure out how to restrict RHO1 = 0, while allowing
estimation of RHO2.
Maybe an answer lies in specifying a different ``position vector''
other than the default: corARMA(..., form = ~ 1 | Subject ...). (See
p226 of Pinheiro & Bates, 2000 for explanation of a position vector.)
But I'm not totally sure if I understand the position vector and I
know I don't know how it works in R.
Then again, there is likely a completely different way to solve this problem.
Any help will be appreciated!
Reid
More information about the R-help
mailing list