Erros do Tipo M e do Tipo S

Erros do Tipo M e do Tipo S

Tipos de erro na testagem de hipóteses

Quando aprendemos sobre a testagem de hipótese nula nas aulas de estatística, geralmente somos apresentados a dois tipos de erro, batizados (sem muita imaginação) de erros do Tipo I e erros do Tipo II. Partimos de uma divisão binária do mundo: em uma versão, a hipótese nula H_0 é verdadeira; noutra, é a hipótese alternativa H_A. Quando conduzimos o teste de hipótese, desconhecemos o verdadeiro estado do mundo e somos obrigados a tomar uma decisão, também binária, com base nos resultados observados: ou obtemos dados incompatíveis com H_0 e a rejeitamos; ou obtemos dados razoavelmente compatíveis com a hipótese nula e decidimos por não a rejeitar – o que nem sempre implica que podemos aceitá-la como verdadeira, contudo.

Nessa divisão binária do verdadeiro estado do mundo e das decisões possíveis, podemos cometer dois tipos de erro: rejeitar H_0 quando ela é verdadeira, ou manter H_0 quando, em verdade, ela é falsa. Temos, respectivamente, o erro do Tipo I e do Tipo II. Nada muito complicado, a princípio! Mas nem sempre o professor lembra ao aluno que nunca chegamos a conhecer o verdadeiro estado do mundo e, portanto, não temos como saber com certeza se cometemos um erro, seja do Tipo I ou do Tipo II, ou se tomamos a decisão correta em um estudo específico. O máximo que podemos fazer é controlar a taxa de erro ao longo de uma série de experimentos ou decisões.

Via de regra, o maior receio do pesquisador é cometer um erro do Tipo I: rejeitar a hipótese nula quando não há nada de interessante acontecendo. Isso implica que o pesquisador gostaria de ser conservador na hora de avaliar as evidências disponíveis para sua hipótese de interesse, assumindo que essa hipótese foi propriamente traduzida para a hipótese estatística H_A. Por esse motivo, a taxa de erro do Tipo I aceitável costuma ser fixada pelo pesquisador a valores relativamente baixos, como 5% ou 1%. Damos à taxa de erro do Tipo I o nome \alpha e ela indica com que frequência iremos rejeitar a hipótese nula em um mundo em que ela é verdadeira.

Em segundo plano, e muitas vezes completamente esquecido, temos o erro do Tipo II: deixar de rejeitar a hipótese nula quando ela é falsa. Controlar a taxa de erro do Tipo II não é tão simples. Isso porque definir a hipótese nula não costuma ser um grande problema: em geral, ela implica apenas que o parâmetro de interesse é igual a zero (a chamada “nil null hypothesis”). A hipótese alternativa, porém, não é tão facilmente definida. Nas ciências sociais, em geral, as teorias são pobres e não produzem uma estimativa pontual ou mesmo um intervalo dos valores aceitáveis para os parâmetros de interesse. Para calcular a taxa de erro do Tipo II, batizado de \beta, o pesquisador precisa definir, minimamente, qual o limite inferior do parâmetro de interesse sob H_A e qual é a variância esperada das estimativas. Afinal de contas, há diversas maneiras como o verdadeiro estado do mundo pode diferir de H_0, e para cada uma teremos uma taxa de erro Tipo II diferente.

No lugar de discutir e computar o erro do Tipo II, é mais comum que o problema da detecção de desvios importantes da hipótese nula seja definido em termos da potência estatística de um teste. A potência é definida como 1-\beta, ou seja, a probabilidade de não cometer um erro do Tipo II, ou, ainda, de rejeitar a hipótese nula quando ela é, de fato, falsa. Os pesquisadores não costumam ser tão obsessivos com relação à potência mínima para seus estudos – muitas vezes esquecendo dela por completo – mas é comum fixar a potência em 80%.

Outra dificuldade em controlar o erro do Tipo II se deve ao fato de ele depender também das taxas do erro do Tipo I fixadas como aceitáveis. Exigir uma taxa de erro do Tipo I muito baixa implica, necessariamente, em um aumento na taxa de erro do Tipo II. Reduzir os dois tipos de erro simultaneamente só é possível se aumentar o tamanho da amostra ou melhorar a qualidade das medidas, de forma a reduzir o erro-padrão das estimativas de interesse.

O problema para a psicologia (e outras ciências sociais)

O framework da testagem da hipótese nula e o uso quase irrestrito dos erros de Tipo I e II para análise de dados na psicologia vai contra as concepções mais ingênuas dos determinantes dos fenômenos psicológicos. Por mais simplista que uma teoria psicológica possa ser, dificilmente um pesquisador da área aceitaria a ideia de que a relação entre os fenômenos poderia ser reduzida a uma equação simples e elegante, como muitas vezes é possível nas ciências naturais. Pelo contrário, qualquer conclusão mais substancial e restritiva sobre um fenômeno psicológico é facilmente criticada por argumentos relativos aos efeitos do contexto, à variabilidade entre sujeitos, à diversidade de operacionalizações possíveis.

Paul Meehl chama essa dependência irrestrita entre fenômenos de crud factor. Em última instância, a sobredeterminação dos fenômenos psicológicos implica que “tudo se correlaciona com tudo”, mesmo que essa correlação seja fraca demais para ser estimada com uma amostra pequena. A ausência de zeros verdadeiros nos parâmetros de qualquer modelo psicológico torna a utilização da testagem de hipótese nula um exercício de futilidade. Com uma amostra grande o suficiente, é possível rejeitar qualquer hipótese nula porque ela sempre será falsa.

O mesmo pode ser esperado de tratamentos psicológicos. Dificilmente o efeito será exatamente igual a zero, mesmo que seja muito pequeno para ser clinicamente significante. E um efeito fraco dificilmente será estável: sob determinadas configurações, o efeito médio pode ser positivo; sobre outros, negativo; e a variação desse efeito dificultará formular qualquer enunciado teórico preciso sobre o fenômeno.

E isso não é tudo: as taxas de erro não dependem exclusivamente da definição da hipótese nula ou da hipótese alternativa. Todo modelo estatístico é baseado em pressupostos sobre o fenômeno e o experimento que, na maioria das vezes, só são aceitáveis aproximadamente. Violações desses pressupostos afetam as estatísticas utilizadas para testagem de hipóteses nulas. E isso não modifica apenas estatísticas de teste: toda medida de incerteza que obtemos a partir de um modelo estatístico é condicionada à certeza de que o modelo gerador dos dados é, definitivamente, aquele utilizado para calculá-la. Em alguns raros casos temos uma leve segurança quanto ao modelo mais adequado; na maioria das vezes, porém, o modelo é puramente heurístico e poderia facilmente ser substituído por outro. Nesses casos, valores-p, intervalos de confiança e erros-padrão são apenas um limite inferior da incerteza sobre o parâmetro.

Em função dos problemas relativos aos erros do Tipo I e II, uma proposta razoavelmente recente sugere deixar o paradigma da testagem de hipótese nula um pouco de lado e substitui-la (ou, pelo menos, complementá-la) pelo uso de estimadores e suas propriedades para resumir os resultados de uma pesquisa. Parece simples, uma vez que não precisamos abandonar o maquinário estatístico necessário para calcular valores-p em testes de hipótese, mas a adoção de uma perspectiva focada na estimação em vez na testagem implica deixar um pouco de lado a ideia binária simplista implicada pelo erros de Tipo I e II. Trocamos uma posição epistemológica que efetua um corte limpo e distinto entre o verdadeiro e falso por uma que permite pensar o processo de investigação de maneira contínua, quantificando o erro em termos da distância esperada entre os valores obtidos em um estudo e o valor verdadeiro hipotético.

Uma alternativa: erros do Tipo S e Tipo M

Mesmo abandonando a visão binária dos resultados de um estudo, é crucial termos alguma noção dos erros incorridos pelas inferências. Mas como tratar os erros sob uma perspectiva contínua? Uma sugestão feita por Andrew Gelman e John Carlin é categorizar os erros inferenciais em duas classes: erros do Tipo M e do Tipo S.

Erros do Tipo M se referem à diferença de magnitude entre a estimativa obtida e o parâmetro de interesse. Trabalhar com uma perspectiva de estimação implica que o erro do Tipo M sempre estará presente, em maior ou menor medida. Mas esse erro tende a ser maior quando a potência estatística do estudo é baixa e as conclusões são formuladas a partir da significância estatística de um teste de hipótese. Nesses casos, uma estimativa estatisticamente significativa estará necessariamente inflando a magnitude do efeito, e é possível calcular, com base no erro do Tipo M, qual será a razão de exagero presente em um estudo. O valor da razão de exagero é dado simplesmente pela razão do valor esperado do estimador condicionado à ocorrência de significância estatística e valor verdadeiro:

\displaystyle E_F = \frac{\mathrm{E}[\hat{\theta} \mid p < \alpha]}{\theta}

Os erros do Tipo S, por sua vez, ocorrem quando a estimativa de um estudo apresenta um sinal na direção contrária do parâmetro verdadeiro. Um estudo com potência estatística suficiente dificilmente incorrerá nesse tipo de erro. Mas, como a maioria das magnitudes de efeito reportadas na literatura estão infladas em função da seleção baseada na significância estatística, a potência de um estudo é quase sempre bem menor do que a reportada. Com uma potência menor que a esperada, a probabilidade de uma estimativa significativa estar na direção errada aumenta, tendo em vista que ela é definida pela razão das áreas sob as caudas da distribuição amostral:

\displaystyle \mathrm{Pr}(\mathrm{Erro Tipo S}) = \mathrm{Pr}(\hat{\theta} < 0 | p < \alpha)

Um exemplo: estimando a reprodutibilidade da ciência psicológica

Em 2015, uma grande equipe de pesquisadores da área da psicologia publicou um artigo que gerou um debate acirrado sobre a reprodutibilidade dos resultados encontrados nas pesquisas da área. Ao reproduzirem 100 estudos amostrados de periódicos de referência, apenas 36% apresentaram um resultado estatisticamente significativo. Como as replicações foram desenhadas de forma mais rigorasa que os estudos originais – diminuindo os graus de liberdade na análise e procurando manter uma potência estatística mínima de 90% – essa taxa de replicação certamente aponta para problemas sobre as pesquisas na psicologia. A discussão não é simples, porém, e os próprios autores debatem sobre o melhor critério para concluir se o resultado de um estudo foi replicado ou não.

Os autores relatam também a diferença na distribuição das estimativa das magnitudes de efeito. Porém, eles permanecem fiéis ao paradigma dos erros de Tipo I e II para refletir sobre os resultados. Como contraponto e exemplo, podemos aplicar a proposta de avaliar os erros do Tipo M e S dos resultados obtidos pelas replicações.

Para facilitar a comparação entre os diferentes estudos, os autores padronizaram todas as estimativas da magnitude de efeito como coeficientes de correlação de Pearson \rho. A transformação não é das mais apropriadas para facilitar nossos cálculos, já que o coeficiente de correlação tem seu domínio limitado a [-1; 1] e sua distribuição, além de assimétrica, depende da distribuição das variáveis utilizadas para computá-lo. Para simplificar os cálculos e permitir a utilização dos resultados baseado em uma aproximação normal, vamos considerar a transformação de Fisher ao coeficiente \rho: F(\rho) = \mathrm{tanh^{-1}}(\rho). Aplicando o inverso da função tangente hiperbólica ao coeficiente de correlação, conseguimos mapeá-lo para o domínio [-\infty; +\infty] e sua distribuição amostral se torna aproximadamente normal, com média \mathrm{tanh^{-1}}(\rho) e desvio-padrão \frac{1}{\sqrt{N-3}}.

A distribuição das magnitudes de efeito padronizadas possui forte assimetria e, por isso, vamos trabalhar com sua mediana. A mediana da magnitude de efeito dos estudos originais foi de cerca \rho = 0.4. Para obter uma potência de 80% para detectar um efeito dessa magnitude, considerando um \alpha = 0.05, é necessária uma amostra de pelo menos 47 observações. Se essa estimativa estiver correta e o pesquisador conseguir obter suas 47 observações, a probabilidade de cometer um erro do Tipo S é ínfima (algo na ordem de 10^{-6}) e a razão de exagero também é baixa, cerca de 9,6%.

r_0 <- 0.0
r_A <- 0.4
N <- 47
alpha <- 0.05
# Calcula o valor crítico para decidir pela significância estatística
# com um alpha = 0.05, considerando que a distribuição nula é normal
# com média = atanh(r_0) e desvio-padrão = 1/sqrt(N-3)
criticalT <- qnorm(alpha/2, atanh(r_0), 1/sqrt(N-3))

# Calcula a probabilidade de erro do Tipo S
pnorm(criticalT, atanh(r_A), 1/sqrt(N-3)) / 
  (pnorm(criticalT, atanh(r_A), 1/sqrt(N-3)) + 
     pnorm(-criticalT, atanh(r_A), 1/sqrt(N-3), lower.tail = F))

# Calcula a razão de exagero usando Monte Carlo

sims  abs(criticalT)])) / r_A

A distribuição das magnitudes de efeito das replicações, todavia, mostra um quadro bem diferente. Como esperado, a maioria das estimativas dos estudos publicados é bem maior do que o valor observado na replicação. A mediana da distribuição dos efeitos nas replicações é de cerca de \rho=0.12 – praticamente um quarto da mediana dos estudos originais! Como as replicações foram desenhadas cuidadosamente para obter uma potência igual ou superior a 80% e as análises foram pré-registradas, diminuindo substancialmente as práticas questionáveis de pesquisa para obter resultados publicáveis, as estimativas das replicações são mais confiáveis.

Mas isso tem consequências graves para os erros do Tipo M e do tipo S. Se o pesquisador desenha seu estudo para ter 80% de potência para uma magnitude de efeito \rho=0.4, mas a magnitude de efeito é apenas de \rho=0.12, sua potência cai para 12%! Isso, por si só, já deveria servir para tomarmos cautela sobre a conclusão de qualquer estudo, mesmo quando seu desenho foi projetado para obtenção de uma alta potência estatística baseada nas estimativas publicadas na literatura.

