Apa hubungan antara kemungkinan profil dan interval kepercayaan?

18

Untuk membuat bagan ini, saya membuat sampel acak dengan ukuran berbeda dari distribusi normal dengan mean = 0 dan sd = 1. Interval kepercayaan kemudian dihitung menggunakan cutoff alfa mulai dari 0,001 hingga 0,999 (garis merah) dengan fungsi t.test (), kemungkinan profil dihitung menggunakan kode di bawah ini yang saya temukan dalam catatan kuliah yang dimasukkan ke baris (saya bisa ' t menemukan tautan saat ini Edit: Ditemukan ), ini ditunjukkan oleh garis biru. Garis hijau menunjukkan kerapatan yang dinormalisasi menggunakan fungsi R kerapatan () dan data ditampilkan oleh plot kotak di bagian bawah setiap bagan. Di sebelah kanan adalah plot ulat dari interval kepercayaan 95% (merah) dan 1/20 dari interval kemungkinan maksimal (biru).

Kode R yang digunakan untuk kemungkinan 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)
                    }
  )

masukkan deskripsi gambar di sini

Pertanyaan khusus saya adalah apakah ada hubungan yang diketahui antara kedua jenis interval ini dan mengapa interval kepercayaan tampaknya lebih konservatif untuk semua kasus kecuali ketika n = 3. Komentar / jawaban tentang apakah perhitungan saya valid (dan cara yang lebih baik untuk melakukan ini) dan hubungan umum antara kedua jenis interval ini juga diinginkan.

Kode 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)
Labu
sumber
Dalam catatan kuliah Anda, mnadalah salah ketik untuk mu, dan tidak mean(dat). Seperti yang saya katakan di komentar untuk pertanyaan Anda yang lain , ini harus jelas dari definisi halaman 23.
Elvis
@ Elvis saya tidak berpikir begitu. mn didefinisikan pada halaman 18 catatan.
Labu
Saya mencoba mengklarifikasi konsep kemungkinan profil. Bisakah Anda berkomentar sedikit lebih jauh tentang apa yang Anda lakukan dalam kode di atas?
Elvis
3
@ Elvis saya juga tidak mengerti. Interval kepercayaan berdasarkan kemungkinan profil harus dibangun dengan bantuan dari persentil, yang muncul di mana-mana. χ2
Stéphane Laurent
1
@ StéphaneLaurent Saya tidak yakin kode asli yang memberikan interval keyakinan. Interval kemungkinan maksimal 1/20. Saya percaya nama untuk interval kepercayaan dalam plot saya adalah interval kepercayaan "tipe-wald" dan garis merah pada plot adalah "kurva kepercayaan" yang dijelaskan pada halaman wikipedia ini
Flask

Jawaban:

10

Saya tidak akan memberikan jawaban yang lengkap (saya kesulitan memahami apa yang sebenarnya Anda lakukan), tetapi saya akan mencoba menjelaskan bagaimana kemungkinan profil dibuat. Saya dapat menyelesaikan jawaban saya nanti.

Kemungkinan penuh untuk sampel ukuran normal n adalah

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

