r - 在 lme 中初始化协方差结构

Question

如何为以下模型初始化非结构化协方差矩阵？

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))

我尝试了一些变体，但每次出现错误时我都会得到:

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

除了治疗平均效应(4 个参数)、随机效应方差 (1) 和残差方差 (1) 之外，要将具有 6 个参数的非结构化相关矩阵拟合到只有 16 个参数的数据集，实际上是很困难的点。如果我尝试使用更大的随机版本的数据集，它工作正常。

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))

现在我们可以进行实验，看看我们可以处理多小的数据集。将上面的东西打包成一个模拟数据并尝试拟合的测试函数，如果拟合失败则返回 TRUE:

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")
}

尝试测试函数N次并报告失败的比例:

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

针对一定数量的主题尝试一下(慢):

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

结果:

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

令我有些惊讶的是，对于 4 个主题，这将在三分之一的时间内起作用； 96% 的时间有 5 个科目:并且总是有 >5 个科目。