Por que o SAS PROC GLIMMIX me fornece pistas aleatórias MUITO diferentes do glmer (lme4) para um binômio glmm

12

Eu sou um usuário mais familiarizado com R e tenho tentado estimar inclinações aleatórias (coeficientes de seleção) há cerca de 35 indivíduos com mais de 5 anos para quatro variáveis ​​de habitat. A variável de resposta é se um local foi habitat "usado" (1) ou "disponível" (0) ("uso" abaixo).

Estou usando um computador com Windows de 64 bits.

No R versão 3.1.0, eu uso os dados e a expressão abaixo. PS, TH, RS e HW são efeitos fixos (distância medida e padronizada para os tipos de habitat). lme4 V 1.1-7.

str(dat)
'data.frame':   359756 obs. of  7 variables:
 $ use     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Year    : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ...
 $ ID      : num  306 306 306 306 306 306 306 306 162 306 ...
 $ PS: num  -0.32 -0.317 -0.317 -0.318 -0.317 ...
 $ TH: num  -0.211 -0.211 -0.211 -0.213 -0.22 ...
 $ RS: num  -0.337 -0.337 -0.337 -0.337 -0.337 ...
 $ HW: num  -0.0258 -0.19 -0.19 -0.19 -0.4561 ...

glmer(use ~  PS + TH + RS + HW +
     (1 + PS + TH + RS + HW |ID/Year),
     family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))

O glmer me fornece estimativas de parâmetros para os efeitos fixos que fazem sentido para mim, e as inclinações aleatórias (que eu interpreto como coeficientes de seleção para cada tipo de habitat) também fazem sentido quando investigo os dados qualitativamente. A probabilidade de log para o modelo é -3050,8.

No entanto, a maioria das pesquisas em ecologia animal não usa R porque, com dados de localização do animal, a autocorrelação espacial pode tornar os erros padrão propensos a erros do tipo I. Enquanto R usa erros padrão baseados em modelo, os erros padrão empíricos (também Huber-branco ou sanduíche) são os preferidos.

Embora R atualmente não ofereça essa opção (que eu saiba - POR FAVOR, me corrija se eu estiver errado), o SAS oferece - embora eu não tenha acesso ao SAS, um colega concordou em me emprestar seu computador para determinar se os erros padrão mudar significativamente quando o método empírico é usado.

Primeiro, desejamos garantir que, ao usar erros padrão baseados em modelo, o SAS produza estimativas semelhantes às de R - para ter certeza de que o modelo seja especificado da mesma maneira nos dois programas. Não me importo se são exatamente iguais - apenas semelhantes. Eu tentei (SAS V 9.2):

proc glimmix data=dat method=laplace;
   class year id;
   model use =  PS TH RS HW / dist=bin solution ddfm=betwithin;
   random intercept PS TH RS HW / subject = year(id) solution type=UN;
run;title;

Eu também tentei várias outras formas, como adicionar linhas

random intercept / subject = year(id) solution type=UN;
random intercept PS TH RS HW / subject = id solution type=UN;

Eu tentei sem especificar o

solution type = UN,

ou comentando

ddfm=betwithin;

Não importa como especificamos o modelo (e tentamos várias maneiras), não consigo que as inclinações aleatórias no SAS se assemelhem remotamente às saídas do R - mesmo que os efeitos fixos sejam semelhantes o suficiente. E quando quero dizer diferente, quero dizer que nem mesmo os sinais são os mesmos. A probabilidade de log -2 no SAS era 71344,94.

Não consigo carregar meu conjunto de dados completo; então eu criei um conjunto de dados de brinquedo apenas com os registros de três indivíduos. O SAS me fornece resultados em alguns minutos; em R demora mais de uma hora. Esquisito. Com este conjunto de dados de brinquedos, agora estou obtendo estimativas diferentes para os efeitos fixos.

Minha pergunta: alguém pode esclarecer por que as estimativas de declives aleatórios podem ser tão diferentes entre R e SAS? Existe algo que eu possa fazer no R ou SAS para modificar meu código para que as chamadas produzam resultados semelhantes? Prefiro alterar o código no SAS, pois "acredito" que meu R calcule mais.

