[R] My own correlation structure with nlme
Spencer Graves
spencer.graves at pdf.com
Sun Nov 26 01:45:30 CET 2006
I haven't seen a reply to this post, so I'll offer some comments,
even though I can't answer the question directly.
1. Given your question, I assume you have consulted Pinheiro and
Bates (2000) Mixed-Effects Models in S and S-Plus (Springer). If you
have not, I'm confident you will find material relevant to your question
there, especially in chapters 6-8.
2. Are you aware that the standard R installation includes a
subdirectory "~R\library\nlme\scripts", which contain copies of R
commands to create all the analyses in the book? In particular,
"ch06.R" and "ch08.R" might be particularly relevant to your question.
If you have not made local copies of these and walked through the code
line by line, I suggest you do so. I've learned a lot doing that.
3. Which versions of R and 'nlme' are you using? Some minor
enhancements to the help files were added a few months ago, and there
might be something helpful in the current examples that wasn't there
before.
4. What have you done to develop simple toy examples to help you
test your understanding of different aspects of the code? This
technique has helped me solve many problems.
5. Are you familiar with the use of the "methods" command to
trace logic flow? Consider for example the following:
> methods("corMatrix")
[1] corMatrix.corAR1* corMatrix.corARMA* corMatrix.corCAR1*
[4] corMatrix.corCompSymm* corMatrix.corIdent* corMatrix.corNatural*
[7] corMatrix.corSpatial* corMatrix.corStruct* corMatrix.corSymm*
[10] corMatrix.pdBlocked* corMatrix.pdCompSymm* corMatrix.pdDiag*
[13] corMatrix.pdIdent* corMatrix.pdMat* corMatrix.reStruct*
Non-visible functions are asterisked
6. Are you familiar with using "getAnywhere" to obtain code that
may be hidden with namespaces? For example, I just learned that
'corAR1' and 'corMatrix.corAR1' are two distinct functions. I found the
latter with this "methods" command, and I got the code to it using
"getAnywhere".
7. Are you familiar with the 'debug' command (in the 'base'
package, not the 'debug' package, which is probably much more powerful
but I haven't used the latter)? This allows a user to trace through
code line by line, examining the contents of objects -- and even
modifying them if you want. This is an extremely powerful way to learn
what a piece of code does.
8. If the above does not produce the answers you seek with a
reasonable additional effort, please submit another post. To increase
your chances of getting replies that are quicker and more helpful than
this one, please include commented, minimal, self-contained,
reproducible code. I'm confident that I could have helped you more with
less effort if your 'pairCorr' code had been included a sample call that
I could copy from your email into R and see if I get the same error
message you got. If I failed to get the same error as you saw, that
suggests a problem in your installation. If I got the same error
message, there is a good chance that I figure out how to get around that
and provide more helpful suggestions in less time than I've been
thinking about this.
Hope this helps.
Spencer Graves
Mohand wrote:
> Dear all,
>
> I am trying to define my own corStruct which is different from the
> classical one available in nlme. The structure of this correlation is
> given below.
> I am wondering to know how to continue with this structure by using
> specific functions (corMatrix, getCovariate, Initialize,...) in order to
> get a structure like corAR1, corSymm which will be working for my data.
>
> Thanks in advance.
>
> Regards
>
> M. Feddag
>
>
>
>
> *Correlation structure _
> _*
>
>
> pairCorr <- function(A, B, lowerTriangle = TRUE){
> n <- length(A)
> m <- length(B)
> if (n != m) stop("A and B must have same length")
> result <- matrix(0, n, n)
> same.A <- outer(A, A, "==")
> same.B <- outer(B, B, "==")
> A.matches.B <- outer(A, B, "==")
> B.matches.A <- outer(B, A, "==")
> result <- result + 0.5*same.A + 0.5*same.B -
> 0.5*A.matches.B - 0.5*B.matches.A
> rownames(result) <- colnames(result) <- paste(A, B, sep = "")
> if (lowerTriangle) result <- result[row(result) > col(result)]
> result
> }
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
More information about the R-help
mailing list