Alternativas para relatar um valor-p: o caso de uma tabela de contingência

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon geração de dados nossa, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Freqüentemente me encontro em discussões com colaboradores sobre os méritos de reportar valores-p, particularmente no contexto de estudos-piloto ou análise exploratória. Nos últimos anos, o Associação Americana de Estatística fez várias declarações fortes sobre a necessidade de considerar abordagens que medem a força da evidência ou incerteza que não depende necessariamente de valores-p. Em 2016, a ASA tentou esclarecer o uso e a interpretação adequados do valor-p, destacando os princípios-chave “que poderiam melhorar a conduta ou a interpretação da ciência quantitativa, de acordo com um amplo consenso na comunidade estatística”. Esses princípios são dignos de nota aqui, caso você não chegue ao artigo original:

  • Os valores p podem indicar a incompatibilidade dos dados com um modelo estatístico especificado.
  • Os valores p não medem a probabilidade de que a hipótese estudada seja verdadeira ou a probabilidade de que os dados foram produzidos apenas por acaso.
  • as conclusões científicas e as decisões de negócios ou políticas não devem se basear apenas no valor de p passar um limite específico.
  • inferência adequada requer relatórios completos e transparência
  • um valor-p, ou significância estatística, não mede o tamanho de um efeito ou a importância de um resultado.
  • por si só, um valor-p não fornece uma boa medida de evidência a respeito de um modelo ou hipótese.

Mais recentemente, a ASA elaborou esse ponto de vista, respondendo àqueles que pensavam que o artigo inicial era muito negativo, uma lista de muitas coisas não façam. Neste novo artigo, a ASA argumenta que “saber o que não fazer com valores-p é realmente necessário, mas não é suficiente”. Também precisamos saber o que devemos Faz. Uma dessas coisas deve ser o tamanho dos efeitos (e alguma medida de incerteza, como um intervalo de confiança ou credibilidade), a fim de avaliar uma intervenção ou exposição.

Aplicando o pensamento de princípios a um pequeno problema

Recentemente, eu estava discutindo a apresentação de resultados para um estudo piloto. Eu estava argumentando que deveríamos transmitir as descobertas de uma maneira que destacasse as tendências gerais sem levar os leitores a tirar conclusões excessivamente fortes, o que os valores-p podem fazer. Então, eu estava argumentando que, em vez de apresentar valores-p, deveríamos exibir tamanhos de efeito e intervalos de confiança e evitar usar o conceito de “significância estatística”.

Geralmente, isso não é um problema; podemos estimar um tamanho de efeito como uma diferença de médias, uma diferença de proporções, uma razão de proporções, uma razão de chances ou mesmo o logaritmo de uma razão de chances. Nesse caso, o resultado foi uma pesquisa do tipo Likert, em que a resposta foi “nenhuma”, “um pouco” e “muito”, e havia três grupos de comparação, então tivemos uma (3 times3 ) tabela de contingência com um fator ordinal (ou seja, ordenado). Nesse caso, não está tão claro qual deve ser a medida do tamanho do efeito.

Uma opção é calcular uma ( chi ^ 2 ) estatística, relate o valor p associado e chame-o de dia. No entanto, desde o ( chi ^ 2 ) não é uma medida de efeito e o valor-p não é necessariamente uma boa medida de evidência, considerei estimar um modelo de probabilidades cumulativas que forneceria uma medida da associação entre grupo e resposta. No entanto, fiquei um pouco preocupado, porque a versão típica desse modelo assume uma probabilidade proporcional, que eu não tinha certeza se seria apropriado aqui. (Já escrevi sobre esses modelos antes, aqui e aqui, se você quiser dar uma olhada.) É possível ajustar um modelo de probabilidades cumulativas sem a premissa de proporcionalidade, mas as estimativas são mais difíceis de interpretar, pois o tamanho do efeito varia por grupo e resposta.

Leia Também  psiconetria 0.7, pré-impressão com meta-análise e curso SEM online

Felizmente, existe uma medida mais geral de associação para tabelas de contingência com pelo menos um, mas possivelmente dois, fatores nominais: V do Cramer. Esta medida que não faz suposições sobre proporcionalidade.

Meu plano é simular dados da tabela de contingência e, neste post, explorarei os modelos de probabilidades acumuladas. Da próxima vez, descreverei o V do Cramer medida de associação.