Sob esse regime de baixa potência estatística, o erro de Tipo S salta para 2,3%. Ou seja: cerca de dois a cada cem estudos que obtém significância estatística apresentam estimativas com sinal contrário ao do efeito verdadeiro. O erro do Tipo M também sofre uma inflação considerável: a razão de exagero é de 2.9. Isso significa que, se fizermos a média de vários estudos publicados (ou seja, que obtiveram significância estatística) que tentam estimar um parâmetro com 12% de potência, o valor será três vezes maior do que o efeito verdadeiro. Coincidência ou não, a mediana dos efeitos publicados é cerca de quatro vezes maior do que a mediana dos efeitos replicados.

r_A <- 0.12

# Calcula a probabilidade de erro do Tipo S
pnorm(criticalT, atanh(r_A), 1/sqrt(N-3)) / 
  (pnorm(criticalT, atanh(r_A), 1/sqrt(N-3)) + 
     pnorm(-criticalT, atanh(r_A), 1/sqrt(N-3), lower.tail = F))

# Calcula a razão de exagero usando Monte Carlo

sims  abs(criticalT)])) / r_A

A imagem abaixo apresenta os erros de Tipo M e Tipo S de maneira mais intuitiva. A distribuição amostral do parâmetro de interesse sob a hipótese nula é representada pela linha pontilhada. As linhas tracejadas curtas nas caudas dessa distribuição indica o valor crítico bicaudal para declarar significância estatística com \alpha = 0.05. A linha cheia indica a distribuição amostral do parâmetro assumindo que o valor verdadeiro da magnitude de efeito é \rho = 0.12. Colorido em azul e vermelho, as áreas sob as caudas indicam o evento estimativa estatisticamente significativa e sua soma é a potência estatística (12%, no caso).

A razão entre a a área em vermelho e a soma das duas áreas coloridas dá a probabilidade da ocorrência de um erro do Tipo S. O erro do Tipo M é resumido pelo razão de exagero, calculada como o valor esperado da estimativa para as áreas coloridas. Como a região em azul é a que mais contribui para o valor esperado, ela é denotada como Erro do Tipo M no gráfico.

library(ggplot2)

plotZ <- data.frame(Z=seq(-0.41, 0.65, length.out=100),
                    Densidade=dnorm(seq(-0.41, 0.65, length.out=100), atanh(r_A), 1/sqrt(N-3)),
                    Densidade_Null=dnorm(seq(-0.41, 0.65, length.out=100), atanh(0), 1/sqrt(N-3)))
shadeZ1 <- subset(plotZ, Z <= criticalT)
shadeZ2 = -criticalT)
ggplot(plotZ, aes(x=Z, y=Densidade)) + geom_line() + geom_line(aes(x=Z, y=Densidade_Null), linetype=3) +
  ylab('Densidade') + 
  ggtitle(expression(paste('Distribuição de ', tanh^{-1}, (rho), ' sob ', H[A]))) +
  geom_segment(aes(x=criticalT, xend=criticalT, y=0, yend=dnorm(criticalT, 0, 1/sqrt(N-3))), linetype=2) +
  geom_segment(aes(x=-criticalT, xend=-criticalT, y=0, yend=dnorm(-criticalT, 0, 1/sqrt(N-3))), linetype=2) +
  geom_ribbon(data=shadeZ1, aes(x=Z, ymin=0, ymax=Densidade), alpha=0.5, fill='indianred')+
  geom_ribbon(data=shadeZ2, aes(x=Z, ymin=0, ymax=Densidade), alpha=0.5, fill='skyblue')+
  geom_vline(xintercept=atanh(r_A), linetype=1)+
  geom_vline(xintercept=atanh(0), linetype=3) + 
  xlab(expression(paste(tanh^{-1}, '(', rho, ')')))+
  annotate('text', 0.5, 1, label='Erro Tipo M') +
  geom_segment(aes(x=0.5, y=0.9, xend=0.37, yend=0.5), arrow=arrow(length=unit(3, 'mm'))) +
  annotate('text', -0.3, 1, label='Erro Tipo S') +
  geom_segment(aes(x=-0.3, y=0.9, xend=-0.35, yend=0.1), arrow=arrow(length = unit(3, 'mm')))

plot of chunk unnamed-chunk-81

E daí?

Muito da discussão sobre os resultados do projeto de reprodutibilidade da psicologia recaiu sobre a pergunta: com uma taxa estimada de 36% de reprodutibilidade, estaria a psicologia perseguindo ruído mais do que produzindo conhecimento? Essa discussão acaba se baseando no framework dos Erros do Tipo I e II, assim como o próprio critério de reprodutibilidade. Por outro lado, avaliando os resultados em termos dos Erros do Tipo M e S, podemos lançar algumas novas luzes sobre a discussão e sobre a utilização da estatística na produção de conhecimento na psicologia.

Os efeitos de interesse na psicologia dificilmente serão exatamente igual a zero. Sob essa condição, a utilização de testagem nula tem utilidade limitada – com potência estatística adequada, qualquer hipótese nula pode ser rejeitada em favor da hipótese alternativa. Um estudo piloto até pode se interessar em fazer uma avaliação preliminar da significância prática de um efeito pouco estudado, mas, de maneira geral, interessa mais poder estimar a magnitude do efeito com algum grau de precisão.

Com a tendência dos valores publicados estarem superestimando os parâmetros de interesse em função da seleção pela significância estatística, estariam as pesquisas conseguindo estimar os efeitos com fidedignidade? Eles possuem magnitude suficiente para permanecerem estáveis sob diversas condições, ou são sensíveis a mudanças nos valores de variáveis com as quais não estariam, em teoria, relacionados? A magnitude dos efeitos e as condições que determinam sua variação são estáveis o suficiente para poderem contribuir para um teoria propriamente científica – ao invés de um punhado de pequeno fatos empíricos sem relação clara entre si?

Sob o paradigma dos erros do Tipo M e S, os resultados encontrados pelo projeto de reprodutibilidade da psicologia podem ser reconsiderados. Em vez da polarização entre o (copo meio vazio) de que “mais da metade dos resultados são falsos ou não foram replicados” e o (copo meio cheio) de que “quase metade dos resultados das pesquisas são verdadeiros ou foram replicados”, podemos concluir que os pesquisadores da área da psicologia estão sistematicamente superestimando os efeitos de interesse por causa da potência estatística medíocre; por estarem decidindo sobre a existência ou não dos efeitos com base exclusivamente na significância estatística; e por estarem manipulando os valores-p para obter significância estatística de maneira óbvia e intencional ou sutilmente e com base em boas intenções.

Anúncios

Distribuição do valor-p

Distribuição do valor-p

A grande maioria dos estudos quantitativos em psicologia faz uso do teste de hipótese nula (NHST, null hypothesis significance testing) como principal ferramenta inferencial. Seu uso é tão amplo que é costume reportar o resultado de testes de hipótese mesmo quando eles não são aplicáveis, como no caso de estudos exploratórios. Não vou entrar aqui no problema que é o uso desregrado de testes de hipótese nula em qualquer tipo de estudo, até porque a literatura sobre as limitações e mau uso do NHST é bastante ampla, ainda que nem sempre faça parte da formação dos pesquisadores.

Vamos nos focar aqui em uma outra questão: qual é a distribuição do valor-p? Um ponto que muitas vezes passa despercebido nas aulas sobre NHST é que o valor-p é uma função dos dados da amostra e, por isso, é uma estatística que possui uma distribuição amostral específica. Isso significa que o valor-p irá necessariamente variar nas diversas repetições de um experimento, mesmo quando o efeito de interesse permanece constante. Como veremos, compreender como o valor-p se distribui em diferentes situações nos permitirá desenvolver algumas intuições sobre as justificativas de seu uso para controlar a taxa de erro em estudos que fazem bom uso do teste estatístico de hipótese nula.

O que é um valor-p?

Primeiro, uma revisão básica da definição do conceito. O valor-p é geralmente definido como a probabilidade de se obter uma estatística de teste igual ou maior que a observada na amostra. Uma estatística de teste é simplesmente uma função dos dados da amostra necessariamente sensível à hipótese que se deseja testar e cuja distribuição é conhecida sob a hipótese nula.

Vamos exemplificar utilizando uma das situações mais comuns na qual é aplicado o NHST: a comparação da média entre dois grupos. Por algum motivo teórico, prevemos que duas populações possuem média diferente com relação a alguma variável de interesse. Esse ponto de partida, apesar de ser muitas vezes varrido para debaixo do tapete, é crucial: de forma alguma o teste de hipótese nula ajudará a corroborar nossa teoria científica se a diferença predita for derivável de outras teorias incompatíveis com a que estamos pondo em risco. Em particular, é interessante que a própria hipótese nula possa ser motivada por alguma teoria substancial, e não simplesmente a negação de nossa teoria de bolso. Como construir boas hipóteses em psicologia, porém, é assunto para outro post.

Vamos supor, então que nossa teoria prevê uma diferença na média de uma variável entre dois grupos, enquanto outra teoria corrente prevê que esses grupos são necessariamente iguais. A hipótese nula, motivada pela teoria concorrente com a de nosso interesse, prediz que H_0 := \mu_A = \mu_B. Nossa hipótese, por outro lado, assume que as médias são diferentes, sem se preocupar com a direção da diferença: H_A := \mu_A \neq \mu_B. Conseguimos obter uma amostra aleatória simples dos grupos de interesse – outro problema que também é costumeiramente ignorado na utilização da NHST – e queremos saber se os dados são compatíveis com a hipótese nula ou com a hipótese alternativa prevista por nossa teoria.

Construindo a estatística de teste

Qual seria a melhor estatística de teste para avaliar qual a hipótese corroborada pelos dados? Podemos expressar a hipótese nula em termos da diferença das médias populacionais: H_0 := \mu_A - \mu_B = 0. Uma boa estatística de teste deve refletir a magnitude dessa diferença; uma solução é encontrar uma forma de estimar as médias populacionais e sua diferença. Seja apelando para a estimação por máxima verossimilhança, mínimos quadrados ou a lei dos grandes números, a média da amostra é o melhor estimator não-enviesado da média populacional. Portanto, nosso estimador para a diferença das médias é igual a \Delta\bar{X} := \bar{X}_A - \bar{X}_B.

O problema, porém, é que esse estimador sofre de erro amostral e, por isso, poderá ser diferente do valor postulado pela hipótese nula mesmo no caso de ela ser verdadeira. Isso significa que, a cada vez que repetimos o experimento de comparar os grupos de interesse, obtemos resultados diferente! Precisamos, portanto, ter uma ideia da magnitude dessa variação para averiguar se o valor observado é compatível com os valores preditos. Por isso, vamos verificar algumas propriedades a respeito dos momentos do nosso estimador:

\displaystyle \mathbb{E}[\Delta\bar{X}] = \mathbb{E}[\bar{X}_A-\bar{X}_B] = \mu_A- \mu_B
\displaystyle \mathrm{var}[\Delta\bar{X}] = \mathrm{var}[\bar{X}_A-\bar{X}_B] = \mathrm{var}[\bar{X}_A]+\mathrm{var}[\bar{X}_B] = \frac{\sigma^2_A}{n_A} + \frac{\sigma^2_B}{n_B}

Como a média amostral é um estimador não-enviesado da média populacional, o valor esperado do nosso estimador é igual a diferença entre a média dos dois grupos. A variância do estimador, por sua vez, é igual à soma da variância dos estimadores das médias de cada grupo e, para amostras independentes e identicamente distribuídas, a variância da média amostral é dada pela razão entre a variância populacional e o tamanho da amostra.

Para criarmos uma estatística de teste facilmente interpretável mesmo quando as variáveis em questão mudam de média e variância, vamos padronizar nosso estimador para que ele possua valor esperado igual a zero e variância igual a um. Para isso, basta subtrairmos a diferença populacional do valor estimado e dividir o resultado pela raiz quadrada da variância do estimador – seu erro-padrão. Vamos batizar essa variável de “Z”.

\displaystyle Z = \frac{\Delta\bar{X} - (\mu_A - \mu_B)}{\sqrt{\frac{\sigma^2_A}{n_A} + \frac{\sigma^2_B}{n_B}}}

Mas chegamos a um impasse. Para calcular o valor de Z, precisamos saber a média e a variância populacional! Para facilitar nossas contas e tornar o exemplo mais claro, vamos supor que, de fato, sabemos a variância populacional da variável aleatória de interesse, mas não conhecemos as médias e, portanto, não podemos calcular sua diferença. É nesse ponto que a definição da hipótese nula se torna útil e a estatística de teste revela seu poder em identificar discrepâncias entre hipóteses estatísticas e dados observados. Se a hipótese nula é verdadeira, a diferença das médias é igual a zero. Não sabemos se é o caso, mas vamos fazer de conta que sim. Nossa estatística Z se torna uma estatística de teste definida sob H_0. Se a hipótese for verdadeira, o valor de Z deve estar próximo de zero; se ela for falsa, o valor distoará de zero.

Quantificando a discrepância observada

Porém: quão grande a diferença precisa ser para indicar uma discrepância significativa com relação ao valor esperado sob a hipótese nula? Poderíamos escolher um valor arbitrário para essa decisão, mas vamos explorar o problema em termos probabilísticos. Qual é a probabilidade de observamos um determinado valor de Z caso a hipótese nula seja verdadeira? Se a variável aleatória de interesse é contínua, Z também será contínua e, por isso, a probabilidade de observar um valor exato específico será sempre igual a zero. Uma saída, então, é computar a área sob a cauda da função de densidade de probabilidade da estatística de teste. Para isso, todavia, precisamos saber a forma funcional dessa distribuição! Na verdade, não precisamos: se estamos dispostos a aceitar um limite superior para essa probabilidade, podemos utilizar a desigualdade de Chebyshev e esquecermos o problema da forma funcional da distribuição de Z:

