Qual é a relação entre a probabilidade do perfil e os intervalos de confiança?

18

Para fazer esse gráfico, gerei amostras aleatórias de tamanho diferente de uma distribuição normal com média = 0 e sd = 1. Os intervalos de confiança foram então calculados usando pontos de corte alfa variando de 0,001 a 0,999 (linha vermelha) com a função t.test (), a probabilidade do perfil foi calculada usando o código abaixo do que encontrei nas notas de aula colocadas on-line (eu posso ' t encontre o link no momento Edit: Found ), isso é mostrado pelas linhas azuis. As linhas verdes mostram a densidade normalizada usando a função R density () e os dados são mostrados pelos gráficos de caixa na parte inferior de cada gráfico. À direita, há um gráfico de lagarta dos intervalos de confiança de 95% (vermelho) e 1/20 dos intervalos de probabilidade máxima (azul).

Código R usado para a probabilidade do perfil:

  #mn=mean(dat)
  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )

insira a descrição da imagem aqui

Minha pergunta específica é se existe uma relação conhecida entre esses dois tipos de intervalos e por que o intervalo de confiança parece ser mais conservador para todos os casos, exceto quando n = 3. Comentários / respostas sobre se meus cálculos são válidos (e uma maneira melhor de fazer isso) e a relação geral entre esses dois tipos de intervalos também são desejados.

Código R:

samp.size=c(3,4,5,10,20,1000)
cnt2<-1
ints=matrix(nrow=length(samp.size),ncol=4)
layout(matrix(c(1,2,7,3,4,7,5,6,7),nrow=3,ncol=3, byrow=T))
par(mar=c(5.1,4.1,4.1,4.1))
for(j in samp.size){


  #set.seed(200)
  dat<-rnorm(j,0,1)
  vals<-seq(.001,.999, by=.001)
  cis<-matrix(nrow=length(vals),ncol=3)
  cnt<-1
  for(ci in vals){
    x<-t.test(dat,conf.level=ci)$conf.int[1:2]
    cis[cnt,]<-cbind(ci,x[1],x[2])
    cnt<-cnt+1
  }


  mn=mean(dat)
  n=length(dat)
  high<-max(c(dat,cis[970,3]), na.rm=T)
  low<-min(c(dat,cis[970,2]), na.rm=T)
  #high<-max(abs(c(dat,cis[970,2],cis[970,3])), na.rm=T)
  #low<--high


  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )


  plot(muVals, likVals, type = "l", lwd=3, col="Blue", xlim=c(low,high),
       ylim=c(-.1,1), ylab="Likelihood/Alpha", xlab="Values",
       main=c(paste("n=",n), 
              "True Mean=0 True sd=1", 
              paste("Sample Mean=", round(mn,2), "Sample sd=", round(sd(dat),2)))
  )
  axis(side=4,at=seq(0,1,length=6),
       labels=round(seq(0,max(density(dat)$y),length=6),2))
  mtext(4, text="Density", line=2.2,cex=.8)

  lines(density(dat)$x,density(dat)$y/max(density(dat)$y), lwd=2, col="Green")
  lines(range(muVals[likVals>1/20]), c(1/20,1/20), col="Blue", lwd=4)
  lines(cis[,2],1-cis[,1], lwd=3, col="Red")
  lines(cis[,3],1-cis[,1], lwd=3, col="Red")
  lines(cis[which(round(cis[,1],3)==.95),2:3],rep(.05,2), 
        lty=3, lwd=4, col="Red")
  abline(v=mn, lty=2, lwd=2)
  #abline(h=.05, lty=3, lwd=4, col="Red")
  abline(h=0, lty=1, lwd=3)
  abline(v=0, lty=3, lwd=1)

  boxplot(dat,at=-.1,add=T, horizontal=T, boxwex=.1, col="Green")
  stripchart(dat,at=-.1,add=T, pch=16, cex=1.1)

  legend("topleft", legend=c("Likelihood"," Confidence Interval", "Sample Density"),
         col=c("Blue","Red", "Green"), lwd=3,bty="n")

  ints[cnt2,]<-cbind(range(muVals[likVals>1/20])[1],range(muVals[likVals>1/20])[2],
                     cis[which(round(cis[,1],3)==.95),2],cis[which(round(cis[,1],3)==.95),3])
  cnt2<-cnt2+1
}
par(mar=c(5.1,4.1,4.1,2.1))


plot(0,0, type="n", ylim=c(1,nrow(ints)+.5), xlim=c(min(ints),max(ints)), 
     yaxt="n", ylab="Sample Size", xlab="Values")