Probabilidades cumulativas não proporcionais

No modelo de chances cumulativas (novamente, dê uma olhada aqui para um pouco mais de descrição desses modelos), assumimos que todas as razões de chances de log são proporcionais. Na verdade, isso pode não ser uma suposição irracional, mas eu queria começar com um conjunto de dados gerado sem assumir explicitamente a proporcionalidade. Na seguinte definição de dados, a distribuição das respostas da pesquisa (Nenhum, um poucoe muito) nos três grupos (1, 2e 3) são especificados exclusivamente para cada grupo:

library(simstudy)

# define the data

def <- defData(varname = "grp", 
            formula = "0.3; 0.5; 0.2", dist = "categorical")

defc <- defCondition(condition = "fgrp == 1", 
            formula = "0.70; 0.20; 0.10", dist = "categorical")
defc <- defCondition(defc, condition = "fgrp == 2", 
            formula = "0.10; 0.60; 0.30", dist = "categorical")
defc <- defCondition(defc, condition = "fgrp == 3", 
            formula = "0.05; 0.25; 0.70", dist = "categorical")

# generate the data

set.seed(99)

dx <- genData(180, def)
dx <- genFactor(dx, "grp", replace = TRUE)
dx <- addCondition(defc, dx, "rating")
dx <- genFactor(dx, "rating", replace = TRUE, 
         labels = c("none", "a little", "a lot"))

dx[]
##       id fgrp  frating
##   1:   1    2 a little
##   2:   2    3 a little
##   3:   3    3    a lot
##   4:   4    2 a little
##   5:   5    2 a little
##  ---                  
## 176: 176    2    a lot
## 177: 177    1     none
## 178: 178    3    a lot
## 179: 179    2 a little
## 180: 180    2 a little

Um gráfico de distribuição baseado nessas 180 observações indica que as probabilidades não são provavelmente proporcionais; o "tell" é a grande protuberância para aqueles em grupo 2 quem responde um pouco.

library(likert)

items <- dx[, .(frating)]
names(items) <- c(frating = "rating")

likert.data <- likert(items = items, grouping = dx$fgrp)
plot(likert.data, wrap = 100, low.color = "#DAECED", 
  high.color = "#CECD7B")

o ( chi ^ 2 ) teste, não surpreendentemente indica que existem diferenças prováveis ​​nas respostas nos três grupos:

chisq.test(table(dx[, .(fgrp, frating)]))
## 
##  Pearson's Chi-squared test
## 
## data:  table(dx[, .(fgrp, frating)])
## X-squared = 84, df = 4, p-value <2e-16

Mas, como estamos tentando fornecer uma imagem mais rica da associação que será menos suscetível a amostras pequenas, aqui está o ajuste cumulativo (proporcional) do modelo de probabilidades usando o método clm função no ordinal pacote.

library(ordinal)

clmFit.prop <- clm(frating ~ fgrp, data = dx)
summary(clmFit.prop)
## formula: frating ~ fgrp
## data:    dx
## 
##  link  threshold nobs logLik  AIC    niter max.grad cond.H 
##  logit flexible  180  -162.95 333.91 5(0)  4.61e-08 2.8e+01
## 
## Coefficients:
##       Estimate Std. Error z value Pr(>|z|)    
## fgrp2    2.456      0.410    5.98  2.2e-09 ***
## fgrp3    3.024      0.483    6.26  3.9e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##                Estimate Std. Error z value
## none|a little     0.335      0.305    1.10
## a little|a lot    2.945      0.395    7.45

Um gráfico das proporções observadas (mostradas pela linha) com as proporções modeladas (mostradas como pontos) indica que o modelo que faz a suposição proporcional pode não estar fazendo um ótimo trabalho:

Leia Também  Lançada a nova empresa consultiva Avonhurst
cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

Se ajustarmos um modelo que não faz a suposição de proporcionalidade e comparamos usando a estatística da AIC (menor é melhor) ou um teste de razão de verossimilhança (pequeno valor p indica que o modelo saturado / não proporcional é melhor), fica claro que o modelo de chances não proporcionais para esse conjunto de dados é mais adequado.