\displaystyle \mathrm{Pr}(|X - \mu_X| \geq |x|) \leq \frac{\sigma^2_X}{x^2}
\displaystyle \mathrm{valor-p} := \mathrm{Pr}(|Z| \geq |z|) \leq \frac{1}{z^2}

Utilizando a desigualdade de Chebyshev, definimos o valor-p como o limite superior da probabilidade de observar um valor igual ou maior (em magnitude) do que aquele computado a partir da amostra.

Mas, como estamos trabalhando com a diferença entre médias, podemos obter uma definição ainda melhor. Sob duas circunstâncias nós conhecemos a distribuição amostral de Z: se X_A e X_B possuem distribuição normal, então, independente do tamanho da amostra, suas médias amostrais também terão distribuição normal; ou, independente da forma funcional da distribuição das variáveis aleatórias, se temos uma amostra grande o suficiente, o Teorema do Limite Central também nos garante que Z possui distribuição normal. Portanto, denotando a função de probabilidade cumulativa da normal padronizada como \Phi, podemos definir um valor exato para o valor-p:

\displaystyle Z \sim \mathrm{Normal}(0, 1; H_0)

library(ggplot2)
plotZ <- data.frame(Z=seq(-3.5, 3.5, length.out=100),
                    Densidade=dnorm(seq(-3.5, 3.5, length.out=100)))
ggplot(plotZ, aes(x=Z, y=Densidade)) + geom_line() +
  xlab('Z') + ylab('Densidade') + ggtitle('Distribuição de Z sob H0')

plot of chunk unnamed-chunk-80

\displaystyle \mathrm{valor-p} := \mathrm{Pr}(|Z| \geq |z|) = \mathrm{Pr}(Z \geq |z|) + \mathrm{Pr}(Z \leq -|z|) = 1-\Phi(|z|) + \Phi(-|z|) = 2\Phi(-|z|)

shadeZ1 <- subset(plotZ, Z <= -1.5)
shadeZ2 <- subset(plotZ, Z >= 1.5)
ggplot(plotZ, aes(x=Z, y=Densidade)) + geom_line() +
  xlab('Z') + ylab('Densidade') + ggtitle('Distribuição de Z sob H0') +
  geom_segment(aes(x=-1.5, xend=-1.5, y=0,yend=dnorm(-1.5)), linetype=2) +
  geom_segment(aes(x=1.5, xend=1.5, y=0,yend=dnorm(1.5)), linetype=2) +
  geom_ribbon(data=shadeZ1, aes(x=Z, ymin=0, ymax=Densidade), alpha=0.5, fill='skyblue')+
  geom_ribbon(data=shadeZ2, aes(x=Z, ymin=0, ymax=Densidade), alpha=0.5, fill='skyblue') +
  geom_text(aes(x=1, y=0.05, label='z = 1.5')) +
  geom_text(aes(x=0, y=.1, label='p = 0.13'))

plot of chunk unnamed-chunk-81

O valor-p é, então, uma medida de discrepância do valor observado para a estatística de teste na amostra obtida com relação aos valores esperados sob a hipótese nula. Essa discrepância é conceituada como a probabilidade de se observar um valor igual ou mais extremos na distribuição da estatística de teste sob a hipótese nula.

Distribuição do valor-p sob a hipótese nula

Vamos utilizar a última definição do valor-p para derivar como ele se distribui quando a hipótese nula é verdadeira. Para isso, vamos utilizar um truque conhecido da estatística: partindo da variável cuja distribuição conhecemos (neste caso, Z), vamos calcular a função de probabilidade cumulativa do valor-p quando a hipótese nula é verdadeira e, em seguida, computamos sua derivada para obter a função de densidade de probabilidade correspondente.

Começamos definindo valor-p explicitamente como uma variável aleatória P baseada na transformação não-linear da variável Z utilizando a função de probabilidade cumulativa para a normal padronizada. Na definição acima, utilizamos uma valor observado z para calcular o valor-p; aqui, consideramos apenas a transformação da variável aleatória.

\displaystyle P = 2\Phi(-|Z|)

Passamos, então, à derivação da função de probabilidade cumulativa da variável P.

\displaystyle F_{P_{H_0}}(p)=\mathrm{Pr}(P \leq p) = \mathrm{Pr}(2\Phi(-|Z|) \leq p)
\displaystyle =\mathrm{Pr}\left(Z \geq - \Phi^{-1} \left(\frac{p}{2}\right)\right) + \mathrm{Pr}\left(Z \leq \Phi^{-1} \left(\frac{p}{2}\right)\right)
\displaystyle =1 - \Phi\left(-\Phi^{-1}\left(\frac{p}{2}\right)\right) + \Phi\left(\Phi^{-1}\left(\frac{p}{2}\right)\right)
\displaystyle =\frac{p}{2} + \frac{p}{2} = p

ggplot(data.frame(P=c(0, 1)), aes(x=P)) + stat_function(fun=function(x) x) +
  ylab('Pr(P < p)') + xlab('p') + ggtitle('Distribuição Cumulativa de P sob H0')

plot of chunk unnamed-chunk-82

Portanto, a função de probabilidade cumulativa do valor-p é igual ao próprio argumento da função, lembrando que seu domínio é restrito ao intervalo [0, 1]. Isso significa que, se a hipótese nula é verdadeira, a probabilidade de obter um valor-p igual ou menor a 0.05 é igual a 0.05; um valor-p de 0.1 ou menor tem probabilidade 0.1, e assim sucessivamente. Derivar a função de densidade de probabilidade é simples, uma vez que a derivada de f(x) = x é simplesmente f'(x) = 1. A distribuição do valor-p, portanto, segue uma distribuição uniforme quando a hipótese nula é verdadeira.

ggplot(data.frame(P=c(0, 1)), aes(x=P)) + stat_function(fun=function(x) 1) +
  ylab(expression(f[P])) + xlab('p') + ggtitle('Distribuição de Densidade de P sob H0')

plot of chunk unnamed-chunk-83

O fato do valor-p possuir uma distribuição uniforme sob a hipótese nula traz uma série de consequências importantes para seu uso inferencial. A primeira delas é que, se estamos lidando com um caso em que a hipótese nula pontual faz sentido e é verdadeira, qualquer valor-p é igualmente provável. O ponto, porém, é que a probabilidade de obter um valor-p dentro de um intervalo específico é igual ao tamanho deste intervalo. Em particular, a probabilidade de obter um valor-p inferior ao nível de significância \alpha eleito para o teste é exatamente igual a \alpha.

Portanto, a probabilidade de decidirmos contra a hipótese nula no caso de ela ser verdadeira é igual ao nível de significância escolhido. Esse Erro de Tipo I não pode ser aplicado diretamente à interpretação de um único valor-p, contudo, justamente porque ele é uma variável aleatória. Assim como posso ter observado um valor-p igual a 0.01, se a hipótese nula for verdadeira, poderia igualmente ter observado um valor-p igual a 0.8. O que dá ao valor-p a possibilidade de decidir contra a hipótese nula é a definição da regra de decisão baseada no tamanho \alpha do teste: (p < \alpha) \implies \neg H_0. Como a massa de probabilidade entre 0 e \alpha é igual a \alpha, a probabilidade de cometer um erro do Tipo I é limitada pelo tamanho de \alpha.

Distribuição do valor-p quando a hipótese nula é falsa

O que acontece com a distribuição do valor-p quando a hipótese nula é falsa? É importante ressaltar que, no caso em questão, a hipótese alternativa, por afirmar apenas a diferença entre médias, será necessariamente verdadeira quando a hipótese nula for falsa. Mas podemos ser mais restritivos na definição da hipótese alternativa, o que significa que a negação da hipótese nula não implicará, automaticamente, na veracidade da hipótese alternativa. Esse é um outro ponto frágil do uso de testes de hipótese para avaliar teorias científicas: com uma hipótese (estatística) ampla o suficiente, a negação da hipótese nula pode ser utilizada para corroborar uma teoria (científica) mesmo sob evidências muito fracas.

Para podermos computar a distribuição dos valores-p em condição em que a hipótese nula é falsa, precisamos primeiro verificar o que acontece com a estatística de teste Z. Em seu cálculo, ainda consideramos o valor proposto pela hipótese nula – diferença igual a zero – ma agora vemos como isso influencia da distribuição da variável aleatória se a diferença verdadeira não for igual a zero.

Como não estamos mais subtraindo o valor verdadeiro da estimativa da diferença entre os dois grupos, a distribuição ficará centra no valor da diferença verdadeira. Ao padronizarmos pelo erro-padrão, a variância continuará igual a um, mas a média será deslocada por um fator igual a \frac{1}{\sqrt{\frac{\sigma^2_A}{n_A} + \frac{\sigma^2_B}{n_B}}} – que, no caso das variâncias serem iguais entre grupos e com o mesmo tamanho de amostra para ambos, pode ser reduzido a \frac{\sqrt{n}}{\sqrt{2\sigma^2}}. Portanto, a distribuição de Z quando a hipótese nula é falsa, sob as mesmas condições de regularidade definidas acima para a hipótese nula, é dada por:

\displaystyle Z \sim \mathrm{Normal}\left(\frac{\Delta\mu\sqrt{n}}{\sqrt{2\sigma^2}}, 1; \mu_A - \mu_B = \Delta\mu\right)

plotZ2 <- data.frame(Z=seq(-3.5, 5.62, length.out=200),                     DensidadeH0=dnorm(seq(-3.5, 5.62, length.out=200)),                     DensidadeHA=dnorm(seq(-3.5, 5.62, length.out=200), 2.12, 1)) ggplot(plotZ2, aes(x=Z, y=DensidadeH0)) + geom_line(linetype=2) + geom_line(aes(x=Z, y=DensidadeHA)) +   xlab('Z') + ylab('Densidade') + geom_text(aes(x=0, y=0.2, label='Z sob H0')) +   geom_text(aes(x=2.12, y=0.2, label='Z sob HA')) +   geom_segment(aes(x=qnorm(1-0.025), xend=qnorm(1-0.025), y=0, yend=dnorm(qnorm(1-0.025))), linetype=2)+   geom_segment(aes(x=qnorm(0.025), xend=qnorm(0.025), y=0, yend=dnorm(qnorm(0.025))), linetype=2)+   geom_ribbon(data=subset(plotZ2, Z>qnorm(1-0.025)), aes(x=Z, ymin=0, ymax=DensidadeHA), alpha=0.5, fill='skyblue') + geom_text(aes(x=2.12, y=0.2, label='Z sob HA')) +
  geom_ribbon(data=subset(plotZ2, Z>qnorm(1-0.025)), aes(x=Z, ymin=0, ymax=DensidadeH0), alpha=0.5, fill='midnightblue') +
  geom_ribbon(data=subset(plotZ2, Z<qnorm(0.025)), aes(x=Z, ymin=0, ymax=DensidadeH0), alpha=0.5, fill='midnightblue')

plot of chunk unnamed-chunk-84

A figura acima ilustra como o valor-p calculado para a distribuição de Z sob a hipótese nula, quando observamos z=1.96, implica em uma massa de probabilidade bastante diferente quando a hipótose alternativa é verdadeira para um determinado valor da diferença entre as médias.

Como não sabemos que a hipótese nula é falsa, continuamos trabalhando como se ela fosse verdadeira e utilizamos a distribuição cumulativa da normal padronizada para calcular o valor-p, conforme já definimos acima. Retomamos a partir da terceira linha, lembrando agora que a distribuição cumulativa de Z não pode mais ser representada diretamente pela distribuição cumulativa normal padrão.

\displaystyle F_{P_{H_A}} = \mathrm{Pr}(P \leq p)=1 - F_Z\left(-\Phi^{-1}\left(\frac{p}{2}\right)\right) + F_Z\left(\Phi^{-1}\left(\frac{p}{2}\right)\right)

Se centralizarmos a variável Z em zero, porém, podemos voltar a utilizar a distribuição cumulativa da normal padronizada. Com isso, obtemos a função de potência para a diferença entre médias.

\displaystyle =1 - \Phi\left(-\Phi^{-1}\left(\frac{p}{2}\right) - \frac{\Delta\mu\sqrt{n}}{\sqrt{2\sigma^2}}\right) + \Phi\left(\Phi^{-1}\left(\frac{p}{2}\right) - \frac{\Delta\mu\sqrt{n}}{\sqrt{2\sigma^2}}\right)