for(i in 1:nrow(ints)){
  segments(ints[i,1],i+.2,ints[i,2],i+.2, lwd=3, col="Blue")
  segments(ints[i,3],i+.3,ints[i,4],i+.3, lwd=3, col="Red")
}
axis(side=2, at=seq(1.25,nrow(ints)+.25,by=1), samp.size)
Frasco
fonte
Em suas notas de aula, mné um erro de digitação mu, e não mean(dat). Como falei nos comentários de sua outra pergunta , isso deve ficar claro na página 23 das definições.
Elvis
@ Elvis, acho que não. mn é definido na página 18 das notas.
Flask
Tentei esclarecer o conceito de probabilidade de perfil. Você pode comentar um pouco mais sobre o que está fazendo no código acima?
Elvis
3
@ Elvis Nem eu entendo. Um intervalo de confiança baseado na probabilidade do perfil deve ser construído com a ajuda dos percentis , que não aparecem em lugar algum. χ2
Stéphane Laurent
1
@ StéphaneLaurent Não sei se o código original está fornecendo intervalos de confiança. Em vez de 1/20 de intervalos de probabilidade máximos. Acredito que o nome para os intervalos de confiança na minha plotagem sejam intervalos de confiança do tipo "wald" e as linhas vermelhas nas plotagens sejam "curvas de confiança" descritas nesta página da wikipedia
Flask

Respostas:

10

Não darei uma resposta completa (tenho dificuldade em entender exatamente o que você está fazendo), mas tentarei esclarecer como a probabilidade de perfil é criada. Posso completar minha resposta mais tarde.

A probabilidade total de uma amostra normal de tamanho é G ( μ , σ 2 ) = ( σ 2 ) - n / 2 exp ( - Σ i ( x i - μ ) 2 / 2 σ 2 ) .n

eu(μ,σ2)=(σ2)-n/2exp(-Eu(xEu-μ)2/2σ2).

μσ2μ

euP(μ)=eu(μ,σ2^(μ))
σ2^(μ)μ
σ2^(μ)=argmaxσ2eu(μ,σ2).

σ2^(μ)=1nk(xk-μ)2.

euP(μ)=(1nk(xk-μ)2)-n/2exp(-n/2).

exp(-n/2)

> data(sleep)
> difference <- sleep$extra[11:20]-sleep$extra[1:10]
> Lp <- function(mu, x) {n <- length(x); mean( (x-mu)**2 )**(-n/2) }
> mu <- seq(0,3, length=501)
> plot(mu, sapply(mu, Lp, x = difference), type="l")

probabilidade de perfil

Link com a probabilidade Tentarei destacar o link com a probabilidade no gráfico a seguir.

Primeiro defina a probabilidade:

L <- function(mu,s2,x) {n <- length(x); s2**(-n/2)*exp( -sum((x-mu)**2)/2/s2 )}

Em seguida, faça um gráfico de contorno:

sigma <- seq(0.5,4, length=501)
mu <- seq(0,3, length=501)

z <- matrix( nrow=length(mu), ncol=length(sigma))
for(i in 1:length(mu))
  for(j in 1:length(sigma))
    z[i,j] <- L(mu[i], sigma[j], difference)

# shorter version
# z <- outer(mu, sigma, Vectorize(function(a,b) L(a,b,difference)))

contour(mu, sigma, z, levels=c(1e-10,1e-6,2e-5,1e-4,2e-4,4e-4,6e-4,8e-4,1e-3,1.2e-3,1.4e-3))

σ2^(μ)

hats2mu <- sapply(mu, function(mu0) mean( (difference-mu0)**2 ))
lines(mu, hats2mu, col="red", lwd=2)

plotagem de contorno de L

Os valores da probabilidade do perfil são os valores obtidos pela probabilidade ao longo da parábola vermelha.

μ^

σ2^(μ)

Você também pode usar a probabilidade do perfil para criar testes de pontuação, por exemplo.

Elvis
fonte
mu no código é uma sequência de valores de baixo a alto, a probabilidade em cada um desses valores é dividida pela probabilidade na média da amostra (mn). Portanto, é uma probabilidade normalizada.
Flask
Eu acho que isso é a mesma coisa, mas não normalizado. Você pode colocá-lo no código R ou plotar a função para alguns dados para que possamos comparar?
Flask
Aqui está. No começo eu pensei que mnera um erro de digitação, agora acho que o código R está errado. Vou checar amanhã - é tarde onde eu moro.
Elvis
Tu podes estar certo. Não entendo como o código consegue normalizá-lo. Ah, entendi, a "normalização" está apenas dividindo pelo máximo?
Elvis
1
Eu acho que é fácil ver quando a razão de verossimilhança é menor que algum limite (por exemplo, 1/20 no máximo) em alguma hipótese nula (por exemplo, zero).
Flask
7

χk2

0,14795%

Estes são resultados clássicos e, portanto, simplesmente fornecerei algumas referências sobre isso:

http://www.jstor.org/stable/2347496

http://www.stata-journal.com/sjpdf.html?articlenum=st0132

http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture11.htm

