Quelle est la relation entre la probabilité de profil et les intervalles de confiance?

18

Pour faire ce graphique, j'ai généré des échantillons aléatoires de taille différente à partir d'une distribution normale avec moyenne = 0 et sd = 1. Les intervalles de confiance ont ensuite été calculés en utilisant des seuils alpha allant de .001 à .999 (ligne rouge) avec la fonction t.test (), la probabilité de profil a été calculée en utilisant le code ci-dessous que j'ai trouvé dans les notes de cours mises en ligne (je peux ' t trouver le lien pour le moment Edit: Found it ), cela est indiqué par les lignes bleues. Les lignes vertes indiquent la densité normalisée à l'aide de la fonction R densité () et les données sont affichées par les boîtes à moustaches au bas de chaque graphique. Sur la droite se trouve un tracé de chenille des intervalles de confiance à 95% (rouge) et 1 / 20ème des intervalles de vraisemblance max (bleu).

R Code utilisé pour la vraisemblance du profil:

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

entrez la description de l'image ici

Ma question spécifique est de savoir s'il existe une relation connue entre ces deux types d'intervalles et pourquoi l'intervalle de confiance semble être plus conservateur pour tous les cas, sauf lorsque n = 3. Des commentaires / réponses sur la validité de mes calculs (et une meilleure façon de le faire) et la relation générale entre ces deux types d'intervalles sont également souhaités.

Code 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)
Ballon
la source
Dans vos notes de cours, mnc'est une faute de frappe mu, et non mean(dat). Comme je vous l'ai dit dans les commentaires de votre autre question , cela devrait être clair à partir de la page 23 des définitions.
Elvis
@Elvis Je ne pense pas. mn est défini à la page 18 des notes.
Flask
J'ai essayé de clarifier le concept de probabilité de profil. Pouvez-vous commenter un peu plus ce que vous faites dans le code ci-dessus?
Elvis
3
@Elvis Je ne comprends pas non plus. Un intervalle de confiance basé sur la probabilité du profil devrait être construit à l'aide des centiles, qui n'apparaissent nulle part. χ2
Stéphane Laurent
1
@ StéphaneLaurent Je ne suis pas sûr que le code d' origine est fournirai des intervalles de confiance. Plutôt 1/20 maximum d'intervalles de vraisemblance. Je crois que le nom des intervalles de confiance dans mon graphique est des intervalles de confiance "de type wald" et les lignes rouges sur les graphiques sont des "courbes de confiance" décrites sur cette page wikipedia
Flask

Réponses:

10