Por fim, para encontrar a função de densidade de probabilidade do valor-p quando a hipótese nula é falsa, basta derivar a função acima, aproveitando-se dos fatos de que a derivada de uma função de probabilidade cumulativa é igual à função de densidade de probabilidade daquela variável; e que a derivada de uma função inversa é dada por: f^{-1\prime}(x) = \frac{1}{f'(f^{-1}(x))}. Com alguma álgebra, obtemos:

\displaystyle f_{P_{H_A}}(p) = \frac{\phi\left(-\Phi^{-1}\left(\frac{p}{2}\right) - \frac{\Delta\mu\sqrt{n}}{\sqrt{2\sigma^2}}\right) + \phi\left(\Phi^{-1}\left(\frac{p}{2}\right) - \frac{\Delta\mu\sqrt{n}}{\sqrt{2\sigma^2}}\right)}{2\phi\left(\Phi^{-1}\left(\frac{p}{2}\right)\right)}

Não é difícil perceber que essa função, mais genérica, é igual a 1 para qualquer valor-p quando \Delta\mu é igual a zero – incluindo, portanto, a conclusão que observamos acima para o caso especial da hipótese nula. Outro resultado importante é que a distribuição do valor-p depende da magnitude da diferença das médias e de seu erro padrão – que, por sua vez, depende do tamanho da amostra e da variância da variável aleatória de interesse nos dois grupos.

Mas é difícil desenvolver uma intuição clara sobre a forma da função encontrada apenas por sua fórmula. Vamos experimentar plotar algumas linhas para observar o que acontece com a distribuição do valor-p à medida em que a magnitude da diferença entre as médias populacionais aumenta. Para simplificar o gráfico, vamos considerar que a variável de interesse possui a mesma variância nos dois grupos e que esse valor é igual a 1. Assim, qualquer diferença entre as médias pode ser diretamente interpretado como uma magnitude de efeito padronizado em termos do número de desvios-padrão que separam as duas médias. Também vamos assumir que o tamanho da amostra é constante e igual para os dois grupos.

pDens2 <- function(p.value,
                   muA,
                   muB=0,
                   sigmaA=1,
                   sigmaB=1,
                   nA=10,
                   nB=10) {
  deltaMu <- muA-muB
  seMu <- sqrt(sigmaA^2/nA + sigmaB^2/nB)
  z <- deltaMu/seMu
  (dnorm(-qnorm(p.value/2) - z) +
    dnorm(qnorm(p.value/2) - z))/(2*dnorm(qnorm(p.value/2)))
}

p=seq(0.01, 1, length.out=100)
musA <- c(0, 0.3, 0.5, 0.8, 1.0)
muDens <- data.frame(sapply(musA, function(x) pDens2(p, x)))
plotP <- data.frame(p=rep(p, times=length(musA)),
                    Densidade=stack(muDens)$values,
                    Magnitude=factor(rep(musA, each=length(p))))
ggplot(plotP, aes(x=p, y=Densidade, color=Magnitude)) +
  geom_line() + geom_vline(xintercept = 0.05, linetype=2)

plot of chunk unnamed-chunk-85

O resultado é bastante evidente: à medida em que a diferença entre as médias aumenta, os valores-p mais baixos aumentam rapidamente de densidade. Isso implica claramente que um valor-p baixo é muito mais comum quando a hipótese nula é falta. A partir daí fica mais fácil compreender a importância da avaliação da potência estatística de um experimento. Se calcularmos a área sob a curva entre 0 e uma valor de significância \alpha de interesse (como 0.05, marcado no gráfico), podemos identificar, a priori qual a probabilidade de obter um valor-p inferior ao nível de significância para cada magnitude de efeito plausível.

Por fim, um gráfico de perspectiva avaliando o efeito contínuo do aumento da diferença entre as médias populacionas sobre a distribuição do valor-p – considerando que o erro-padrão permanece constante.

m <- seq(0, 0.5, length.out=20)
p <- seq(0.01, .5, length.out=20)
z <- outer(p, m, FUN = pDens2)
persp(p, m, z, col='skyblue', border='steelblue',
      theta=60, phi=20, zlim=c(0,max(z)),
      xlab='p', ylab='Magnitude', zlab='Densidade')

plot of chunk unnamed-chunk-86

Uma proposta unificada para estimação da fidedignidade de uma escala

Uma proposta unificada para estimação da fidedignidade de uma escala

O problema

Qualquer procedimento de mensuração de uma variável incorre, invariavelmente, em alguma forma de erro. Nas ciências naturais, a proposta de modelos para os erros de medida levou à descoberta de uma série de modelos estatísticos cruciais, como a distribuição normal, e permitiu qualificar a avaliação dos resultados de experimentos empíricos em relação a valores preditos pelas teorias.

Para a psicologia, o processo de mensuração costuma esbarrar em mais um problema: muitas das variáveis de interesse do pesquisador são latentes, ou seja, não podem ser diretamente observadas – seja por circunstância do processo de medida ou por definição. Portanto, qualquer medida em psicologia precisa estabelecer com clareza a relação entre a variável latente de interesse e as variáveis indicadoras utilizadas para inferir seu valor.

Uma propriedade psicométrica costumeiramente reportada para avaliar a qualidade de um instrumento é a fidedignidade (denotada como \rho_{xx}), definida na Teoria Clássica dos Testes como a proporção da variância da escala atribuível à variância da variável latente. Apesar da definição possuir propriedades atraentes, como ser uma razão interpretável como proporção e não depender da unidade original da escala, estimá-la a partir de um conjunto de dados não é uma tarefa simples.

Sob alguns pressupostos restritivos, a fidedignidade pode ser estimada a partir da correlação entre aplicações repetidas da escala para uma mesma amostra. O pressuposto de estabilidade temporal dos escores verdadeiros e da irrelevância do efeito de memória, todavia, são difíceis de serem aceitos para a maioria das escalas psicométricas. Além disso, é rara a situação em que o pesquisador dispõe de duas aplicações de um instrumento, o que levou ao desenvolvimento de estimadores baseados em uma única amostra.

A pesquisa sobre estimadores internos à amostra gerou um enxurrada de alternativas, baseadas em pressupostos diferentes e com propriedades dependentes desses pressupostos. O mais famoso desses coeficientes é o \alpha de Cronbach, proposto originalmente por Guttman, e rotineiramente reportado na avaliação de propriedades psicométricas de escalas. Uma motivação para interpretação do \alpha como fidedignidade do escore composto da escala a partir do pressuposto de testes paralelos já foi tema do blog, além de alguns ditirambos por seu uso indiscriminado (para outras críticas, cf., p.e., Green & Yang, 2008). Muitos desses estimadores são denominados de medidas de “consistência interna” da escala. Porém, ao contrário da noção de fidedignidade, a noção de consistência interna não é bem definida e, por isso, seu uso é desencorajado (McDonald, 1999).

Uma proposta de solução

Frente à profusão e limitação de estimadores para a fidedignidade de uma escala – ou mesmo diversas confusões sobre sua interpretação – como o pesquisador pode proceder? Recentemente, cruzei com um artigo da autoria de Alonso, Laenen, Molenberghs, Geys & Vangeneugden (2010) que propõe uma família de estimadores que busca unificar, a partir de uma série de desiderata, diversos coeficientes da fidedignidade de uma escala. Esses desiderata são:

  1. 0 \leq R \leq 1;
  2. Se R = 0, a escala não traz nenhum informação sobre o escore verdadeiro;
  3. Se R = 1, a escala é uma função linear determinística do escore verdadeiro;
  4. Sob o pressuposto de testes paralelos, R = \frac{\sigma^2_{\tau}}{\sigma^2_{\tau} + \sigma^2_{\epsilon}}.

O apelo dessas propriedades é imediatamente óbvio: além de manter a definição de fidedignidade da TCT como caso específico, preserva a interpretação intuitiva de um coeficiente sem dimensão como proporção de variância explicada.

A proposta de Alonso et al. também amplia a noção de fidedignidade para o caso de escalas com mais de um fator latente. Além disso, possui o mérito de diferenciar a fidedignidade da escala, tomada como um um conjunto multivariado de itens, da fidedignidade de um escore baseado em alguma transformação linear dos itens da escala – em geral, a soma ou média dos itens.

Para este post, vamos tentar compreender, por meio de simulações e comparações com os estimadores clássicos de fidedignidade, três estimadores diferentes:

\displaystyle R_T = 1 - \frac{\mathrm{tr}(\Theta)}{\mathrm{tr}(\Sigma)}
\displaystyle R_{\Lambda} = 1 - \frac{\mathrm{det}(\Theta)}{\mathrm{det}(\Sigma)}
\displaystyle \rho(a) = 1 - \frac{\mathbf{a}'\Theta\mathbf{a}}{\mathbf{a}'\Sigma\mathbf{a}}

O primeiro estimador multivariado da fidedignidade, R_T, é baseado apenas nas diagonais da matriz de covariância dos erros (\Theta) e na matriz de covariância das variáveis (\Sigma). Por isso, ela não é influenciada pela presença de correlações entre erros nem pelos elementos para fora da diagonal da covariância dos indicadores.

O segundo estimador multivariado, R_{\Lambda}, é baseado na determinante das matrizes de covariância dos erros e das variáveis, e leva em consideração os elementos para fora da diagonal.

O último estimador, \rho(a), avalia a fidedignidade de um escore composto a partir de uma combinação linear dos indicadores da escala. A maneira mais comum de se combinar indicadores é por meio da soma ou média dos itens ou um subconjunto deles.

Mas por que se importar?

Da maneira como as medidas de fidedignidade são rotineiramente reportadas e interpretadas em estudos psicométricos, não há nenhum problema imediato. Quando se começa a fazer uso desses coeficientes para análises posteriores, porém, as dificuldades se tornam imediatamente aparentes.

Um primeiro problema é relacionado ao cálculo do erro-padrão da medida como uma estimativa da incerteza sobre o escore efetivo de um sujeito em uma avaliação. Uma estimativa rápida para o erro de medida é \sigma_x\sqrt{1 - \rho_{xx}}. Se subestimarmos demais a fidedignidade da escala, o erro será superestimado e dificultará a comparação entre sujeitos em função da sobreposição dos intervalos de confiança; no caso contrário, subestimar o erro levará a decidir erroneamente pela diferença entre escores próximos.

Um segundo problema, imediatamente relacionado com o primeiro, é a utilização da estimativa da fidedignidade para desatenuar a correlação entre escores de escalas psicométricas e outras variáveis. Na presença de erro de medida, a correlação entre duas variáveis será enviesada para baixo, ou seja, será atenuada. A estimativa de fidedignidade pode ser utilizada para estimar a correlação verdadeira, \rho_{xy} = \frac{\rho_{x'y'}}{\sqrt{\rho_{x'x'}\rho_{y'y'}}}. Subestimar a fidedignidade irá, invariavelmente, superestimar a correlação verdadeira. Como o \alpha de Cronbach sistematicamente subestima a fidedignidade em escalas congenéricas, um estudo poderá reportar correlações enviesadas para cima se utilizar esse coeficiente para corrigir as correlações.

Fazer uso de estimadores eficientes para a fidedignidade e interpretá-los corretamente, portanto, é crucial para estimar corretamente a incerteza de uma escala psicométrica e avaliar seu impacto na prática profissional e na realização de pesquisas.

Algumas definições

Para avaliar a proposta de Alonso et al., é preciso recordar a definições de alguns conceitos importantes para a história da teoria clássica de testes. Os primeiros dizem respeito a alguns pressupostos que definem diferentes classes de itens ou testes: testes paralelos, teste \tau-equivalentes e teste congenéricos.

Testes paralelos

O pressuposto de testes paralelos pode ser definido da seguinte maneira:

\displaystyle \mathbf{X}_i = \boldsymbol{\tau} + \boldsymbol{\epsilon}_i
\displaystyle \sigma^2_{\epsilon_i}=\sigma^2_{\epsilon}, \forall i\in {1, \dots,N}

Ou seja, os valores observados do teste i são centrados no escore verdadeiro e dispersos em função de um erro aleatório. Há apenas um construto avaliado pelos testes, sem nenhum viés multiplicativo e com a mesma variância dos erros para todos os itens. Também estamos pressupondo que não há covariância entre os erros e escores verdadeiros; e que não há nenhum viés aditivo. Esse último pressuposto pode ser flexibilizado, pois o acréscimo de constantes à equação não altera a covariância entre os itens.

Testes τ-equivalentes

O pressuposto de igualdade de variâncias entre itens permite derivar algumas propriedades interessantes a respeito da fidedignidade de um conjunto de indicadores, como discutimos no post passado, e serviu de base para vários avanços importantes na TCT. Um primeiro passo na flexibilização desses pressupostos é permitir que a variância dos erros possa diferir entre os itens. O viés multiplicativo do escore verdadeiro continua sendo o mesmo para todos os itens, por isso o nome:

\displaystyle \mathbf{X}_i = \boldsymbol{\tau} + \boldsymbol{\epsilon}_i
\displaystyle \mathrm{var}[\epsilon_i] = \sigma^2_{\epsilon_i}
\displaystyle \mathrm{cov}[\epsilon_i, \epsilon_j] = 0, \forall i,j \in {1 \dots N} \land i \neq j

A última propriedade também deveria estar presente na definição dos testes paralelos: a correlação existente entre indicadores se deve apenas ao fator latente e, depois de levá-lo em consideração, os erros entre os itens não possuem nenhuma covariância adicional – ou seja, os erros são independentes.

Testes congenéricos

Relaxando o pressuposto de equivalência do fator multiplicativo do escore verdadeiro, avançamos mais um passo na flexibilização do modelo original dos testes paralelos. Todos os outros pressupostos relativos aos testes \tau-equivalentes permanecem os mesmos, inclusive a independência dos erros, mas agora acrescentamos um viés multiplicativo específico a cada item:

\displaystyle\mathbf{X}_i = \lambda_i\boldsymbol{\tau} + \boldsymbol{\epsilon}_i

Avaliando a proposta de Alonso et al.

Para compreendermos melhor o que Alonso et al. estão propondo, vamos começar definindo algumas funções auxiliares. A primeira delas, compRel, calcula os coeficientes de fidedignidade propostos no artigo e retomados acima a partir dos dados da população. A segunda, alonsoEmpRel calcula os mesmos coeficientes, mas agora a partir de um conjunto de dados simulados com auxílio do pacote lavaan e utilizando um modelo fatorial para computar a matriz de variâncias residuais. Para as simulações, consideramos as variáveis latentes como possuindo expectância zero e variância um. Com o objetivo de obter estimativas pouco influenciadas pelo erro amostral, vamos simular amostras grandes, com 5.000 casos.

tr <- function(M){
  sum(diag(M))
}

compRel <- function(errorVar,
                    loadings,
                    factorVar=matrix(1),
                    a=rep(1, length(errorVar))){
  sigma <- diag(errorVar)
  Var <- loadings %*% factorVar %*% t(loadings) + sigma
  if (is.vector(a)){
    rho1 <- 1 - (t(a)%*%sigma%*%a)/(t(a)%*%Var%*%a)
  } else {
    rho1 <- t(as.matrix(apply(a, 1, function(x)
      1 - (t(x) %*% sigma %*% x/(t(x)%*%Var%*%x)))))
    colnames(rho1) <- paste('rho', 1:nrow(a), sep='')
  }
  cbind(data.frame(Rt=(1 - tr(sigma)/tr(Var)),
                   Rl=(1 - det(sigma)/det(Var))),
        rho1)
}

alonsoEmpRel <- function(lavaanModel,
                         a=rep(1, nrow(inspect(lavaanModel)$lambda))){
  pars <- inspect(lavaanModel, what='est')
  implied <- lavaanModel@implied$cov[[1]]
  resVar <- pars$theta
  if (is.vector(a)){
    rho1 <- 1 - (t(a)%*%resVar%*%a)/(t(a)%*%implied%*%a)
  } else {
    rho1 <- t(as.matrix(apply(a, 1, function(x)
      1 - (t(x) %*% resVar %*% x/(t(x)%*%implied%*%x)))))
    colnames(rho1) <- paste('rho', 1:nrow(a), sep='')
  }
  cbind(data.frame(Rt = (1 - tr(resVar)/tr(implied)),
                   Rl = (1 - det(resVar)/det(implied))),
        rho1)
}

Testes paralelos

Com auxílio do lavaan, simulamos cinco mil observações de quatro indicadores de uma variável latente. Como iniciamos com os testes paralelos, a carga fatorial do modelo é igual a um para os quatro itens e a variância do erro é a igual a 0.25. Utilizamos a função reliability do pacote semTools para calcular três coeficientes clássicos de fidedignidade com base no modelo proposto: o \alpha de Cronbach, o \omega_T de McDonald e a Variância Média Extraída (Average Variance Extracted, AVE), proposta por Fornell & Larcker (1981).

Como sabemos os valores das variâncias do fator latente e dos erros, podemos calcular, manualmente, a fidedignidade para um único item da escala – 0.8 – ou para um escore composto pela soma simples dos quatro itens – 0.94.

library(lavaan)
library(semTools)

parallel <- '
f =~ 1*i1 + 1*i2 + 1*i3 + 1*i4

f ~ 0*1
f ~~ 1*f

i1 ~~ 0.25*i1
i2 ~~ 0.25*i2
i3 ~~ 0.25*i3
i4 ~~ 0.25*i4
'
parallelData <- simulateData(parallel, model.type=cfa, sample.nobs=5000)

out <- rbind(compRel(rep(0.25, 4), rep(1, 4)),
             alonsoEmpRel(cfa('f=~i1+i2+i3+i4', parallelData)))
rownames(out) <- c('População', 'Estimado')
knitr::kable(out, digits=3)
Rt Rl rho1
População 0.800 0.941 0.941
Estimado 0.792 0.938 0.938
knitr::kable(reliability(cfa('f=~i1+i2+i3+i4', parallelData))[c(1, 2, 5), ], digits=3)
f total
alpha 0.938 0.938
omega 0.938 0.938
avevar 0.792 0.792

Para escalas compostas por itens paralelos, a fidedignidade multivariada do conjunto específico de indicadores é exatamente igual à fidedignidade da soma dos itens. O \alpha e o \omega possuem o mesmo valor, também idênticos ao R_{\Lambda} e \rho(1). Conforme o quarto desideratum proposto por Alonso et al., os estimadores sugeridos podem ser reduzidos à definição clássica de fidedignidade. O coeficiente R_T, por sua vez, corresponde à fidedignidade de um item específico e é equivalente à AVE.

Testes τ-equivalentes

Passamos, em seguida, à avaliação dos estimadores para testes \tau-equivalentes, considerando as seguintes variâncias para os quatro itens da escala: 0.1., 0.2, 0.5 e 1. Como as variâncias do erro são diferentes para cada indicador, não temos apenas um valor para a fidedignidade dos itens, mas quatro: 0.91, 0.83, 0.67 e 0.5 – diminuindo à medida em que aumenta a variância residual. Podemos argumentar em favor de um valor da fidedignidade média dos itens, utilizando a média das variâncias residuais junto com a fórmula de fidedignidade para testes paralelos. Essa definição nos dá a fidedignidade média dos itens como 0.69. Considerando também a média das variâncias residuais para calcular a fidedignidade verdadeira da soma dos itens da escala, obtemos 0.89.

tau <- '
f =~ 1*i1 + 1*i2 + 1*i3 + 1*i4
f ~ 0*1
f ~~ 1*f
i1 ~~ 0.1*i1
i2 ~~ 0.2*i2
i3 ~~ 0.5*i3
i4 ~~ 1*i4
'
tauData <- simulateData(tau, model.type=cfa, sample.nobs = 5000)
out <- rbind(compRel(c(0.1, 0.2, 0.5, 1.0), rep(1, 4)),
             alonsoEmpRel(cfa('f=~i1+i2+i3+i4', tauData)))
rownames(out) &lt;- c('População', 'Estimado')
knitr::kable(out, digits=3)
Rt Rl rho1
População 0.69 0.947 0.899
Estimado 0.69 0.945 0.899
knitr::kable(reliability(cfa('f=~i1+i2+i3+i4', tauData))[c(1, 2, 5), ], digits=3)
f total
alpha 0.899 0.899
omega 0.899 0.899
avevar 0.690 0.690

Aqui começamos a perceber as primeiras diferenças entre os estimadores propostos por Alonso et al. e os estimadores tradicionais que estamos considerando. Um primeiro ponto de interesse é que a fidedignidade multivariada do conjunto de itens, R_{\Lambda}, é maior do que a fidedignidade da escore composto pela soma dos itens. Por um lado, esse resultado faz sentido: ao utilizarmos um combinação linear das variáveis em vez de toda a matriz de dados, estamos jogando informação fora em favor da simplicidade de trabalhar com uma única variável.

Outro aspecto importante é que o \alpha e o \omega possuem o mesmo valor, que corresponde à fidedignidade do escore da escala, \rho(1). Portanto, se formos interpretar ao pé da letra, os dois estimadores tradicionais da fidedignidade dizem respeito ao escore produzido pela escala e não à escala como um todo. Utilizar pesos diferentes para construir o escore também significa que o \alpha já não é mais um estimador adequado da fidedignidade.

Por fim, a AVE corresponde, novamente, ao estimador R_T, referendando sua interpretação como fidedignidade (média) dos itens tomados de forma isolada.

Testes congenéricos

Para os testes congenéricos, mantemos as variâncias residuais das simulações acima, mas acrescentamos cargas fatoriais diferentes de um: 0.6, 1.2, 0.5 e 0.3. Novamente, a fidedignidade de cada item varia em função de sua carga fatorial e variância residual: 0.98, 0.97, 0.67 e 0.08. Como definimos que a variância da variável latente é igual a um, podemos simplesmente utilizar a média dos quadrados das cargas fatoriais (sua comunalidade) e a média das variâncias residuais, em conjunto com a fórmula da fidedignidade para testes paralelos, para estimarmos a fidedignidade média dos itens como 0.54. Da mesma forma, para um escore composto pela soma dos itens, a fidedignidade verdadeira é 0.79.

con <- '
f =~ 0.6*i1 + 1.2*i2 + 0.5*i3 + 0.3*i4
f ~ 0*1
f ~~ 1*f
i1 ~~ 0.1*i1
i2 ~~ 0.2*i2
i3 ~~ 0.5*i3
i4 ~~ 1*i4
'
conData <- simulateData(con, model.type=cfa, sample.nobs = 5000)

out <- rbind(compRel(c(0.1, 0.2, 0.5, 1.0), c(0.6, 1.2, 0.5, 0.3)),
             alonsoEmpRel(cfa('f=~i1+i2+i3+i4', conData)))
rownames(out) &lt;- c('População', 'Estimado')
knitr::kable(out, digits=3)
Rt Rl rho1
População 0.543 0.919 0.79
Estimado 0.545 0.920 0.79
knitr::kable(reliability(cfa('f=~i1+i2+i3+i4', conData))[c(1, 2, 5), ], digits=3)
f total
alpha 0.720 0.720
omega 0.790 0.790
avevar 0.545 0.545

O \alpha se mostra inadequado como estimativa da fidedignidade do escore da escala, possuindo um viés que subestima o valor verdadeiro. O coeficiente \omega é equivalente ao \rho(1), indicando sua adequação para estimar a fidedignidade de um escore, mas novamente subestima a fidedignidade do conjunto de itens como um todo, tal qual mensurado pelo R_{\Lambda}. O coeficiente R_T, por fim, continua igual à AVE e a fidedignidade média dos itens constitutivos da escala.

Testes multifatoriais (e congenéricos)

Para finalizar, consideramos uma situação comum no desenvolvimento de escalas: a presença de mais de um fator latente responsável pela covariância entre os itens. A proposta de mais de um fator pode se dar tanto pela construção da escala a partir de um construto multidimensional ou pode se mostrar necessária a partir da avaliação empírica das correlações entre itens. Independente de sua motivação, vamos considerar uma escala multifatorial composta por duas variáveis latentes positivamente correlacionadas, cada uma mensurada por três itens com as mesmas cargas fatoriais e variâncias residuais.

A presença de dois fatores torna mais complicada a noção de fidedignidade. Como construir um escore? A resposta mais usual é calcular escores separados para cada dimensão da escala, o que reduz o problema a avaliar a fidedignidade do escore de subescalas. Mas, com isso, perdemos a economia de resumir toda a escala a um único escore. Nada nos impede, porém, de construir um escore único, mesmo que não suponhamos a existência de um fator geral que explicaria a covariância entre fatores. Nesse caso, contudo, como avaliar a fidedignidade?

As complicações do caso geral (e, certamente, mais comum!) mostram a relevância do trabalho de Alonso et al. Ao considerar a natureza multivariada da escala, a fidedignidade continua com sua interpretação de variância explicada, mas extendida para duas interpretações possíveis do que constitui a variância no caso de um vetor de variáveis aleatórias: a fidedignidade da variância total, definida pelo traço das matrizes de covariância e que constitui o R_T, ou a fidedignidade da variância generalizada, baseada na determinante das matrizes de covariância e empregada pelo R_{\Lambda}.

Neste caso, vamos considerar a fidedignidade dos escores baseados em cada fator latente e um escore baseado na soma de todos os itens da escala.

multi <- '
f1 =~ 0.6*i1 + 1.2*i2 + 0.5*i3
f2 =~ 0.6*i4 + 1.2*i5 + 0.5*i6

f1 ~ 0*1
f1 ~~ 1*f1
f2 ~ 0*1
f2 ~~ 1*f2
f1 ~~ 0.5*f2

i1 ~~ 0.1*i1
i2 ~~ 0.2*i2
i3 ~~ 0.5*i3

i4 ~~ 0.1*i4
i5 ~~ 0.2*i5
i6 ~~ 0.5*i6

'
multiData <- simulateData(multi, model.type=cfa, sample.nobs = 5000)
out <- rbind(compRel(errorVar=rep(c(0.1, 0.2, 0.5), 2),
                     loadings = matrix(c(0.6, 1.2, 0.5, rep(0, 3), rep(0, 3), 0.6, 1.2, 0.5), 6, 2),
                     factorVar=matrix(c(1, 0.5, 0.5, 1), 2, 2, byrow=T),
                     a=rbind(c(1, 1, 1, 0, 0, 0),
                             c(0, 0, 0, 1, 1, 1),
                             c(1, 1, 1, 1, 1, 1))),
             alonsoEmpRel(cfa('f1=~i1+i2+i3 \n f2=~i4+i5+i6', multiData),
                          a=rbind(c(1, 1, 1, 0, 0, 0),
                                  c(0, 0, 0, 1, 1, 1),
                                  c(1, 1, 1, 1, 1, 1))))
rownames(out) <- c('População', 'Estimado')
knitr::kable(out, digits=3)
Rt Rl rho1 rho2 rho3
População 0.719 0.992 0.869 0.869 0.908
Estimado 0.719 0.992 0.864 0.873 0.908
knitr::kable(reliability(cfa('f1=~i1+i2+i3 \n f2=~i4+i5+i6', multiData))[c(1, 2, 5), ], digits=3)
f1 f2 total
alpha 0.791 0.804 0.805
omega 0.864 0.873 0.908
avevar 0.712 0.726 0.719

Comecemos pelo aspecto que permaneceu constante ao longo de todas as variações examinadas. O coeficiente R_T é, novamente, equivalente à variância média extraída total. Uma das principais vantagens da utilização de ambos os coeficientes é o fato de sua estimativa não ser influenciada por correlações negativas entre fatores. Sua interpretação continua sendo a fidedignidade média dos itens da escala, considerando todos os fatores presentes.

O coeficiente multivariado de fidedignidade baseado na determinante, R_{\Lambda} é consideravelmente mais alto que qualquer fidedignidade basedo em escores derivados da soma dos itens. Esse fato é relevante porque mostra que mesmo uma escala bastante fidedigna pode gerar escores sumários com propriedades psicométricas mais pobres.

Na tabela acima, as colunas rho1, rho2 e rho3 indicam, respectivamente, a fidedignidade da soma dos três primeiros, dos três últimos e de todos os itens. Como os itens possuem os mesmos coeficientes para os dois fatores, a fidedignidade verdadeira é a mesma para as duas subescalas; e a fidedignidade do escore total é superior à fidedignidade das subescalas. As estimativas da fidedignidade dos diferentes escores são as mesmas tanto para o estimador proposto por Alonso et al., \rho(1) quanto para \omega de McDonald.

O \alpha, porém, subestima a fidedignidade dos subescores e do escore geral, sendo especialmente ruim para esse último. Apesar de não aparecer nos resultados acima, o \alpha é ainda pior se os fatores são ortogonais ou se os fatores são negativamente correlacionados e a direção dos itens não é invertida. Também não simulamos acima o caso de covariâncias entre os erros – nesse caso, a estimativa do \alpha pode ser inflacionada.

Conclusões

  1. As diferentes medidas de fidedignidade tradicionais não são, estritamente falando, uma estimativa da fidedignidade da escala tomada de maneira multivariada, mas a fidedignidade de um escore composto pela soma simples das respostas aos indicadores individuais.
  2. No caso (improvável) da escala ser composta por itens paralelos ou \tau-equivalentes, o \alpha de Cronbach é uma boa estimativa da fidedignidade do escore baseado na soma, sendo inclusive equivalente à fidedignidade multivariada no caso específico de testes paralelos.
  3. Em todos os casos mais flexíveis (e verossímeis), o \alpha tende a subestimar a fidedignidade do escore da escala, especialmente no caso de escalas multifatoriais e com itens congenéricos. Nesses caso, é ideal utilizar o \omega ou os coeficientes propostos por Alonso et al.
  4. As três definições de fidedignidades propostas por Alonso et al. investigadas acima avaliam aspectos diferentes da variância explicada pelos escores verdadeiros:
    • O R_T avalia a média da proporção de variância explicada por item da escala;
    • O R_{\Lambda} avalia a proporção de variância multivariada de um conjunto específico de itens explicada pela covariância dos escores verdadeiros;
    • O \rho(\mathbf{a}) avalia a fidedignidade de um escore composto por transformações lineares dos valores dos itens.

Por isso, na avaliação das propriedades psicométricas de uma escala, o ideal é apresentar as três medidas, facilitando a avaliação da aplicabilidade do instrumento em diferentes contextos. Um psicometrista pode estar interessado na fidedignidade multivariada do conjunto de indicadores utilizados para avaliar a escala como um todo; um pesquisador pode estar interessado na fidedignidade média por item para avaliar a priori a fidedignidade de escalas derivadas; e quem utiliza a escala em situações práticas e depende do uso de estatísticas sumárias, como a soma dos itens, está particularmente preocupado com a fidedignidade do escore composto.

Uma intuição sobre o α de Cronbach

Uma intuição sobre o α de Cronbach

Nos estudos de avaliação das propriedades psicométricas de escalas psicológicas, uma das estatísticas mais reportadas é o \alpha de Cronbach. Sua interpretação nos estudos costuma variar entre a fidedignidade ou confiabilidade do instrumento; a consistência interna da escala; ou mesmo a dimensionalidade dos indicadores utilizados. Apesar da pluralidade com a qual o \alpha é interpretado – em termos do seu significado para um estudo psicométrico – o entendimento de seu valor numérico costuma ser varrido para debaixo do tapete pelo uso de faixas de referência, p.e.: \alpha > 0.8 = “bom”, \alpha < 0.7 = “ruim”, etc. etc.

Para a compreensão de uma estatística, porém, convém desenvolver pelo menos uma intuição das quantidades que ela relaciona e sumariza. Não há outro caminho senão retomarmos a equação que define o coeficiente \alpha:

\displaystyle \alpha = \frac{n}{n-1} \left( 1 - \frac{\sum_{i=1}^n \mathrm{var}[X_i]}{\mathrm{var}[Y]} \right)

Em que n é o número de itens ou indicadores de um instrumento, \mathrm{var}[X_i] é a variância de um indicador específico X_i e \mathrm{var}[Y] é a variância do total da escala, obtido pela soma de todos os n itens, i.e. \sum_{i=1}^n X_i. Essa é a definição dada por Cronbach em seu artigo original, ainda que com alguma diferença na notação utilizada para indicar a variância.

Para esse post, vamos nos focar em compreender melhor o que significa a razão que se encontra entre parênteses na equação. Para isso, porém, é necessário relembrar algumas propriedades essenciais sobre a variância de variáveis aleatórias.

Revisando algumas propriedades da variância

A variância é uma das medidas mais utilizadas para avaliar a dispersão de uma variável aleatória, ou seja, quão distantes os valores que a variável X pode tomar estão do seu valor esperado, \mathbb{E}[X]. O valor esperado, por sua vez, pode ser interpretado como a média aritmética dos valores da variável aleatória ao longo de uma série longa de ensaios. Para a variância, a “distância” com relação a esse valor esperado é conceituada em termos da diferença quadrática esperada:

\mathrm{var}[X] = \mathbb{E}[(X - \mathbb{E}[X])^2]

Expandindo o binômio e rearranjando os termos, obtemos uma equação mais direta que facilitará as computações posteriores:

\displaystyle \mathrm{var}[X] = \mathbb{E}[(X - \mathbb{E}[X])^2] \\  = \mathbb{E}[X^2 - 2X\mathbb{E}[X] + \mathbb{E}[X]^2] \\  = \mathbb{E}[X^2] - 2\mathbb{E}[X]\mathbb{E}[X] + \mathbb{E}[X]^2 \\  = \mathbb{E}[X^2] - \mathbb{E}[X]^2

Voltando rapidamente ao contexto de escalas psicológicas: o escore de uma escala é definido, na maioria das vezes, como a soma ou a média dos itens utilizados para avaliação do construto. Para facilitar a notação e mantermos a coerência com a equação utilizada por Cronbach, vamos assumir que o escore final é definido simplesmente pela soma dos itens, conforme sugerimos acima. Sem entrar na definição de fidedignidade, uma propriedade psicométrica que o \alpha se propõe mensurar, podemos questionar como a variância do escore se comporta em função da variância dos itens específicos.

Em outras palavras, estamos interessados no que acontece com a variância da soma de variáveis aleatórias e, em particular, como a variância da soma se relaciona com as variâncias das variáveis individuais. Para facilitar a exposição, vamos considerar o caso da soma de duas variáveis aleatórias, X e Y. É importante notar que nenhum pressuposto é feito com relação à forma funcional dessas variáveis e que, portanto, os resultados obtidos não dependem, p.e., da normalidade das distribuições. Retomando a equação acima para o caso da somas das duas variáveis, expandindo os binômios e reorganizando os termos:

\displaystyle \mathrm{var}[X + Y] = \mathbb{E}[(X + Y)^2] - \mathbb{E}[X+Y]^2 \\ = \mathbb{E}[X^2 + 2XY + Y^2] - \mathbb{E}[X+Y]\mathbb{E}[X+Y] \\ = \mathbb{E}[X^2] + 2\mathbb{E}[XY] + \mathbb{E}[Y^2] - \left( (\mathbb{E}[X] + \mathbb{E}[Y]) (\mathbb{E}[X] + \mathbb{E}[Y])\right)\\ = \mathbb{E}[X^2] + \mathbb{E}[Y^2] + 2\mathbb{E}[XY] - \mathbb{E}[X]^2 - 2\mathbb{E}[X]\mathbb{E}[Y] - \mathbb{E}[Y]^2 \\ = \mathbb{E}[X^2] - \mathbb{E}[X]^2 + \mathbb{E}[Y^2] - \mathbb{E}[Y]^2 + 2\left(\mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]\right) \\ = \mathrm{var}[X] + \mathrm{var}[Y] + 2\left(\mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]\right)