Jika adalah parameter yang Anda minati, dan σ 2 adalah parameter gangguan, solusi untuk membuat kesimpulan hanya pada μ adalah untuk menentukan profil kemungkinan L P ( μ ) = L ( μ ,μσ2μ di mana ^ σ 2 (μ)adalah MLE untukμtetap: ^ σ 2 (μ)=argmaxσ2L(μσ

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")

kemungkinan profil

Tautkan dengan kemungkinan Saya akan mencoba menyorot tautan dengan kemungkinan dengan grafik berikut.

Pertama-tama tentukan kemungkinan:

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

Kemudian lakukan plot kontur:

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)

plot kontur L

Nilai-nilai kemungkinan profil adalah nilai yang diambil oleh kemungkinan sepanjang parabola merah.

μ^ adalah sama.

σ2^(μ) , tetapi selama Anda hanya berurusan dengan segmen pendek itu, itu hampir linier, dan perbedaannya akan sangat kecil .

Anda juga dapat menggunakan kemungkinan profil untuk membuat tes skor, misalnya.

Elvis
sumber
Mu dalam kode adalah urutan nilai dari rendah ke tinggi, kemungkinan pada masing-masing nilai ini dibagi dengan kemungkinan pada mean sampel (mn). Jadi itu adalah kemungkinan yang dinormalisasi.
Labu
Saya pikir ini adalah hal yang sama tetapi tidak dinormalisasi. Bisakah Anda memasukkannya ke dalam kode R atau memplot fungsi untuk beberapa data sehingga kami dapat membandingkan?
Labu
Ini dia. Awalnya saya pikir mnitu salah ketik, sekarang saya pikir kode R semuanya salah. Saya akan mengeceknya besok - sudah larut saya hidup.
Elvis
Kamu mungkin benar. Saya tidak mengerti bagaimana kode berhasil dinormalisasi. Oh, saya mengerti, "normalisasi" hanya membaginya dengan maksimum?
Elvis
1
Saya pikir ini adalah untuk membuatnya mudah untuk melihat ketika rasio kemungkinan kurang dari beberapa ambang batas (misalnya 1/20 maks) pada beberapa hipotesis nol (misalnya nol).
Labu
7

χk2

0.14795%

Ini adalah hasil klasik dan oleh karena itu saya hanya akan memberikan beberapa referensi tentang ini:

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

Kode R berikut menunjukkan bahwa, bahkan untuk sampel kecil, interval yang diperoleh dengan kedua pendekatan serupa (saya menggunakan kembali contoh Elvis):

Perhatikan bahwa Anda harus menggunakan kemungkinan profil yang dinormalisasi.

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

Jika kami menggunakan ukuran sampel yang lebih besar, interval kepercayaan lebih dekat:

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

TITIK PENTING:

Perhatikan bahwa untuk sampel tertentu, berbagai interval kepercayaan yang berbeda mungkin berbeda dalam hal panjang atau lokasi, yang sebenarnya penting adalah cakupannya. Dalam jangka panjang, mereka semua harus memberikan cakupan yang sama, secara independen pada seberapa besar perbedaan mereka untuk sampel tertentu.

Prokofiev
sumber
@Prokoflev jika ada beberapa hubungan sederhana antara interval kepercayaan yang dihitung dengan fungsi Rt.test () dan oleh yang dihitung dengan kode fungsi kemungkinan di atas, bisa Anda posting. Saya terutama tertarik pada kasus n = 3. Sayangnya saya memiliki sedikit latar belakang dalam matematika sehingga banyak makalah membawa saya ke lubang kelinci mencari nama untuk simbol dan apa yang mereka wakili dll ketika beberapa baris kode (paling mudah adalah R) dapat menjelaskannya kepada saya.
Labu
@Flask Apakah Anda tertarik untuk mendapatkan interval kepercayaan untuk parameter distribusi normal atau kerangka kerja yang lebih umum?
Prokofiev
@ Prokoflev khusus untuk rata-rata distribusi normal seperti yang ditunjukkan pada contoh saya dalam pertanyaan. Saya terutama bertanya-tanya mengapa interval kepercayaan lebih konservatif kecuali dalam kasus n = 3.
Labu
95%
1
Saya mulai percaya saya harus mengalikan interval kemungkinan dengan beberapa kuantil dari distribusi normal atau chisquare untuk mendapatkan interval kepercayaan yang sesuai ..
Flask
1

χ2nHairmSebuahlsayazed

  1. Log-kemungkinan profil adalah aproksimasi kuadratik
  2. Ada transformasi parameter yang membuat kemungkinan log profil lebih kurang kuadrat.

Kuadratik penting karena mendefinisikan distribusi normal dalam skala log. Semakin kuadratik, semakin baik perkiraan dan CI yang dihasilkan. Pilihan Anda untuk 1/20 cutoff untuk interval kemungkinan setara dengan lebih dari 95% CI dalam batas asimptotik, yang mengapa interval biru umumnya lebih lama daripada yang merah.

Sekarang, ada masalah lain dengan kemungkinan profil yang perlu mendapat perhatian. Jika Anda memiliki banyak variabel yang Anda profilingkan, maka jika jumlah titik data per dimensi rendah, kemungkinan profil bisa sangat bias dan optimis. Kemungkinan profil marjinal, bersyarat, dan dimodifikasi kemudian digunakan untuk mengurangi bias ini.

Jadi, jawaban untuk pertanyaan Anda adalah YA ... hubungannya adalah normalitas asimptotik dari kebanyakan penduga kemungkinan maksimum, seperti yang dimanifestasikan dalam distribusi chi-kuadrat dari rasio kemungkinan.


sumber
" Jika Anda memiliki banyak variabel yang Anda profilingkan, maka jika jumlah titik data per dimensi rendah, kemungkinan profil bisa sangat bias dan optimis " Optimis dibandingkan dengan apa?
Labu
@Flask Dengan optimis maksud saya akan terlalu sempit untuk memberikan probabilitas cakupan nominal saat memperlakukannya sebagai interval kepercayaan.
Saya mengerti, terima kasih, tetapi dalam kasus khusus saya itu sebenarnya pesimistis? Saya bingung tentang hal ini apakah kita berbicara tentang interval kemungkinan atau interval kepercayaan yang berasal dari kemungkinan.
Labu
@Flask Saya pikir Anda interval tampak pesimis karena Anda membandingkan interval kemungkinan 1/20 (kemungkinan relatif 5%) dengan 95% CI. Seperti yang dinyatakan oleh orang lain di sini, Anda benar-benar ingin membandingkannya dengan interval kemungkinan relatif 15% untuk memiliki apel ke apel ... setidaknya tanpa gejala. Interval kemungkinan Anda saat ini mempertimbangkan lebih banyak opsi sebagai masuk akal.
Saya telah merinci masalah aktual yang ingin saya terapkan pada apa yang saya pelajari di sini . Saya khawatir bahwa dalam kasus di mana distribusi pengambilan sampel tidak diketahui (tetapi tidak mungkin tidak normal) dan rumit bahwa dua persyaratan Anda mungkin tidak berlaku. Namun kemungkinan profil yang saya hitung tampaknya normal dan masuk akal. Apakah distribusi sampling rata - rata harus didistribusikan secara normal?
Labu