[R-sig-ME] Custom nlme::CorStruct Class based on CorCompSymm

[Mychaleckyj, Josyf C. (Joe) (jcm6t)]

I am working on a custom CorStruct class based on CorCompSymm. There will be a single parameter passed into the class for optimization, like rho.The grouping will be ~ 1 | Subject. Initially i plan to use this with gls but ultimately would like to extend to lme also.

The complication in this case is that the calculation of the correlation matrices requires another fixed parameter that varies between pairs of repeated observations within each subject - this will be passed into the functions as a numerical lower.tri() vector.  Within corMatrix I need to be sure that I am applying the correct parameter value from the vector to the correct pair of observations for the within-Subject correlation matrix.

For example: in a Subject with 3 observations, fixparam[1,2] = fixparam[2,1] = 0.25; fixparam[1,3] = param[3,1] = 0.125, etc.
(this is a genetics application but I’m keeping genetics out of the discussion).

From P&B discussion of groupedData I expected the Subjects to be default ordered as max within-group response. This is true in getGroups but my simple tests using the Orthodont data suggest that the subject list and observations are simply ordered by first occurrence in the data.frame of observations as exposed within CorMatrix. I can’t determine the ordering within Subject from my test case below but assuming as per df order.

I’ve posted a scrambled unbalanced test case below to show my reasoning. My specific questions:

1. Am i correct in my inference as to how Subjects and observations are accessed in the CorCompSymm class on which I am basing my custom class ? Hopefully the  print method is not reordering under the covers.

2. Can i rely on this default ordering to be preserved ?

3. I would feel better if there were some way to specify the rowname or a unique identifier for each observation to index the fixed parameter vector, but not alter the RE structure with a covariate. Any hidden tricks ?


Thanks for your help,

Joe.


# Test case:
# 1.  Create a simple unbalanced scrambled subset of Orthodont to test for precedence

> Orthodont.sub
Grouped Data: distance ~ age | Subject
    distance age Subject    Sex
11      24.0  12     M03   Male
106     25.0  10     F11 Female
105     24.5   8     F11 Female
107     28.0  12     F11 Female
80      26.5  14     F04 Female
78      24.5  10     F04 Female
77      23.5   8     F04 Female
25      22.0   8     M07   Male
10      22.5  10     M03   Male
108     28.0  14     F11 Female

> table(Orthodont.sub$Subject)
M07 M03 F04 F11
  1   2   3   4

# as expected based on P & B
> getGroups(Orthodont.sub)
 [1] M03 F11 F11 F11 F04 F04 F04 M07 M03 F11
Levels: M07 < M03 < F04 < F11

> cs4<- corCompSymm(value = 0.3,form = ~ 1 | Subject)
> cs4<- Initialize(cs4, data = Orthodont.sub)

# Subject order in list M03, F11, F04, M07
> corMatrix(cs4)
$M03
     [,1] [,2]
[1,]  1.0  0.3
[2,]  0.3  1.0

$F11
     [,1] [,2] [,3] [,4]
[1,]  1.0  0.3  0.3  0.3
[2,]  0.3  1.0  0.3  0.3
[3,]  0.3  0.3  1.0  0.3
[4,]  0.3  0.3  0.3  1.0

$F04
     [,1] [,2] [,3]
[1,]  1.0  0.3  0.3
[2,]  0.3  1.0  0.3
[3,]  0.3  0.3  1.0

$M07
     [,1]
[1,]    1

# Now check that the same thing happens with Orthodont coerced to simple data frame
> Orthodont.sub<-as.data.frame(Orthodont.sub)
> Orthodont.sub
    distance age Subject    Sex
11      24.0  12     M03   Male
106     25.0  10     F11 Female
105     24.5   8     F11 Female
107     28.0  12     F11 Female
80      26.5  14     F04 Female
78      24.5  10     F04 Female
77      23.5   8     F04 Female
25      22.0   8     M07   Male
10      22.5  10     M03   Male
108     28.0  14     F11 Female

# double check
> table(Orthodont.sub$Subject)
M07 M03 F04 F11
  1   2   3   4

# all appears to be the same as cs4 - does not depend on groupedData object
> cs5<- corCompSymm(value = 0.3,form = ~ 1 | Subject)
> cs5<- Initialize(cs5, data = Orthodont.sub)
> corMatrix(cs5)
$M03
     [,1] [,2]
[1,]  1.0  0.3
[2,]  0.3  1.0

$F11
     [,1] [,2] [,3] [,4]
[1,]  1.0  0.3  0.3  0.3
[2,]  0.3  1.0  0.3  0.3
[3,]  0.3  0.3  1.0  0.3
[4,]  0.3  0.3  0.3  1.0

$F04
     [,1] [,2] [,3]
[1,]  1.0  0.3  0.3
[2,]  0.3  1.0  0.3
[3,]  0.3  0.3  1.0

$M07
     [,1]
[1,]    1






	[[alternative HTML version deleted]]