Portanto, a variância da soma de duas variáveis aleatórias é sempre igual à soma das variâncias dessas variáveis, mais duas vezes a quantidade indicada entre parênteses, \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]. Essa quantidade é a covariância entre as variáveis aleatórias. Isso nos leva a três quadros possíveis para os resultados:

  1. As variáveis são independentes, o que significa que sua covariância é igual a 0 e, portanto, a variância das somas é igual à soma das variâncias, \mathrm{var}[X+Y] = \mathrm{var}[X] +\mathrm{var}[Y];
  2. As variáveis possuem uma correlação linear positiva, e a variância da soma é maior que a soma das variâncias, \mathrm{var}[X+Y] > \mathrm{var}[X] +\mathrm{var}[Y];
  3. As variáveis possuem uma correlação linear negativa, e a variância da soma é menor que a soma das variâncias, \mathrm{var}[X+Y] < \mathrm{var}[X] +\mathrm{var}[Y].

Retornando ao α

Voltando ao caso das escalas psicométricas, vamos nos focar no caso em que, se existe alguma covariância entre os itens de um instrumento, ela será necessariamente positiva. A equação que define o \alpha certamente assume isso, e é o motivo pelo qual o coeficiente só deve ser calculados com os itens negativos invertidos – caso contrário, a covariância negativa irá reduzir indevidamente o valor do coeficiente. Também vamos considerar que o instrumento é composto por um número grande de itens, o que implica que o primeiro termo da equação \frac{n}{n-1} \approx 1 e pode ser ignorado.