http://en.wikipedia.org/wiki/Likelihood-ratio_test

http://en.wikipedia.org/wiki/Likelihood_function#Profile_likelihood

O código R a seguir mostra que, mesmo para amostras pequenas, os intervalos obtidos com as duas abordagens são semelhantes (estou reutilizando o exemplo de Elvis):

Observe que você precisa usar a probabilidade de perfil normalizada.

data(sleep)
x <- sleep$extra[11:20]-sleep$extra[1:10]
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(0,3, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(0.5,1.5))$root,uniroot(Rpt,c(1.51,3))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

Se usarmos um tamanho de amostra maior, os intervalos de confiança serão ainda mais próximos:

set.seed(123)
x <- rnorm(100)
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(-0.5,0.5, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(-0.4,0))$root,uniroot(Rpt,c(0,0.4))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

UM PONTO IMPORTANTE:

Observe que, para amostras específicas, diferentes tipos de intervalos de confiança podem diferir em termos de comprimento ou localização, o que realmente importa é a cobertura. A longo prazo, todos eles devem fornecer a mesma cobertura, independentemente de quanto diferem para amostras específicas.

Prokofiev
fonte
@Prokoflev se houver alguma relação simples entre os intervalos de confiança calculados com a função R t.test () e os calculados pelo código da função de probabilidade acima, você pode publicá-lo. Estou especialmente interessado no caso n = 3. Infelizmente, como tenho pouco conhecimento em matemática, muitos trabalhos me levam até a toca do coelho, procurando os nomes dos símbolos e o que eles representam etc., quando algumas linhas de código (o mais fácil é R) poderiam me explicar.
Flask
@Flask Você está interessado em obter intervalos de confiança para os parâmetros de uma distribuição normal ou uma estrutura mais geral?
Prokofiev
@Prokoflev especificamente para a média de uma distribuição normal, como mostrado no meu exemplo na pergunta. Estou especialmente me perguntando por que os intervalos de confiança são mais conservadores, exceto no caso n = 3.
Flask
95%
1
Estou começando a acreditar que eu deveria estar multiplicando os intervalos de probabilidade por algum quantil da distribuição normal ou do qui-quadrado para obter o intervalo de confiança correspondente.
Flask
1

χ2normumaeuEuzed

  1. A probabilidade de log do perfil é aproximadamente quadrática
  2. Existe uma transformação de parâmetro que torna a probabilidade do log do perfil aproximadamente quadrática.

O quadrático é importante porque define uma distribuição normal em escala logarítmica. Quanto mais quadrático, melhor a aproximação e os ICs resultantes. Sua escolha do ponto de corte de 1/20 para os intervalos de probabilidade é equivalente a mais de um IC de 95% no limite assintótico, motivo pelo qual os intervalos azuis geralmente são mais longos que os vermelhos.

Agora, há outro problema com a probabilidade de perfil que precisa de alguma atenção. Se você tiver várias variáveis ​​nas quais está criando um perfil, se o número de pontos de dados por dimensão for baixo, a probabilidade do perfil poderá ser muito tendenciosa e otimista. As probabilidades de perfil marginal, condicional e modificado são usadas para reduzir esse viés.

Portanto, a resposta para sua pergunta é SIM ... a conexão é a normalidade assintótica da maioria dos estimadores de probabilidade máxima, como se manifesta na distribuição qui-quadrado da razão de verossimilhança.


fonte
" Se você possui muitas variáveis ​​em que está criando um perfil, se o número de pontos de dados por dimensão for baixo, a probabilidade do perfil pode ser muito tendenciosa e otimista " Otimista em comparação com o que?
Flask
@Flask Por otimista, quero dizer que será muito estreito para fornecer a probabilidade de cobertura nominal ao tratá-lo como um intervalo de confiança.
Entendo, obrigado, mas no meu caso específico é realmente pessimista? Estou confuso quanto a esse ponto sobre se estamos falando sobre os intervalos de probabilidade ou de confiança derivados das probabilidades.
Flask
@ Atlas Eu acho que os intervalos parecem pessimistas porque você está comparando o 1/20 de intervalo de probabilidade (5% de probabilidade relativa) com um IC de 95%. Como afirmado por outros aqui, você realmente gostaria de compará-lo a um intervalo de probabilidade relativa de 15% para ter maçãs com maçãs ... pelo menos assintoticamente. Seu intervalo de probabilidade, tal como está, está considerando mais opções como plausíveis.
Eu detalhei o problema real ao qual desejo aplicar o que estou aprendendo aqui . Eu me preocupo que, no caso em que a distribuição da amostra seja desconhecida (mas provavelmente não seja normal) e complexa, seus dois requisitos possam não se sustentar. No entanto, as probabilidades de perfil que calculei parecem normais e razoáveis. Será que a distribuição amostral da média deve ser normalmente distribuída?
Flask