Какова связь между вероятностью профиля и доверительными интервалами?

18

Для построения этой диаграммы я сгенерировал случайные выборки разного размера из нормального распределения со средним значением = 0 и sd = 1. Затем были рассчитаны доверительные интервалы с использованием альфа-срезов в диапазоне от 0,001 до 0,999 (красная линия) с помощью функции t.test (), вероятность профиля была рассчитана с использованием кода, приведенного ниже, который я нашел в заметках к лекциям, помещенных в строку (я могу '). На данный момент ссылка не найдена. Редактировать: Найдено ), это показано синими линиями. Зеленые линии показывают нормализованную плотность с использованием функции R density (), а данные отображаются в виде прямоугольников в нижней части каждой диаграммы. Справа - гусеничный график с 95% доверительными интервалами (красный) и 1/20 максимальных правдоподобных интервалов (синий).

R код, используемый для вероятности профиля:

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

введите описание изображения здесь

Мой конкретный вопрос заключается в том, существует ли известная связь между этими двумя типами интервалов и почему доверительный интервал представляется более консервативным для всех случаев, кроме случаев, когда n = 3. Комментарии / ответы о том, верны ли мои расчеты (и лучший способ сделать это), и общие отношения между этими двумя типами интервалов также желательны.

Код 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)
колба
источник
В тебе лекция, mnэто опечатка mu, а не mean(dat). Как я уже говорил вам в комментариях к вашему другому вопросу , это должно быть ясно со страницы определений 23.
Элвис
@ Элвис Я так не думаю. Mn определено на странице 18 примечаний.
Настой
Я попытался уточнить понятие профиля вероятности. Можете ли вы прокомментировать немного больше того, что вы делаете в приведенном выше коде?
Элвис
3
@ Элвис Я тоже не понимаю. Доверительный интервал, основанный на вероятности профиля, должен быть построен с помощью процентилей , которые нигде не появляются. χ2
Стефан Лоран
1
@ StéphaneLaurent Я не уверен , что исходный код является обеспечение доверительных интервалов. Скорее всего 1/20 максимальных интервалов вероятности. Я считаю, что доверительные интервалы на моем графике называются доверительными интервалами типа «Вальда», а красные линии на графиках - это «доверительные кривые», описанные на этой странице википедии
Настой

Ответы:

10

Я не буду давать полный ответ (мне трудно понять, что именно вы делаете), но я постараюсь уточнить, как построена вероятность профиля. Я могу завершить свой ответ позже.

Полная вероятность для нормальной выборки размера : л ( μ , сг 2 ) = ( σ 2 ) - п / 2 ехр ( - Σ я ( х я - μ ) 2 / 2 σ 2 ) .n

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

Если - ваш интересующий параметр, а - неприятный параметр, решение сделать вывод только для - это определить вероятность профиля где - MLE для исправленного : μ ц л Р ( ц ) = L ( ц , ^ сг 2 ( ц ) ) ^ сг 2 ( ц ) ц ^ сг 2 ( ц ) = Argmax сг 2 л ( ц , сг 2 ) .σ2μ

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

Проверяется, что

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

Следовательно, вероятность профиля

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

Вот некоторый код R для вычисления и построения вероятности профиля (я удалил постоянный член ):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")

вероятность профиля

Связь с вероятностью Я постараюсь выделить ссылку с вероятностью на следующем графике.

Сначала определите вероятность:

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

Затем сделайте контурный сюжет:

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)

контурный участок L

Значения профиля вероятности являются значениями, взятыми по вероятности вдоль красной параболы.