Comecemos com o caso limite: se todos os indicadores de uma escala são independentes entre si (ou, por qualquer outra razão, são linearmente não-correlacionados), a variância de sua soma é igual à soma das variâncias de cada item em particular – ou seja, o quadro (1) nos desfechos acima.

Flexibilizemos um pouco os pressupostos: em vez de completa independência, vamos supor que há pelo menos um par de indicadores com covariância diferente de 0. Qual será a variância da soma desses indicadores? Como assumimos acima que qualquer covariância existente é estritamente positiva, a variância da soma é necessariamente maior do que a soma das variâncias de cada item isolado – agora, o quadro (2) nos desfechos elencados.

Portanto, a soma das variâncias dos indicadores nos permite estabelecer um limite inferior para a variância da soma. Ou seja, do ponto de vista da construção de uma escala, no pior cenário possível, os itens não tem correlação nenhuma entre si e a variância de sua soma será igual à soma de suas variâncias. Em qualquer outra situação, assumindo que as covariâncias são positivas, a variância da soma é sempre maior que a soma das variâncias.

Retomemos a equação que define o \alpha: onde podemos encontrar um termo que corresponda a esse limite inferior? Justamente no numerador da razão colocada entre parênteses na representação acima, \sum_{i=1}^n \mathrm{var}[X_i]. Esse numerador, por ser constituído pela soma das variâncias de cada um dos n itens do instrumento, é igual ao valor mínimo que a variância da soma dos itens pode obter.

O denominador da razão é a variância da soma, \mathrm{var}[Y]. No pior quadro possível – completa independência entre todos os pares de indicadores – esse valor é igual ao do denominador, e o resultado da divisão é igual a 1. Como demonstrado acima, \mathrm{var}[Y] \geq \sum_{i=1}^n \mathrm{var}[X_i], e, portanto, o resultado dessa razão nunca pode ultrapassar 1. Esse limite superior para a razão explica porque ela é deduzida de 1 no termo colocado entre parênteses, e estabelece um limite inferior para o coeficiente: 0. Com uma amostra limitada, o valor do \alpha pode ser negativo, mas isso é apenas um artefato de pequenas amostras.

À medida em que aumentam as covariâncias entre pares de indicadores ou aumenta o número de itens com uma mesma covariância, o denominador fica cada vez maior, o que implica que o resultado da razão diminui e, como o deduzimos de 1, o valor do \alpha aumenta. Como a razão nunca pode resultar em um número negativo, uma vez que a variância sempre é positiva ou igual a zero, o coeficiente \alpha se aproxima de 1 quando o número de itens ou a covariância entre indicadores aumenta. Ao contrário do limite inferior, que pode ocorrer no caso (improvável) da ausência de covariância entre todas as variáveis, o limite superior só pode ser obtido assintoticamente, se desconsiderarmos o primeiro termo da equação.

Entrando em cena a Teoria Clássica dos Testes

Até agora, nossa intuição permitiu estabelecer os limites inferior e superior do \alpha de Cronbach, ignorando o primeiro termo da equação. Essa intuição, porém, ainda não permite apreciar com cuidado a magnitude do valor do coeficiente. Um \alpha de 0.5, por exemplo, significa apenas que a soma das covariâncias entre pares de indicadores é igual à metade da soma das variâncias.

Em nenhum lugar, porém, estabelecemos uma explicação para a covariância entre os pares de indicadores. Pelo contrário, poderíamos conceber um conjunto de itens que possui alta covariância entre seus pares, mas que cada covariância é determinada pela relação única entre cada par, sem implicar na existência de um fator comum. Nesse caso, obteríamos um \alpha alto mesmo sem pressupormos a existência de um fator latente! Ou, então, alguns poucos pares com forte covariância poderiam ser responsáveis por um aumento no coeficiente.

Do ponto de vista da construção de escalas, porém, esses casos extremos não nos interessam – ainda que não possam ser esquecidos na explicação da correlação entre itens. Em seu lugar, assumimos alguns pressupostos a respeito do que causa a covariância entre os indicadores. Para não complicarmos muito, vamos partir de algumas definições da Teoria Clássica dos Testes (TCT). Em particular, vamos considerar que o valor de um indicador ou item corresponde ao valor verdadeiro mais um termo de erro aleatório:

\displaystyle X_{ij} = T_j + \epsilon_{ij}

Ou seja, para o sujeito j, o valor da variável X_i é igual ao escore verdadeiro T_j mais o erro aleatório daquele sujeito para aquele item, \epsilon_{ij}. O que o escore verdadeiro significa nessa situação é um ponto de debate espinhoso, e vamos simplesmente assumir que o modelo descreve bem o problema. Com base nessa definição, vamos explorar como podemos compreender a variância e covariância dos indicadores.

Comecemos com a variância:

\displaystyle \mathrm{var}[X_i] = \mathrm{var}[T + \epsilon_i]

Primeiro pressuposto adicional: para facilitar as contas, vamos assumir que T e \epsilon são independentes, e, portanto, \mathrm{cov}[T, \epsilon] = 0:

\displaystyle = \mathrm{var}[T] + \mathrm{var}[\epsilon_i]

Portanto, a variância dos itens é simplesmente a soma da variância do escore verdadeiro e a variância residual referente ao item i. Passemos à covariância entre um par de itens qualquer:

\displaystyle \mathrm{cov}[X_i, X_j] = \mathbb{E}[X_iX_j] - \mathbb{E}[X_i]\mathbb{E}[X_j] \\ = \mathbb{E}[(T + \epsilon_i)(T + \epsilon_j)] - \mathbb{E}[T + \epsilon_i]\mathbb{E}[T + \epsilon_j] \\ = \mathbb{E}[T^2 + T\epsilon_j + T\epsilon_i + \epsilon_i\epsilon_j] - (\mathbb{E}[T] + \mathbb{E}[\epsilon_i])(\mathbb{E}[T] + \mathbb{E}[\epsilon_j])