Estou realmente preocupado com essas diferenças e quero chegar ao fundo desse problema!

Minha saída de um conjunto de dados de brinquedo que usa apenas três dos 35 indivíduos no conjunto de dados completo para R e SAS é incluída como jpegs.

Saída R Saída SAS 1 Saída SAS 2 Saída SAS 3


EDITAR E ATUALIZAR:

Como o @JakeWestfall ajudou a descobrir, as inclinações no SAS não incluem os efeitos fixos. Quando adiciono os efeitos fixos, eis o resultado - comparando inclinações R com inclinações SAS para um efeito fixo, "PS", entre programas: (Coeficiente de seleção = inclinação aleatória). Observe o aumento da variação no SAS.

R vs SAS para PS

Nova
fonte
Percebo que isso IDnão é um fator em R; verifique e veja se isso muda alguma coisa.
Aaron saiu de Stack Overflow
Vejo que você está ajustando os dois usando a Aproximação de Laplace para a probabilidade de log. Quais são as respectivas pontuações de probabilidade de log?
usεr11852 diz Reinstate Monic
1
Você verificou se está modelando a variável dependente na mesma direção?
Peter Flom - Restabelece Monica
1
A propósito, o que Peter está entendendo é que, por padrão, com dados binomiais rotulados como 0s e 1s, Rmodelará a probabilidade de uma resposta "1", enquanto o SAS modelará a probabilidade de uma resposta "0". Para tornar o modelo SAS a probabilidade de "1", você precisa escrever sua variável de resposta como use(event='1'). É claro que, mesmo sem fazer isso, acredito que ainda deveríamos esperar as mesmas estimativas das variações de efeito aleatório, bem como as mesmas estimativas de efeito fixo, embora com seus sinais revertidos.
Jake Westfall 10/09
1
@EricaN Uma coisa que você acabou de me lembrar é que você deve comparar os efeitos aleatórios de R com os do SAS usando a ranef()função em vez de coef(). O primeiro fornece os efeitos aleatórios reais, enquanto o segundo fornece os efeitos aleatórios mais o vetor de efeitos fixos. Portanto, isso explica muito por que os números ilustrados em sua postagem diferem, mas ainda há uma discrepância substancial que não posso explicar totalmente.
Jake Westfall

Respostas:

11

Parece que eu não deveria ter esperado que as pistas aleatórias fossem semelhantes entre os pacotes, de acordo com Zhang et al 2011. Em seu artigo Sobre a adaptação de modelos de efeitos mistos lineares generalizados para respostas binárias usando diferentes pacotes estatísticos , eles descrevem:

Abstrato:

O modelo de efeitos mistos lineares generalizados (GLMM) é um paradigma popular para estender modelos para dados de seção transversal a um ajuste longitudinal. Quando aplicado à modelagem de respostas binárias, pacotes de software diferentes e até procedimentos diferentes dentro de um pacote podem fornecer resultados bastante diferentes. Neste relatório, descrevemos as abordagens estatísticas subjacentes a esses diferentes procedimentos e discutimos seus pontos fortes e fracos quando aplicados para ajustar respostas binárias correlatas. Em seguida, ilustramos essas considerações aplicando esses procedimentos implementados em alguns pacotes de software populares a dados de estudo simulados e reais. Nossos resultados de simulação indicam uma falta de confiabilidade para a maioria dos procedimentos considerados, o que traz implicações significativas para a aplicação prática de tais pacotes de software populares.

Espero que o @BenBolker e a equipe considerem minha pergunta como um voto para o R incorporar erros padrão empíricos e a capacidade da Quadratura Gauss-Hermite para modelos com vários termos de inclinação aleatória para glmer, pois eu prefiro a interface R e gostaria de poder aplicar algumas análises adicionais nesse programa. Felizmente, mesmo que R e SAS não tenham valores comparáveis ​​para inclinações aleatórias, as tendências gerais são as mesmas. Obrigado a todos por sua contribuição, eu realmente aprecio o tempo e a consideração que você dedicou a isso!