Вы можете использовать профиль правдоподобия просто как одномерное классическое правдоподобие (cf @ Prokofiev's answer). Например, MLE - это то же самое.μ^

Для вашего доверительного интервала результаты будут немного отличаться из-за кривизны функции , но пока вы имеете дело только с ее коротким сегментом, она почти линейна, и разница будет очень мала.σ2^(μ)

Вы также можете использовать вероятность профиля, например, для построения тестов.

Элвис
источник
mu в коде - это последовательность значений от низкого до высокого, вероятность каждого из этих значений делится на вероятность в среднем по выборке (mn). Так что это нормализованная вероятность.
Настой
Я думаю, что это то же самое, но не нормализовано. Можете ли вы положить его в коде R или иным образом построить функцию для некоторых данных, чтобы мы могли сравнить?
Настой
Вот. Сначала я думал, что mnэто опечатка, но теперь я думаю, что код R все неправильно. Я дважды проверю это завтра - уже поздно, где я живу.
Элвис
Вы можете быть правы. Я не понимаю, как коду удается его нормализовать. О, я понял, "нормализация" просто делится на максимум?
Элвис
1
Я думаю, чтобы было легко увидеть, когда отношение правдоподобия меньше некоторого порога (например, 1/20 макс.) При некоторой нулевой гипотезе (например, ноль).
Настой
7

В общих рамках интервалы вероятности профиля являются приблизительными доверительными интервалами. Доказательство этого результата по существу аналогично доказательству того, что статистика отношения правдоподобия (асимптотически) приблизительно распределена как распределение . Идея состоит в том, чтобы инвертировать гипотезу, полученную из статистики отношения правдоподобия.χk2

Например, интервал вероятности профиля уровня имеет приблизительную достоверность .95 %0.14795%

Это классические результаты, и поэтому я просто приведу некоторые ссылки на это:

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

Следующий код R показывает, что даже для небольших выборок интервалы, полученные с помощью обоих подходов, аналогичны (я повторно использую пример Элвиса):

Обратите внимание, что вы должны использовать нормализованный профиль вероятности.

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

Если мы используем больший размер выборки, доверительные интервалы будут еще ближе:

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

ВАЖНАЯ ТОЧКА:

Обратите внимание, что для конкретных образцов различные виды доверительных интервалов могут различаться с точки зрения их длины или местоположения, что действительно важно, так это их охват. В конечном итоге все они должны обеспечивать одинаковое покрытие, независимо от того, насколько они различаются для конкретных образцов.

Прокофьев
источник
@Prokoflev, если есть какая-то простая связь между доверительными интервалами, вычисленными с помощью функции R t.test (), и теми, которые вычисляются с помощью приведенного выше кода функции правдоподобия, вы можете опубликовать его. Меня особенно интересует случай n = 3. К сожалению, я немного разбираюсь в математике, поэтому многие статьи ведут меня по кроличьей норе, просматривая названия символов, их символы и т. Д., Когда несколько строк кода (проще всего R) могут мне это объяснить.
Настой
@Flask Вы заинтересованы в получении доверительных интервалов для параметров нормального распределения или более общей структуры?
Прокофьев
@ Prokoflev специально для среднего нормального распределения, как показано в моем примере в вопросе. Мне особенно интересно, почему доверительные интервалы являются более консервативными, за исключением случая n = 3.
Настой
@Flask Какой уровень доверия вас интересует? ? 95%
Прокофьев
1
Я начинаю верить, что я должен умножать интервалы правдоподобия на некоторый квантиль нормального или квадратного распределения, чтобы получить соответствующий доверительный интервал.
Настой
1

Я не буду давать слишком математический ответ, но я хотел бы обратиться к вашему центральному вопросу о взаимосвязи между КИ и интервалами вероятности профиля. Как отмечали другие респонденты, КИ могут быть построены по профильной вероятности, используя приближение к отношению правдоподобия. Точность этого подхода зависит от того, является ли одна из двух вещей приблизительно верной: п о р м л я г е гχ2normalized

  1. Профиль журнала правдоподобия является приблизительно квадратичным
  2. Существует преобразование параметра, которое делает логарифмическую вероятность профиля приблизительно квадратичной.

Квадратик важен, потому что он определяет нормальное распределение в логарифмическом масштабе. Чем оно квадратичнее, тем лучше аппроксимация и результирующие КИ ». Ваш выбор 1/20-й отсечки для интервалов правдоподобия эквивалентен 95% -ному доверительному интервалу в асимптотическом пределе, поэтому синие интервалы обычно длиннее красных.

Теперь есть еще одна проблема с вероятностью профиля, которая требует некоторого внимания. Если у вас много переменных, для которых вы выполняете профилирование, то, если число точек данных на измерение мало, вероятность профиля может быть очень предвзятой и оптимистичной. Маргинальные, условные и модифицированные вероятности профиля затем используются для уменьшения этого смещения.

Итак, ответ на ваш вопрос - ДА ... связь является асимптотической нормальностью большинства оценок максимального правдоподобия, что проявляется в распределении хи-квадрат отношения правдоподобия.


источник
« Если у вас много переменных, для которых вы выполняете профилирование, то, если число точек данных в измерении мало, вероятность профиля может быть очень предвзятой и оптимистичной » Оптимистично по сравнению с чем?
Настой
@Flask Под оптимизмом я подразумеваю, что он будет слишком узким, чтобы обеспечить номинальную вероятность покрытия, рассматривая его как доверительный интервал.
Понятно, спасибо, но в моем конкретном случае это на самом деле пессимистично? Я смущен в этом вопросе относительно того, говорим ли мы об интервалах вероятности или доверительных интервалах, полученных из вероятностей.
Настой
@Flask Я думаю, что ваши интервалы кажутся пессимистичными, потому что вы сравниваете 1/20-й интервал правдоподобия (5% относительной вероятности) с 95% ДИ. Как утверждают другие здесь, вы действительно хотели бы сравнить его с 15% -ным относительным интервалом вероятности наличия яблок и яблок ... по крайней мере, асимптотически. Ваш вероятностный интервал в его нынешнем виде рассматривает больше вариантов как правдоподобных.
Я подробно изложил реальную проблему, к которой я хотел бы применить то, чему я учусь здесь . Я беспокоюсь о том, что в случае, когда распределение выборки неизвестно (но, вероятно, не нормально) и сложно, что ваши два требования могут не выполняться. Все же вероятности профиля, которые я вычислил, кажутся нормальными и разумными. Это то, что выборочное распределение среднего должно быть нормально распределено?
Настой