R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

(g)lmer on epäluotettava?

Zhang et al. ovat julkaisseet vastikään artikkelin, jossa he vertaavat eräiden tilasto-ohjelmistojen sekamallien sovittamiseen kehitettyjä ominaisuuksia. Artikkelissa esitetyt tulokset eivät ylipäätään ole kovin mairittelevia, sillä parhaiksi havaittujenkin ohjelmistojen (SAS:n NLMIXED) tyypin I -virhetaso on noin kaksinkertainen ilmoitettuun verrattuna. Erityisen huonosti tuntuu kuitenkin pärjäävän R:n laajennuspaketin lme4 funktion glmer(). Lisäksi R sai pyyhkeitä myös epäluotettavasta mallin paramatrien arvioinnista.

Yritetäänpä toistaa artikkelin tulokset ajantasaisella 32-bittisellä R:llä (R 2.14.2., lme4 0.999902344-0, Matrix 1.0-4, splines 2.14.2). Asiasta on keskusteltu myös ainakin kahdessa viestiketjussa (1, 2) R-sig-ME -listalla.

1. Davisin aineiston analyysi

Davisin aineiston analyysi voidaan toistaa R:ssä seuraavalla koodilla:

# HSAUR contains respiratory data
data("respiratory", package = "HSAUR")
 
# Correct an error in the data
respiratory$status[rownames(respiratory)=="428"]<-"poor"
 
# Convert pretest to covariate "baseline"
resp <- subset(respiratory, month &gt; "0" )
resp$baseline <- rep(subset(respiratory, month == "0")$status, each=4)
 
# Invert the factor levels to correspond to Davis et al.
resp$treatment2<-factor(resp$treatment, levels=c("treatment","placebo"))
resp$sex2<-factor(resp$sex, levels=c("male", "female"))
 
# Fit the models
library(lme4)
 
# Laplace approximation
print(m_glmer_4.L <- glmer(status ~ centre + treatment2 + sex2 + age + baseline + (1|subject), family=binomial,data=resp), cor=FALSE)
 
# Gauss-Hermite approximation
print(m_glmer_4.L <- glmer(status ~ centre + treatment2 + sex2 + age + baseline + (1|subject), family=binomial,data=resp, nAGQ=16), cor=FALSE)

Artikkelissa esitetyt SAS:n NLMIXED-ohjelmalla saadut mallin parametrien estimaattein arvot ja tässä R:llä lasketut ovat hyvin lähellä toisiaan, vakiotermiä lukuunottamatta. Tässä lasketut R:llä tuotetut parametrien estimaattien arvot eroavat selvästi artikkelissa esitetyistä, eikä ole vielä selvillä, miten ne oikeastaan on laskettu. On toki myönnettävä, että esimerkiksi iän vaikutuksen estimaatti on R:llä laskettuna noin 33% pienempi kuin SAS:n NLMIXED-ohjelmalla laskettu estimaatti, joten suhteellinen virhe voi pienillä estimaatin arvoilla olla suurikin.

            SAS     R
Intercept   2.25    0.70
Center      0.97    0.95
Treatment  -2.10   -2.18
Sex        -0.24   -0.24
Age        -0.03   -0.02
Baseline    2.94    2.94

R:ssä funktio glmer() tuottaa hyvin samanlaiset estimaatit, käyttipä mallin sovittamiseen Laplace- tai Gauss-Hermian -menetelmää.

2. Simulaatiot

Otetaan tähän esimerkiksi yksi simulaatio, jossa mallin molempien parametrien estimaattien (Beta) pitäisi olla 1, koska simulaatio on toteutettu siten. Oletetaan kovarianssimatriisin (Cor) luomisessa käytettävän parametrin (Tau) arvoksi simulaatioissa 2, ja ryhmäkohtaisen havaintojen (N) lukumääräksi 500. Toistetaan aineiston simulointi ja mallin sovitus 1000 kertaa (R), ja vedetään toistoista saadut tulokset lopuksi yhteen. Simulaatioissa käytetty koodi on seuraava (muokattu aiemmin mainituista R-sig-ME:n viestiketjuista, erityisesti Ben Bolker ja Dave Fournier):

# Parameters of the simulation
N<-500
Beta<-c(1,1)
Tau<-2
Cor<-matrix(c(1,0.25,0.25,1),nrow=2)
R<-1000
 
# Simulation function for the data
library("lme4")
simfun <- function(n,beta=c(1,1),tau=2.0, Cor=matrix(c(1,0.25,0.25,1),nrow=2)) {
    D <- expand.grid(id=factor(1:n),t=1:3)
    D$x <- rnorm(3*n)  # x_{it} ~ N(0,1)
    X <- model.matrix(~x,data=D)
    Z <- model.matrix(~id-1+id:x,data=D)
    b <- MASS::mvrnorm(n,mu=c(0,0),Sigma=tau^2*Cor)
    D$Y <- rbinom(nrow(D), prob=plogis(X %*% beta + Z %*% c(b)), size=1)
    D
}
 
# Simulations
m1<-matrix(ncol=2, nrow=R, data=NA)
for(i in 1:R) {
  set.seed(i)
  dat=simfun(N)
  fm <- glmer(Y~  x + ( x|id), data=dat,family=binomial)
  m1[i,] <- cbind(as.vector(fixef(fm))[2], as.vector(sqrt(vcov(fm,useScale = FALSE)[2,2])))
}
 
# Calculate the confidence intervals based on normal approximation
cilower<-m1[,1]-m1[,2]*1.96
ciupper<-m1[,1]+m1[,2]*1.96
 
# How many times is the known parameter outside of the confidence interval?
sum(cilower>1 | ciupper<1)/R

Artikkelissa SAS:n NLMIXED-ohjelmalle ilmoitettu tyypin I-virheen frekvenssi oli noin 5-10%, ja liki samaan tulokseen on ADMB-työkalun osalta päätynyt Dave Fournier. Yllä suoritetun simulaation perusteella R:n glmer()-funktion tyypin I-virheen frekvenssi sen sijaan todella näyttää olevan paljon korkeampi, noin 25%.

Mitä tästä sitten pitäisi päätellä? Ainakin se voitaneen varmuudella sanoa, että glmer()-funktion antamat luottamusvälit ja p-arvot ovat liian optimistisia, mikä johtaa paljon oletettua korkeampaan väärien positiivisten tulosten määrään. Funktion glmer() antamien p-arvojen tulkinnassa kannattanee siis olla hieman varovainen.


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *