I recently took some time to figure out how to write a new method for nlme to enable structuring the variance-covariance matrix of the random effects in a specific way. My goal here was to be able to run Dave Kenny's social relations model (Kenny, 1994) using multilevel modeling and the approach described by Snijders and Kenny (1999). Taking this approach requires "tricking" the software in a way through the use of dummy variables and constraints on the variance-covariance matrix.
Figuring out how to write a new method was more challenging than I had initially expected. There are many twists and turns in lme and it took quite a bit of time to reverse engineer the software to figure out what was going on. Unfortunately, there isn't great documentation on the web for this process.
As part of my process, I created my own replication of one of the existing methods--pdCompSymm. I went through and commented each part of the different functions that are called, explaining my interpretation of what is going on. As you can see, there are some places where I'm just off and don't really know what's going on. I also converted some of the C code in nlme for running pdCompSymm into R code (this is the pdFactor.pdCompSymm function).
In the end, I was able to figure out enough of it to succeed in my goal of creating a new method for the social relations model through multilevel modeling in R. You can find this on my github page. I've called it pdSRM and it has some comments at the top that explain how to use it.
One lesson learned from this is that it is challenging--but not impossible!--to specify a structure for the variance-covariance matrix using nlme that is not already in the generic methods that are provided. I also learned a ton about how lme is working behind the scenes. This took a bunch of time, but did pay off in the end.