clmFit.sat <- clm(frating ~ 1, nominal = ~ fgrp, data = dx)
summary(clmFit.sat)
## formula: frating ~ 1
## nominal: ~fgrp
## data:    dx
## 
##  link  threshold nobs logLik  AIC    niter max.grad cond.H 
##  logit flexible  180  -149.54 311.08 7(0)  8.84e-11 4.7e+01
## 
## Threshold coefficients:
##                            Estimate Std. Error z value
## none|a little.(Intercept)     0.544      0.296    1.83
## a little|a lot.(Intercept)    1.634      0.387    4.23
## none|a little.fgrp2          -4.293      0.774   -5.54
## a little|a lot.fgrp2         -0.889      0.450   -1.98
## none|a little.fgrp3          -2.598      0.560   -4.64
## a little|a lot.fgrp3         -1.816      0.491   -3.70
anova(clmFit.prop, clmFit.sat)
## Likelihood ratio tests of cumulative link models:
##  
##             formula:       nominal: link: threshold:
## clmFit.prop frating ~ fgrp ~1       logit flexible  
## clmFit.sat  frating ~ 1    ~fgrp    logit flexible  
## 
##             no.par AIC logLik LR.stat df Pr(>Chisq)    
## clmFit.prop      4 334   -163                          
## clmFit.sat       6 311   -150    26.8  2    1.5e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

É possível que o ajuste inadequado seja apenas uma ocorrência rara. Abaixo está um gráfico que mostra o resultado médio ( ( pm 1 sd )) para o modelo 1000 se encaixa em 1000 conjuntos de dados usando o mesmo processo de geração de dados. Parece que os resultados iniciais não foram uma aberração - o modelo de chances proporcionais se encaixa em uma estimativa tendenciosa, particularmente para grupos 1 e 2. (O código para fazer esta simulação é mostrado no adendo.)

Premissa proporcional cumprida

Aqui, o processo de geração de dados é modificado para que a suposição de proporcionalidade seja incorporada.

def <- defData(varname = "grp", formula = ".3;.5;.2", 
               dist = "categorical")
def <- defData(def, varname = "z", formula = "1*I(grp==2) + 2*I(grp==3)", 
               dist = "nonrandom")

baseprobs <- c(0.7, 0.2, 0.1)

dx <- genData(180, def)
dx <- genFactor(dx, "grp", replace = TRUE)
dx <- genOrdCat(dx, adjVar = "z", baseprobs, catVar = "rating")
dx <- genFactor(dx, "rating", replace = TRUE,
          labels = c("none", "a little", "a lot")
)

É assim que as probabilidades proporcionais - não existem protuberâncias óbvias, apenas uma mudança geral para a direita à medida que saímos do grupo 1 para 3:

Quando ajustamos o modelo proporcional e o comparamos ao modelo saturado, não vemos razão para rejeitar a suposição de proporcionalidade (com base nas estatísticas da AIC ou da LR).

clmFit.prop <- clm(frating ~ fgrp, data = dx)
summary(clmFit.prop)
## formula: frating ~ fgrp
## data:    dx
## 
##  link  threshold nobs logLik  AIC    niter max.grad cond.H 
##  logit flexible  180  -176.89 361.77 4(0)  3.04e-09 2.7e+01
## 
## Coefficients:
##       Estimate Std. Error z value Pr(>|z|)    
## fgrp2    1.329      0.359    3.70  0.00022 ***
## fgrp3    2.619      0.457    5.73    1e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##                Estimate Std. Error z value
## none|a little     0.766      0.299    2.56
## a little|a lot    2.346      0.342    6.86
clmFit.sat <- clm(frating ~ 1, nominal = ~ fgrp, data = dx)
anova(clmFit.prop, clmFit.sat)
## Likelihood ratio tests of cumulative link models:
##  
##             formula:       nominal: link: threshold:
## clmFit.prop frating ~ fgrp ~1       logit flexible  
## clmFit.sat  frating ~ 1    ~fgrp    logit flexible  
## 
##             no.par AIC logLik LR.stat df Pr(>Chisq)
## clmFit.prop      4 362   -177                      
## clmFit.sat       6 365   -177    0.56  2       0.75

E aqui está um gráfico que resume um segundo conjunto de 1000 iterações, este usando a suposição de probabilidades proporcionais. As estimativas parecem não imparciais:

