线性混合模型中畸形拟合(Singular fit)的判断问题

线性混合模型中畸形拟合(Singular fit)的判断问题

在线性混合模型中,经常遇到畸形拟合的问题,如下:

> ####### singular info
> Model1 = lmer(data = DF, Y ~ A * B + (1 +  A |Sub) )
boundary (singular) fit: see ?isSingular

这时可以通过函数isSingular()来判断:

> ####### singular info
> Model1 = lmer(data = DF, Y ~ A * B + (1 +  A |Sub) )
boundary (singular) fit: see ?isSingular
> isSingular(Model1)
[1] TRUE

但有时(情况较少)建模时输出的信息和通过函数的结果不一致:

> Model2 = lmer(data = DF,Y ~ A * B + (1 +  A |Sub) + (1|Item))
boundary (singular) fit: see ?isSingular
> isSingular(Model2)
[1] FALSE

这是因为isSingular()中还有一个参数tol(tolerance, 容忍度):

isSingular(x = Model2,tol = 1e-05)

容忍度默认为0.00001,但建模时的容忍度要低一些。

这里还要说明的是,如果建模中显示畸形拟合,会在模型信息中储存:

> Model1@optinfo$conv$lme4$messages
[1] "boundary (singular) fit: see ?isSingular"
> Model2@optinfo$conv$lme4$messages
[1] "boundary (singular) fit: see ?isSingular"

因此是否畸形拟合,可以从模型信息中获取,也可以通过函数判断来获取。

我模拟了从0.00001到0.0001(步长为0.00001)不同容忍度下的对含有两个固定因子(包含交互作用)和两个随机因子的所有可能模型,它们的模型信息和函数检测的信息之间的一致性。

结果如下:

> TestFit
# A tibble: 10 x 2
   tolerance consistency
       <dbl>       <dbl>
 1  0.00001        0.971
 2  0.00002        0.971
 3  0.00003        0.971
 4  0.00004        1    
 5  0.00005        1    
 6  0.00006        1    
 7  0.00007        1    
 8  0.00008        1    
 9  0.00009        1    
10  0.0001         1    

看出,容忍度不小于0.00004时,模型信息和函数检测的结果就完全一致了,因此用户如果在用isSingular()判断是否畸形拟合时要小心。

后来我分别测了容忍度为0.0001,0.001,和0.01时,一致性都为1。从这个结果来看,将容忍度设为0.0001或许是可行的。


  • 数据信息:

执行下面信息,可获取数据信息:

DF = read_csv('https://raw.githubusercontent.com/usplos/Eye-movement-related/master/DataforShiny.csv')
DF[c('A','B')] = lapply(DF[c('A','B')], factor) # 变量因子化
  • 计算一致性的代码:
ModelInfo = LMMRun_Parallel(df = DF,DV = 'Y',IV = c('A','B'),Cluster = c('Sub','Item'),Ifrun = T,output = 'AA',Family = NULL) # 计算所有可能的模型
ModelInfoT = ModelInfo %>% arrange(-Converge,Singular,-Nchar) %>% select(-c(R2.M,R2.C,AIC,BIC)) # 对上一步的结果筛选变量

TestFit = tibble(tolerance = double(length = 10),
                    Fittedness = double(length = 10))
k = 1
for(tt in seq(1,10,1)/1e5){
  print(tt)
  ModelWarningSingular = logical(length = nrow(ModelInfoT)) # 储存模型信息
  ModelTestSingular = logical(length = nrow(ModelInfoT)) # 储存函数检测信息
  
  for(rr in 1:nrow(ModelInfoT)){
    print(ModelInfoT$formula[rr] %>% as.character())
    ModelT = lmer(data = DF, 
                  formula = ModelInfoT$formula[rr] %>% as.character() %>% as.formula())
    Message = ModelT@optinfo$conv$lme4$messages # 模型信息
    ModelTestSingular[[rr]] = isSingular(ModelT,tol = tt) # 函数检测信息
    if(length(Message) == 0){
      ModelWarningSingular[rr] = F
    }else if(str_sub(Message[[1]],start = 1,end = 8) == 'boundary'){
      ModelWarningSingular[rr] = T
    }else{ # 这种情况是模型无法收敛,但不出现畸形拟合
      ModelWarningSingular[rr] = F
    } 
  }
  print(k)
  TestFit$tolerance[k] = tt
  TestFit$Fittedness[k] = mean(ModelWarningSingular == ModelTestSingular)
  k=k+1
  Sys.sleep(2)
  cat('\014')
}
names(TestFit)[[2]] = 'consistency'
TestFit

编辑于 2019-10-27

文章被以下专栏收录