Aqui, novamente, retomamos o pressuposto de que T e \epsilon são independentes e acrescentamos mais dois pressupostos para facilitar as contas: primeiro, que o erro aleatório não possui nenhum viés e seu valor esperado é igual a zero para todos os itens, \mathbb{E}[\epsilon] = 0; e assumimos também que não existe nenhuma covariância entre erros de diferentes itens, \mathrm{cov}[\epsilon_i, \epsilon_j] = 0. Da posse desses pressupostos, fica fácil simplificar a expressão acima:

\displaystyle = \mathbb{E}[T^2] - \mathbb{E}[T]^2 \\ =\mathrm{var}[T]

Ou seja, a covariância entre qualquer par de indicadores, no modelo da Teoria Clássica dos Testes, é igual à variância do escore verdadeiro. Da posse dessas duas conclusões, podemos passar ao cálculo da variância das somas e a soma das variâncias dos itens de um instrumento sob o modelo da TCT. Comecemos com o numerador da razão que estamos explorando no coeficiente:

\displaystyle \sum_{i=1}^n \mathrm{var}[X_i] = \sum_{i=1}^n \mathrm{var}[T + \epsilon_i] \\ = \sum_{i=1}^n \mathrm{var}[T] + \mathrm{var}[\epsilon_i]

Vamos assumir que os itens da escala são paralelos, ou seja, que a parcela da variância do escore verdadeiro e a variância residual são iguais para todos os itens. Apesar de ser um pressuposto forte e pouco verossímil, ele facilitará os cálculos posteriores. Sob esse novo pressuposto:

\displaystyle = n\mathrm{var}[T] + n\mathrm{var}[\epsilon]

O resultado é óbvio sob o pressuposto de itens paralelos: a soma das variâncias é simplesmente o número de itens n vezes a variância dos itens – fica claro que a variância dos itens será necessariamente igual. Continuemos com o denominador, lembrando que a variância da soma de variáveis aleatórias é igual à soma das variâncias de cada variável mais duas vezes a soma da covariância de cada par.

\displaystyle \mathrm{var}\left[\sum_{i=1}^n X_i\right] = \sum_{i=1}^n \mathrm{var}[X_i] +2\sum_{i \neq j}\mathrm{cov}\left[ X_iX_j\right]

Como estamos supondo que os itens são paralelos e, portanto, sua covariância é igual à variância do escore verdadeiro, como demonstramos acima, basta multiplicarmos \mathrm{var}[T] pelo número de pares não redundantes. O número de pares é igual ao total de elementos para fora da diagonal da matriz de covariância, n(n-1), e os casos não redundantes é igual ao número de elementos acima ou abaixo da diagonal, \frac{n(n-1)}{2}.

\displaystyle=\sum_{i=1}^n \mathrm{var}[X_i] + 2\frac{n(n-1)}{2}\mathrm{var}[T] \\ = n\mathrm{var}[T] + n\mathrm{var}[\epsilon] + n(n-1)\mathrm{var}[T] \\ = n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]

Esse resultado aponta para uma propriedade interessante da soma dos indicadores: a parcela da variância da soma derivada da variância do escore verdadeiro aumenta quadraticamente em função do número de itens, enquanto que o variância atribuível ao erro aleatório aumenta apenas linearmente em função do número de itens. Isso significa que, à medida em que aumentamos o número de itens de um instrumento, aumentamos a proporção da variância da soma atribuível ao escore verdadeiro, tornando o escore total mais confiável.

Agora que reformulamos a variância da soma e a soma das variância em termos da variância do escore verdadeiro, da variância do erro aleatório e do número de itens, vamos retomar a equação que define o \alpha de Cronbach para chegarmos ao final de nossa intuição alinhando a interpretação do coeficiente com a noção de fidedignidade na TCT. Comecemos retomando a definição do \alpha, mas substituindo os termos da razão pelas equivalências encontradas acima.

\displaystyle \alpha = \frac{n}{n-1} \left( 1 - \frac{n\mathrm{var}[T] + n\mathrm{var}[\epsilon]}{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]} \right) \\ = \frac{n}{n-1} \left(\frac{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon] - n\mathrm{var}[T] - n\mathrm{var}[\epsilon]}{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]} \right)\\ = \frac{n}{n-1} \left(\frac{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon] - n\mathrm{var}[T] - n\mathrm{var}[\epsilon]}{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]} \right) \\ = \frac{n}{n-1} \left(\frac{(n-1)n\mathrm{var}[T]}{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]} \right) \\ = \frac{n^2\mathrm{var}[T]}{n^2\mathrm{var}[T] + n\mathrm{var}[\epsilon]} \\ = \frac{n^2\mathrm{var}[T]}{\mathrm{var}[Y]} = \rho_{xx}

Das equações acima, fica evidente que o \alpha expressa a proporção da variância da soma dos indicadores atribuível à variação dos escores verdadeiros. Essa propriedade do \alpha é exatamente a noção de reliability na Teoria Clássica dos Testes, geralmente indicada por \rho_{xx} para representar que ela pode ser hipoteticamente estimada a partir da correlação de dois testes paralelos. Como essa estimativa é obtida a partir de um único teste sob o conjunto de pressupostos expostos acima, ela é uma medida interna de fidedignidade.

Conclusões

  1. O \alpha é uma estimativa da fidedignidade (reliability) da soma (ou média) de um conjunto de indicadores de um escala, sob um conjunto de pressupostos a respeito do comportamento dos itens e definindo fidedignidade como proporção de variância explicada pelo escore verdadeiro;
  2. O \alpha estima a variância do escore verdadeiro a partir da razão da soma das variâncias e da variância das somas;
  3. A fidedignidade, como mensurada pelo \alpha, é um atributo do escore total de uma escala;
  4. O número de itens da escala influencia na estimativa de fidedignidade – uma escala problemática nesses termos pode ser melhorada aumentando o número de indicadores, conforme o gráfico abaixo, que apresenta o valor de \alpha em função do \mathrm{log}_{10}(n) e da razão entre a variância dos escores verdadeiros e a variância do erro para um único item.
alpha <- function(logN, varT=1, varE){
  (10^logN)^2*varT/
    ((10^logN)^2*varT + (10^logN)*varE)
}
curve(alpha(x, varE=4),
      from=0, to=3, col='skyblue',
      ylab=expression(alpha),
      xlab=expression(log[10](n)))
curve(alpha(x, varE=2), col='deepskyblue', add=T)
curve(alpha(x, varE=1), col='steelblue', add=T)
curve(alpha(x, varE=0.5), col='navyblue', add=T)
curve(alpha(x, varE=0.25), col='midnightblue', add=T)
legend(x=2, y=0.65, lty = 1, bty = 'n',
       legend=c(4, 2, 1, 1/2, 1/4),
       title=expression(paste(symbol(var)*tau %/% symbol(var)*epsilon)),
       col = c('midnightblue', 'navyblue', 'steelblue', 'deepskyblue', 'skyblue'))

Valor de α em função de log10(n) e var T/var e

Perigos dos rituais estatísticos

Perigos dos rituais estatísticos

O psicanalista Alfredo Jerusalinsky publicou um artigo no Zero Hora sobre a adoção de protocolos para para detecção precoce de riscos no desenvolvimento infantil pelo SUS por força de lei (PL 5501/2013). Apesar de não mencionado no artigo original, o protocolo sugerido ao longo da tramitação do projeto é o IRDI (Indicadores Clínicos de Risco para o Desenvolvimento Infantil), desenvolvido a partir de fundamentos psicanalíticos por uma ampla equipe, da qual Jerusalinsky fez parte.

A despeitos dos méritos (e problemas, como o risco de sobrediagnóstico) da criação de políticas para avaliação precoce do desenvolvimento, a proposta esbarra num problema caro ao Propagação de incerteza: o abuso de rituais estatísticos no lugar de um pensamento propriamente quantitativo. Nesse caso, o problema incide sobre a classe de publicações acadêmicas caracterizadas como “validação” de um instrumento psicológico.

Como a psicologia lida com diversos constructos que não são operacionalmente definidos, a interpretação desejada do escore de um teste exige que hipóteses sobre o constructo tenham sido derivadas da rede nomológica na qual ele se insere, e que essas hipóteses tenham sido corroboradas pelos dados obtidos. No lugar da derivação de hipóteses, todavia, alguns rituais foram instituídos como necessários para um estudo de validação considerado aceitável: normalmente, engloba a realização de análise fatorial exploratória dos indicadores, o cálculo de índices denominados de “fidedignidade” ou “consistência interna”, como o α de Cronbach; e correlações das mais variadas. O problema não é utilizar esses procedimentos em particular, mas o modo como são aplicados. Sem uma direção certa, definida pela cuidadosa proposição de hipóteses a partir das teorias que envolvem o constructo, qualquer resultado (preferencialmente, estatisticamente significativo) pode ser tomado como evidência que corrobora a validade da interpretação do instrumento.

Mas instrumentos para auxiliar no diagnóstico de Transtornos do Espectro Autista não sofrem do problema da ausência de uma definição operacional. Pelo contrário: apesar de podermos ter interesse em classificar diagnósticos em um contínuo, há um desfecho discreto bem definido — ou, pelo menos, definido de acordo com manuais diagnósticos.

Em um dos estudos publicados sobre o IRDI, porém, nenhum diagnóstico bem definido foi utilizado para avaliar a acurácia preditiva do instrumento. Em seu lugar, os autores preferiram comparar os resultados do IRDI, utilizando o corte arbitrário (ou, pelo menos, não justificado no artigo) da ausência de dois ou mais indicadores, com outro instrumento desenvolvido pela mesma equipe para identificar problemas de desenvolvimento ou “risco psíquico”. Como esse instrumento não possuia evidências claras para identificar o que se propunha, as análises todas acabam sendo um exercício de futilidade, justamente por não permitirem estabelecer relações claras entre o instrumento avaliado e as condições que ele deveria diagnosticar a partir do melhor critério disponível.

Mas o ritual foi seguido: foi feita uma análise fatorial (ou seria uma análise de componentes principais?), que não é relatada no texto ou sequer justificada; foram computados diversos riscos relativos e seus intervalos de confiança, ainda que o procedimento caia no erro comum de seleção a partir dos resultados, apresentando o desempenho de cada indicador individual que se mostrou mais “significativo”. Os resultados não contribuem para sabermos quão acurado é o IRDI para identificação precoce de TEA (ou outro desfecho melhor definido do que “risco psíquico”), mas o estudo foi publicado e pode agora servir de referência para a validade do instrumento.

Essa avaliação pode parecer excessivamente cáustica sobre o valor preditivo do IRDI. Afinal, trata-se de um único estudo criticado; a noção de TEA não estava definida em sua publicação; a amostra do estudo é bastante ampla e envolve diversas cidades brasileiras. Outros estudos sobre o instrumento, porém, não melhoram o quadro. Um estudo mais recente, por exemplo, propôs testar diretamente a capacidade do IRDI em identificar casos de TEA. A construção da amostra já traz dificuldades para a interpretação dos resultados: o IRDI se propõe avaliar indicadores para crianças entre zero e 18 meses de idade; mas, provavelmente pelas pressões concretas que afetam qualquer pesquisa, a amostra foi constituída de crianças entre três e sete anos. A dificuldade de interpretação é óbvia: como avaliar o poder preditivo de um instrumento se ele é aplicado de maneira retrospectiva, com base na memória dos cuidadores? Apesar da possibilidade do instrumento realmente indicar uma boa discriminação entre casos e controles, não é possível generalizar para a situação pretendida, ou seja, a avaliação precoce.

Nesse estudo, a validade do IRDI é avaliada, novamente, pelos rituais usuais: correlação com outra escala “já validada”, o CARS-BR (Childhood Autism Rating Scale); cálculo do α de Cronbach. O artigo possui o mérito de identificar a existência de um desfecho bem definido e avaliar o desempenho do instrumento em sua identificação. Porém, as estimativas de sensibilidade e especificidade são baseadas no desempenho da classificação com os próprios dados utilizados para treinamento. Como é conhecido no campo da Aprendizagem Estatística, as estimativas de erro baseadas exclusivamente nos dados da amostra superestimam a habilidade do classificador em diferenciar casos quando aplicado para dados fora da amostra — o que é conhecido como erro de generalização.

Outro problema é a utilização do CARS-BR como “padrão ouro” para avaliação de crianças com suspeita de TEA. A despeito do grande esforço no processo de adaptação cultural para o Brasil1, o estudo sugere ter validado o instrumento para língua portuguesa a partir dos já referidos rituais: correlações e o bendito α. Não apresenta, na proposta validação, se os pontos de cortes originais são condizentes com a nova realidade, se o instrumento é, de fato, capaz de distinguir casos. Ora, como dissemos, esses instrumentos possuem um critério operacional (razoavelmente) bem definido, e não precisam utilizar procedimentos usualmente aplicados para constructos sem definição operacional.

O que é necessário, no fim das contas, é que os estudos reconheçam que estão frente a um problema de decisão e a pesquisa precisa refletir isso. Primeiramente, mapear o espaço de decisão de maneira fundamentada — o que não parece difícil, p.e., se nos atermos a uma classificação binária, mas pode ser nuançada para levar em consideração a variabilidade dentro dos TEA. O espaço de decisão precisa estar relacionado com algum critério considerado definitivo, como o diagnóstico de TEA a partir das descrições de manuais diagnósticos na idade apropriada. Outro aspecto crucial é definir uma função de perda relevante para a decisão a ser tomada: uma opção natural é a perda do tipo 0-1, adaptada para dar pesos diferentes aos falsos positivos e falsos negativos, de acordo com o que for considerado mais grave. Por fim, a construção de um classificador também precisa ser feita de maneira cuidadosa: mesmo com apenas um preditor (p.e., o escore na escala), a utilização de polinômios pode permitir a identificação de limites de decisão mais flexíveis.

