Initialise covariance structure in lme

Question

How can I initialise a unstructured covariance matrix for the following model?

y<-data.frame(response=c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17),
               treatment=factor(rep(1:4,4)),
               subject=factor(rep(1:4,each=4))
               )
fit<-lme(response~-1+treatment,y,random=~1|subject,
         correlation=corSymm(form=~1|subject))

I tried some variants but I get every time I get the error:

Error in lme.formula(response ~ -1 + treatment, y, random = ~1 |  : 
  nlminb problem, convergence error code = 1
  message = function evaluation limit reached without convergence (9)

Answer 1

It's practically difficult to fit an unstructured correlation matrix with 6 parameters in addition to a treatment mean effect (4 parameters), a random-effects variance (1), and a residual variance (1) to a data set with only 16 points. If I try with a larger, randomized version of your data set, it works fine.

nSubj <- 20
respVec <- c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17)
set.seed(101)
y<-data.frame(response=sample(respVec,size=4*nSubj,replace=TRUE),
               treatment=factor(rep(1:4,nSubj)),
               subject=factor(rep(1:nSubj,each=4))
               )
library(nlme)
fit<-lme(response~-1+treatment,y,random=~1|subject,
         correlation=corSymm(form=~1|subject),
         control=lmeControl(msVerbose=TRUE))

Now we can experiment and see how small a data set we can get away with. Package the stuff above into a test function that simulates data and tries a fit, returning TRUE if the fit fails:

testFun <- function(nSubj) {
    y<-data.frame(response=sample(respVec,size=4*nSubj,replace=TRUE),
               treatment=factor(rep(1:4,nSubj)),
               subject=factor(rep(1:nSubj,each=4))
               )
    fit <- try(lme(response~-1+treatment,y,random=~1|subject,
         correlation=corSymm(form=~1|subject)),silent=TRUE)
    inherits(fit,"try-error")
}

Try the test function N times and report the proportion of failures:

testFun2 <- function(nSubj,N) {
   mean(replicate(N,testFun(nSubj)))
}

Try it out for a range of numbers of subjects (slow):

set.seed(101)
testRes <- sapply(4:20,testFun2,N=50)

Results:

##  [1] 0.64 0.04 0.00 0.00  ... 0.00

Somewhat to my surprise, this will work a third of the time with 4 subjects; 96% of the time with 5 subjects: and always with >5 subjects.