Nova
fonte
desculpe: o que é um "erro padrão padrão"? Você quer dizer erros padrão de componentes de variação? Ou você quer dizer erros padrão em sanduíche?
precisa
desculpe ... significou SEs empíricos / sanduíche. Eu editei minha resposta.
Nova
@BenBolker Isso já foi incorporado?
Lepidopterist
Não. Eu continuo tentando descobrir como eu vou apoiar o desenvolvimento como esta, uma vez que não é tecnicamente parte do meu programa de pesquisa ...
Ben Bolker
4

Uma mistura de resposta e comentário / mais perguntas:

Eu ajustei o conjunto de dados 'brinquedo' com três opções diferentes de otimizador. (* Nota 1: provavelmente seria mais útil, para fins comparativos, criar um pequeno conjunto de dados subamostrando a partir de cada ano e ID, em vez de subamostrando as variáveis ​​de agrupamento. Como é, sabemos que o GLMM não funcionará particularmente bem com um número tão pequeno de níveis de variáveis ​​de agrupamento.Você pode fazer isso através de algo como:

library(plyr)
subdata <- ddply(fulldata,c("year","id"),
    function(x) x[sample(nrow(x),size=round(nrow(x)*0.1)),])

Código de montagem em lote:

Ntoy <- readRDS("Newton_toy.RDS")
library(lme4)
fitfun <- function(opt) {
    tt <- system.time(fit1 <- glmer(use ~  ps + th + rs + hw +
                                    (1 + ps + th + rs + hw |id/year),
                                    family = binomial, data = Ntoy,
                                    control=glmerControl(optimizer=opt),
                                    verbose=100))
    return(list(time=tt,fit=fit1))
}

opts <- c("nloptwrap","nlminbwrap","bobyqa")
## use for() instead of lapply so we can checkpoint more easily
res <- setNames(vector("list",length(opts)),opts)
for (i in opts) {
    res[[i]] <- fitfun(i)
    save("res",file="Newton_batch.RData")
}

Então eu li os resultados em uma nova sessão:

load("Newton_batch.RData")
library(lme4)

Tempo decorrido e desvio:

cbind(time=unname(sapply(res,function(x) x$time["elapsed"])),
          dev=sapply(res,function(x) deviance(x$fit)))
##                time      dev
## nloptwrap  1001.824 6067.706
## nlminbwrap 3495.671 6068.730
## bobyqa     4945.332 6068.731

Esses desvios estão consideravelmente abaixo do desvio relatado pelo OP do R (6101.7) e ligeiramente abaixo dos relatados pelo OP do SAS (6078.9), embora a comparação dos desvios entre os pacotes nem sempre seja sensata.

Fiquei realmente surpreso que o SAS tenha convergido em apenas cerca de 100 avaliações de funções!

Os tempos variam de 17 minutos ( nloptwrap) a 80 minutos ( bobyqa) em um Macbook Pro, consistente com a experiência do OP. O desvio é um pouco melhor para nloptwrap.

round(cbind(sapply(res,function(x) fixef(x$fit))),3)
##             nloptwrap nlminbwrap bobyqa
## (Intercept)    -5.815     -5.322 -5.322
## ps             -0.989      0.171  0.171
## th             -0.033     -1.342 -1.341
## rs              1.361     -0.140 -0.139
## hw             -2.100     -2.082 -2.082

As respostas parecem bastante diferentes com nloptwrap- embora os erros padrão sejam bastante grandes ...

round(coef(summary(res[[1]]$fit)),3)
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -5.815      0.750  -7.750    0.000
## ps            -0.989      1.275  -0.776    0.438
## th            -0.033      2.482  -0.013    0.989
## rs             1.361      2.799   0.486    0.627
## hw            -2.100      0.490  -4.283    0.000

(o código aqui fornece alguns avisos sobre year:id que devo rastrear)

Continua ... ?

Ben Bolker
fonte
seria mais útil se eu lhe enviasse o conjunto de dados completo? O único problema é que a convergência leva cerca de 9 horas com o conjunto de dados completo; portanto, sua sugestão sobre amostragem é boa. Tentei transformar os dados com uma transformação de log, mas o gráfico residual em caixa ainda é feio - você acha que o gráfico residual explica parte do problema com esses dados? Finalmente - seus resultados no SAS foram semelhantes aos do R?
Nova