除了处理平均效应(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 个科目。