O problema dos rituais automatizados nas pesquisas de validação de instrumentos psicológicos não é monopólio da psicanálise, certamente. E não vejo problema com a inundação de artigos que reproduzem o automatismo metodológico, de valor (muitas vezes questionável) puramente acadêmico. Mas quando essa atitude automatizada se encontra com uma posição ambígua com relação à ciência e proposição de políticas públicas, os procedimentos deixam de ser puramente desleixo e passam a se tornar um problema de má-fé.


  1.  O foco na tradução e retrotradução de instrumentos como um dos principais mecanismos de adaptação cultural incorre, mais uma vez, nos problema dos automatismos de procedimentos. Afinal de contas, avaliar a invariância de mensuração de um instrumento entre grupos diferentes é uma questão a ser verificada empiricamente, e não pode ser garantida apenas por uma tradução fidedigna dos indicadores. Isso não significa, é claro, que o cuidado na tradução não deva ser observado; ele não deve, contudo, se sobrepor à avaliação quantitativa. 

Erro de medida e replicação

Andrew Gelman e Eric Loken publicaram, na revista Science, um pequeno artigo sobre o impacto do erro de medida sobre as inferências baseadas no uso de testes de hipótese nula. O argumento criticado por eles pode ser derivado a partir algumas premissas bem estabelecidas:

  1. A presença de ruído aleatório aditivo em uma variável – como o que pode ser introduzido pelo erro de medida – enviesa a estimativa das correlações dessa variável com outras para baixo. Em outras palavras, a presença de erro de medida subestima a correlação entre variáveis;
  2. O valor-p, principal estatística utilizada nos testes de hipótese nula, depende da magnitude do efeito, da variância do estimador e do tamanho da amostra. Quanto menor o efeito verdadeiro, sob a mesma condição de variância do estimador e tamanho da amostra, maior o valor-p.

A partir dessas duas premissas, a conclusão usual feita por pesquisadores é que, se uma estimativa atenuada pelo erro de medida atinge significância estatística, podemos ter confiança de que a efeito verdadeiro é ainda maior. O problema, argumentam Gelman e Loken, é que a maioria das pesquisa sofre seleção em função do valor-p encontrado.
Os valores-p também sofrem do que os autores chamam de Garden of Forking Paths: mudança nos procedimentos de análise aplicados com base em informações retiradas dos próprios dados. Essa flexibilidade nas análises, ou graus de liberdade do pesquisador, torna difícil a interpretação dos valores-p, pois os valores reportados não tem mais uma relação clara com o modelo original utilizado para sua computação.

Simulando alguns cenários

Quão grave é a situação? Para ter uma ideia, vamos rodar algumas simulações. Primeiro, consideremos o cenário usual de uma pesquisa correlacional nas ciências sociais:

  1. A correlação entre variáveis costuma ser relativamente pequena e, por isso, vamos considerar um coeficiente de regressão com magnitude 0.2.
  2. A variabilidade também costuma ser um problema em variáveis sociais, e, por isso, vamos considerar duas variáveis normalmente distribuídas com desvio-padrão igual a 4 para a variável independente, e 4.1 para a variável dependente.
  3. O erro da medida é modelado como ruído normal aditivo, iniciando com um cenário sem erro e indo até um cenário em que o erro corresponde a 50% da variância total.
  4. Consideramos também tamanho de amostra variável, indo de 10 a 1000, conforme o código abaixo.
  5. A “média sob seleção” aqui significa o valor esperado do estimador condicionado à ocorrência de significância estatística (p < 0.05).
  6. Para cada combinação de tamanho de amostra e nível de ruído, rodamos 2.000 simulações para obter os valores esperados e variância.
library(ggplot2)

simCor <- function(noise, n){
  X <- rnorm(n, 0, 4)
  # Gera Y como uma função linear de X
  Y <- 0.5 + 0.2*X + rnorm(n, 0, 4)
  # Acrescenta ruído à variável X
  X_noise <- X + rnorm(n, 0, noise)
  # Retorna as estimativas do coeficiente de X sob ruído
  summary(lm(Y~X_noise))$coefficients[2, c(1, 2, 4)]
}

noise <- seq(0, 4, length.out=5)
sampleSize <- c(10, 50, 100, 500, 1000)
out <- array(NA, c(length(sampleSize), length(noise), 3))
for (n in 1:length(sampleSize)){
  for (nLevel in 1:length(noise)){
    teste <- t(replicate(2000, simCor(noise[nLevel], sampleSize[n])))
    out[n, nLevel, 1] <- mean(teste[, 1])
    out[n, nLevel, 2] <- mean(teste[teste[, 3] < 0.05, 1])
    out[n, nLevel, 3] <- var(teste[, 1])
  }
}
sim <- data.frame(Amostra=as.factor(rep(sampleSize, length(noise))),
                  Ruído= 1 - 16/(16 + as.vector(t(replicate(length(sampleSize), noise)))^2),
                  Média=as.vector(out[, , 1]),
                  Variância=as.vector(out[, , 3]),
                  Média.Sob.Seleção=as.vector(out[, , 2]))
sim <- tidyr::gather(sim, Estatística, Valor, -Amostra, -Ruído)
ggplot(sim, aes(x=Ruído, y=Valor, color=Amostra)) +
  geom_point() +
  geom_line() +
  geom_hline(yintercept = 0.2, linetype=2) +
  facet_wrap(~Estatística)

Variação da estimativa de efeito em função da amostra e nível de ruído

A simulação revela algumas coisas interessantes, que não são óbvias na leitura do artigo. O primeiro painel apresenta o valor esperado da estimativa do efeito em função do nível de ruído na variável independente. A linha horizontal indica o valor verdadeiro do parâmetro. Como esperado, a presença de ruído na variável atenua o valor esperado da estimativa. O valor da atenuação é diretamente proporcional à proporção de ruído na variável medida – a relação linear é evidente.

Quando a estimativa é baseada na seleção em função da significância estatística, caso ilustrado no segundo quadro, o problema identificado no artigo de Gelman e Loken se mostra evidente. A não ser em casos com amostras consideráveis (500 e 1000 casos), a estimativa de efeito quase sempre superestima o valor do parâmetro. O caso fica ainda mais grave nas amostras pequenas, mas o efeito do ruído continua sendo de atenuar o valor médio das estimativas.

Por fim, o último quadro apresenta a variância do estimador em função do ruído. Curiosamente, a presença de ruído diminui a variância em vez de aumentá-la. O motivo, porém, é simples: com uma magnitude de efeito consideravelmente pequena, o ruído puxa a estimativa cada vez mais para próximo de zero. A variância do estimador também reduz drasticamente com o aumento da amostra.

Conclusões…

O grande vilão da história toda é, no fim das contas, a seleção com base na significância estatística. A presença de ruído enviesa as estimativas, é verdade, assim como o tamanho da amostra dá lugar a um aumento na variância, fenômenos bem conhecidos. Mas não é nem o ruído ou o tamanho da amostra a maior fonte de viés, e sim a insistência de reportar efeito apenas com base em sua significância estatística.

Um fetiche pelo N

Um fetiche pelo N

A large amount of data is in no way synonymous with a large amount of information. In some settings at least, if a modest amount of poor quality data is likely to be modestly misleading, an extremely large amount of poor quality data may be extremely misleading (Cox & Donnely, 2011).

Em qualquer curso introdutório à estatística ou qualquer livro didático sobre o tema, somos forçosamente apresentados a diversos procedimentos de amostragem. A motivação é óbvia: se boa parte da estatística tradicional se ocupa da inferência de valores de uma população bem definida a partir de uma amostra, é crucial definir como obter dados que tornem a inferência confiável.

Em minha experiência, porém, esse material, apesar de compor um dos tópicos obrigatórios da introdução à estatística na formação de pesquisadores em psicologia, uma vez ministrado, é esquecido sem graves prejuízos. Por que relembrar? Afinal, queremos produzir ciência, e não duplicar o trabalho que o IBGE já realiza. E, para produzir ciência, temos teorias e hipóteses a serem testadas, e apenas um interesse marginal em características de populações concretas.

De fato, uma pergunta de pesquisa em psicologia não versará necessariamente sobre uma população concretamente definida. Nossas pesquisas abordam construtos teóricos, sua estrutura e relações. Como nos interesa estabelecer relações de dependência entre determinadas variáveis, a preocupação do pesquisador costuma passar diretamente para os procedimentos estatísticos de análise que serão utilizados para avaliar se a correlação predita está lá, de fato.

O problema, porém, é que esquecemos que a validade das inferências sobre a relação entre essas variáveis depende dos procedimentos de amostragem. Em outras palavras, além do processo gerador de dados pressuposto pela teoria, temos mais uma camada de mecanismos que levam à composição dos dados analisados. Se esse mecanismo não é considerado – ou, como acontece, é substituído por um mecanismo muito mais simples, como a amostragem aleatória simples – quaisquer inferências serão puramente imaginárias.

A preocupação em estabelecer correlações entre construtos teóricos e o conveniente esquecimento das aulas sobre amostragem redundam na utilização de amostras de conveniência. Alunos do curso de psicologia, universitários imediatamente disponíveis, compartilhamento de um formulário pela internet… A não ser no caso de pesquisas experimentais, esses procedimentos todos retiram do pesquisador a capacidade de definir um plano amostral claro e obscurecem o tipo de inferência que é possível de ser feita.

Uma amostra de conveniência pode convir muito bem para um estudo piloto sobre uma questão ainda pouco investigada. O problema maior com essas amostras é quando o pesquisador acredita poder aplicar procedimentos estatísticos que foram desenvolvidos especificamente para casos idealizados, como amostragem ou atribuição aleatória simples, para poder concluir sobre uma mal-definida significância da relação encontrada entre variáveis.

A falta de uma boa compreensão do papel do plano amostral nas inferências que um estudo pode fazer é substituída por um fetiche: o tamanho da amostra. A questão do tamanho amostral, ou “N” para os íntimos, é uma preocupação constante na realização e avaliação de pesquisas psicológicas, como uma espécie de garantia intuitiva a respeito da validade dos achados. O problema, todavia, como expressado na epígrafe, é que uma grande quantidade de dados de pouca qualidade podem ser extremamente enganoso – ainda mais quando acreditamos que o tamanho da amostra serve de alguma garantia contra erros.

Sem dúvida, o tamanho da amostra nos permite realizar inferências mais precisas, desde que a análise leve em consideração possíveis vieses nos procedimentos amostrais. Se, como de costume, a análise estatística assume uma forma de amostragem idealizada e completamente inadequada para como os dados foram coletados, não só a estimativa de incerteza sobre os resultados será muito menor do que é de fato, mas as próprias estimativas das quantidades de interesse poderão estar enviesadas. E pior: não saberemos quão grave ou em que direção o viés afeta as conclusões.

Como os modelos estatísticos idealizam a realidade, alguma discrepância sempre é esperada para além do erro teórico predito. Porém, se pegarmos casos simples como a predição da proporção de votos em uma eleição, e verificarmos o que ocorreu nas últimas eleições presidenciais, tanto no Brasil quantos nos EUA, o erro causa assombro. Se não podemos confiar em um modelo relativamente simples para uma informação consideravelmente bem delimitada, por que aceitaríamos procedimentos ainda inferiores para a construção de conhecimento?

Isso não é um problema enquanto o conhecimento fica restrito à academia. Afinal, a quem interessa os resultados de pesquisas senão aos pares imediatos do pesquisador, que poderão ter uma nova referência para citar? Com teorias fracas e pequenos fatos estilizados construídos a partir de uma empiria meio caolha, a maior parte das pesquisas em psicologia é pouco informativa sobre o conhecimento que pretende construir e ainda menos interessante para os psicólogos que estão praticando.

Contudo, esse fetiche pelo N, substituto de um bom entendimento do papel da amostragem sobre as inferências estatísticas, acaba chegando à prática psicológica pelos mais variados meios. Um deles, talvez dos mais perniciosos, é na validação de testes psicológicos. Basta folhear alguns manuais de testes, mesmo aqueles com parecer favorável do SATEPSI, para verificar que as propriedades psicométricas e, ainda mais grave, as normas são baseadas, via de regra, em amostras de conveniência. O crivo de avaliação de testes submetidos ao SATEPSI, ingrata surpresa, inicia as considerações sobre as amostras questionando sobre o tal do “N”. As únicas considerações seguintes são:

  1. Há “controle” de variáveis importantes?
  2. O tamanho da amostra é suficiente?

Ah, se nossas teorias fossem capaz de elencar todas as variáveis intervenientes e moderadoras do desfecho de interesse… Infelizmente, não é o caso, e toda evidência aponta para o contrário: o crud factor proposto por Meehl é justamente a constatação de que, nas ciências sociais e humanas, parece haver correlações “de fundo” entre variáveis aparentemente não relacionadas. Portanto, ainda não conseguimos “controlar por variáveis relevantes”, e provavelmente não conseguiremos em um futuro próximo, o que torna procedimentos de randomização necessários.

O caso é ainda mais grave no estabelecimento de normas para os testes. As normas de um teste são cruciais para estabelecer o lugar de um sujeito avaliado dentro de um contínuo que a escala se propõe a medir. Mas esse contínuo só é válido se, como se diz, estamos comparando laranjas com laranjas. Podemos de fato conceber que o escore de QI total do WISC-III, por exemplo, obtido por uma criança proveniente de uma família de baixa renda do Norte do país pode ser diretamente comparado à amostra de referência, obtida em Pelotas-RS? Não é preciso refletir muito para concluir que a amostra utilizada para estimar as quantidades de referência não é representativa do Brasil, mas qualquer psicólogo que utilizar o WISC-III para avaliar a inteligência de crianças irá comparar, invariavelmente, com informações que talvez sejam válidas para uma cidade.

Em última instância, amostras fortuitas introduzem vieses fortuitos; quanto maior elas forem, maiores os vieses possíveis e menor a segurança que podemos ter na validade das conclusões.