Leia Também  Como adicionar uma vinheta a um pacote no RStudio

Suspeito que, em muitos casos, as respostas do tipo Likert se pareçam mais com o segundo caso que com o primeiro, de modo que o modelo de chances proporcionais cumulativas possa muito bem ser útil para caracterizar a associação entre grupo e resposta. Mesmo que a suposição não seja razoável, o viés pode não ser terrível e a estimativa ainda pode ser útil como uma medida de associação. No entanto, podemos preferir uma medida livre de quaisquer suposições, como V do Cramer. Eu vou falar sobre isso na próxima vez.

Referências:

Ronald L. Wasserstein e Nicole A. Lazar (2016) A Declaração da ASA sobre p-Valores: Contexto, Processo e Propósito, The American Statistician, 70: 2, 129-133.

Ronald L. Wasserstein, Allen L. Schirm e Nicole A. Lazar (2019) Mudando para um mundo além de "p <0,05", The American Statistician, 73: sup1, 1-19.

Adendo: código para análise replicada

library(parallel)
RNGkind("L'Ecuyer-CMRG")  # to set seed for parallel process

dat.nonprop <- function(iter, n) {
  
  dx <- genData(n, def)
  dx <- genFactor(dx, "grp", replace = TRUE)
  dx <- addCondition(defc, dx, "rating")
  dx <- genFactor(dx, "rating", replace = TRUE,
            labels = c("none", "a little", "a lot")
  )

  clmFit <- clm(frating ~ fgrp, data = dx)
  
  dprob.obs <- data.table(iter, 
      prop.table(dx[, table(fgrp, frating)], margin = 1))
  
  setkey(dprob.obs, fgrp, frating)
  setnames(dprob.obs, "N", "p.obs")
  
  dprob.mod <- data.table(iter, fgrp = levels(dx$fgrp),
      predict(clmFit, newdata = data.frame(fgrp = levels(dx$fgrp)))$fit)
  
  dprob.mod <- melt(dprob.mod, id.vars = c("iter", "fgrp"), 
                    variable.name = "frating", value.name = "N")
  
  setkey(dprob.mod, fgrp, frating)
  setnames(dprob.mod, "N", "p.fit")
  
  dprob <- dprob.mod[dprob.obs]
  dprob[, frating := factor(frating, 
                        levels=c("none", "a little", "a lot"))]
  
  dprob[]
  
}

def <- defData(varname = "grp", formula = ".3;.5;.2", 
            dist = "categorical")

defc <- defCondition(condition = "fgrp == 1", 
            formula = "0.7;0.2;0.1", dist = "categorical")
defc <- defCondition(defc, condition = "fgrp == 2", 
            formula = "0.1;0.6;0.3", dist = "categorical")
defc <- defCondition(defc, condition = "fgrp == 3", 
            formula = "0.05;0.25;0.70", dist = "categorical")

res.nonp <- rbindlist(mclapply(1:1000, 
                        function(iter) dat.nonprop(iter,180)))

sum.nonp <- res.nonp[, .(mfit = mean(p.fit), sfit = sd(p.fit), 
              mobs = mean(p.obs), sobs = sd(p.obs)), 
              keyby = .(fgrp, frating)]

sum.nonp[, `:=`(lsd = mfit - sfit, usd = mfit + sfit)]

ggplot(data = sum.nonp, aes(x = frating, y = mobs)) +
  geom_line(aes(group = fgrp), color = "grey60") +
  geom_errorbar(aes(ymin = lsd, ymax = usd,  color = fgrp), 
                width = 0) +
  geom_point(aes(y = mfit, color = fgrp)) +
  theme(panel.grid = element_blank(),
        legend.position = "none",
        axis.title.x = element_blank()) +
  facet_grid(fgrp ~ .) +
  scale_y_continuous(limits = c(0, 0.85), name = "probability") +
  scale_color_manual(values = c("#B62A3D", "#EDCB64", "#B5966D"))

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'};

(função (d, t) {
var s = d.createElement
s.src = '//cdn.viglink.com/api/vglnk.js';
var r = d.getElementsByTagName
} (documento, 'script'));

Para Deixe um comentário para o autor, siga o link e comente no blog: geração de dados nossa.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da dados / R.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br