Je ne donnerai pas une réponse complète (j'ai du mal à comprendre ce que vous faites exactement), mais j'essaierai de clarifier comment la probabilité de profil est construite. Je pourrai terminer ma réponse plus tard.

La probabilité totale pour un échantillon normal de taille est L ( μ , σ 2 ) = ( σ 2 ) - n / 2 exp ( - Σ i ( x i - μ ) 2 / 2 σ 2 ) .n

L(μ,σ2)=(σ2)n/2exp(i(xiμ)2/2σ2).

μσ2μ

LP(μ)=L(μ,σ2^(μ))
σ2^(μ)μ
σ2^(μ)=argmaxσ2L(μ,σ2).

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

LP(μ)=(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")

probabilité de profil

Lien avec la probabilité Je vais essayer de mettre en évidence le lien avec la probabilité avec le graphique suivant.

Définissez d'abord la probabilité:

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

Faites ensuite un tracé de contour:

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)

tracé de contour de L

Les valeurs de la vraisemblance du profil sont les valeurs prises par la vraisemblance le long de la parabole rouge.

Vous pouvez utiliser la vraisemblance de profil comme une vraisemblance classique univariée (cf. la réponse de @ Prokofiev). Par exemple, le MLEμ^

σ2^(μ)

Vous pouvez également utiliser la probabilité de profil pour créer des tests de score, par exemple.

Elvis
la source
mu dans le code est une séquence de valeurs de faible à élevée, la probabilité à chacune de ces valeurs est divisée par la probabilité à la moyenne de l'échantillon (mn). Il s'agit donc d'une probabilité normalisée.
Flask
Je pense que c'est la même chose mais pas normalisé. Pouvez-vous le mettre en code R ou tracer autrement la fonction pour certaines données afin que nous puissions comparer?
Flask
C'est ici. Au début, je pensais que mnc'était une faute de frappe, maintenant je pense que le code R est tout faux. Je vérifierai demain - il est tard si je vis.
Elvis
Vous pouvez avoir raison. Je ne comprends pas comment le code parvient à le normaliser. Oh, je comprends, la "normalisation" divise simplement par le maximum?
Elvis
1
Je pense que c'est pour le rendre facile à voir lorsque le rapport de vraisemblance est inférieur à un certain seuil (par exemple 1 / 20e max) à une hypothèse nulle (par exemple zéro).
Flask
7

χk2

0.14795%

Ce sont des résultats classiques et je vais donc simplement fournir quelques références à ce sujet:

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

Le code R suivant montre que, même pour de petits échantillons, les intervalles obtenus avec les deux approches sont similaires (je réutilise l'exemple d'Elvis):

Notez que vous devez utiliser la probabilité de profil normalisée.

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

Si nous utilisons un échantillon plus grand, les intervalles de confiance sont encore plus proches:

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

UN POINT IMPORTANT:

Notez que pour des échantillons spécifiques, différents types d'intervalles de confiance peuvent différer en termes de longueur ou d'emplacement, ce qui compte vraiment, c'est leur couverture. À long terme, ils devraient tous offrir la même couverture, indépendamment de leur différence pour des échantillons spécifiques.

Prokofiev
la source
@Prokoflev s'il existe une relation simple entre les intervalles de confiance calculés avec la fonction R t.test () et ceux calculés par le code de fonction de vraisemblance ci-dessus, pouvez-vous le poster. Je suis particulièrement intéressé par le cas n = 3. Malheureusement, j'ai peu de connaissances en mathématiques, donc de nombreux articles me conduisent dans le terrier du lapin en recherchant les noms des symboles et ce qu'ils représentent, etc. quand quelques lignes de code (le plus simple est R) pourraient me l'expliquer.
Flask
@Flask Êtes-vous intéressé à obtenir des intervalles de confiance pour les paramètres d'une distribution normale ou un cadre plus général?
Prokofiev
@Prokoflev spécifiquement pour la moyenne d'une distribution normale comme le montre mon exemple dans la question. Je me demande surtout pourquoi les intervalles de confiance sont plus conservateurs sauf dans le cas n = 3.
Flask
95%
1
Je commence à croire que je devrais multiplier les intervalles de vraisemblance par un quantile de la distribution normale ou chisquare pour obtenir l'intervalle de confiance correspondant.
Flask
1

χ2normuneljeze

  1. La log-vraisemblance du profil est approximativement quadratique
  2. Il existe une transformation de paramètre qui rend la log-vraisemblance du profil approximativement quadratique.

Le quadratique est important car il définit une distribution normale à l'échelle logarithmique. Plus elle est quadratique, meilleure est l'approximation et les IC résultants. Votre choix de 1 / 20e de coupure pour les intervalles de vraisemblance équivaut à plus d'un IC à 95% dans la limite asymptotique, c'est pourquoi les intervalles bleus sont généralement plus longs que les rouges.

Maintenant, il y a un autre problème avec la probabilité de profil qui nécessite une certaine attention. Si vous avez de nombreuses variables sur lesquelles vous profilez, alors si le nombre de points de données par dimension est faible, la probabilité de profil peut être très biaisée et optimiste. Les probabilités de profil marginales, conditionnelles et modifiées sont ensuite utilisées pour réduire ce biais.

Donc, la réponse à votre question est OUI ... le lien est la normalité asymptotique de la plupart des estimateurs du maximum de vraisemblance, comme en témoigne la distribution chi carré du rapport de vraisemblance.


la source
" Si vous avez beaucoup de variables sur lesquelles vous profilez, alors si le nombre de points de données par dimension est faible, la probabilité de profil peut être très biaisée et optimiste " Optimiste par rapport à quoi?
Flask
@Flask Par optimiste, je veux dire qu'il sera trop étroit pour fournir la probabilité de couverture nominale lors du traitement comme un intervalle de confiance.
Je vois, merci, mais dans mon cas spécifique, c'est en fait pessimiste? Je suis confus sur ce point quant à savoir si nous parlons des intervalles de vraisemblance ou des intervalles de confiance dérivés des probabilités.
Flask
@Flask Je pense que vos intervalles semblent pessimistes parce que vous comparez 1/20 ème intervalle de vraisemblance (probabilité relative de 5%) avec un IC à 95%. Comme indiqué par d'autres ici, vous voudriez vraiment le comparer à un intervalle de probabilité relative de 15% d'avoir des pommes aux pommes ... au moins asymptotiquement. Votre intervalle de vraisemblance tel qu'il se présente envisage plus d'options comme plausbile.
J'ai détaillé le problème réel auquel je souhaite appliquer ce que j'apprends ici . Je crains que dans le cas où la distribution d'échantillonnage soit inconnue (mais pas probablement pas normale) et complexe que vos deux exigences puissent ne pas tenir. Pourtant, les probabilités de profil que j'ai calculées semblent être normales et raisonnables. Est-ce que la distribution d'échantillonnage de la moyenne devrait être normalement distribuée?
Flask