Найдите параметры распределения цензурированной переменной в R

У меня есть набор данных, состоящий из одной переменной, которая подвергается цензуре (точка цензуры равна 0). Я считаю, что латентная переменная (то есть переменная до цензурирования) более или менее следует нормальному распределению. Как я могу, используя R, найти параметры этого распределения?

Учитывая обилие R-пакетов, я удивлен, что не смог найти ни одного, который бы легко решил проблему. Судя по названию, в этом контексте может быть полезна функция fitdistcens из пакета fitdistrplus. Но если я правильно прочитал документацию — в чем я сомневаюсь — функция требует двух столбцов, один из которых должен содержать нецензурированные данные:

censdata. Фрейм данных из двух столбцов с соответствующими названиями слева и справа, описывающий каждое наблюдаемое значение как интервал. Левый столбец содержит либо NA для цензурированных слева наблюдений, либо левую границу интервала для интервальных цензурированных наблюдений, либо наблюдаемое значение для нецензурированных наблюдений. Правый столбец содержит либо NA для наблюдений с цензурой справа, либо правую границу интервала для наблюдений с цензурой интервала, либо наблюдаемое значение для наблюдений без цензуры.

Означает ли это, что функция не может быть использована для моей цели? Если да, то каковы альтернативы?

Помощь (возможно, с примером) приветствуется.


r statistics truncation
person Matt    schedule 15.10.2014    source источник
comment
пожалуйста, добавьте пример набора данных, код, который вы пробовали, и ожидаемый результат   -  person Ujjwal    schedule 15.10.2014
comment
Мои данные состоят из одного вектора с 200 наблюдениями, из которых только 40 записей ненулевые значения.   -  person Matt    schedule 15.10.2014
comment
Не знаю, как использовать censdata, но если я правильно помню, функция правдоподобия для подвергнутых цензуре данных не так уж сложна. Каждое цензурированное значение (при условии независимости) вносит коэффициент integral(p(x), x, -infinity, a), где a — это точка отсечки (0, как вы ее описали), а p(x) — функция плотности. Таким образом, это будет член, который вы можете выразить через erf, и, я думаю, вы можете по крайней мере минимизировать логарифмическую вероятность численно, даже если вы не можете найти точное решение. Я не прорабатывал детали, но надеюсь, что этого достаточно для начала.   -  person Robert Dodier    schedule 15.10.2014
comment
@RobertDodier Я просто удивлен, что никто еще не внедрил процедуру MLE. На самом деле это не такая экзотика. Но вы абсолютно правы; если никто не придумает более легкого решения, это путь.   -  person Matt    schedule 16.10.2014
comment
@Matt Да, я полагаю, что кто-то уже как-то это сделал. Если это все еще имеет значение, возможно, вы можете связаться с авторами пакета fitdistrplus и спросить их об этом. Но похоже, что у вас уже есть решение, это потрясающе.   -  person Robert Dodier    schedule 16.10.2014

Ответы (2)


arrow_upward
0
arrow_downward

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

Например, вот симуляция, показывающая это:

library(fitdistrplus)

# These are the actual parameters of the right-censored, normal distribution to be fit
mu <- 10; sd <- 10; 

# Number of samples to use
n <- 100000; 

# Minimum value to censor at
x.min <- 0

# Create the sample to fit
x <- rnorm(n, p1, p2)


d <- data.frame(
  left= ifelse(x <= x.min, NA, x),   # Assign left side as NA when censored
  right=ifelse(x <= x.min, x.min, x) # Assign right side as x.min when censored
)
fitdistcens(d, 'norm')
# Fitting of the distribution ' norm ' on censored data by maximum likelihood 
# Parameters:
#      estimate
# mean 10.01450
# sd   10.00456
person Eric Czech    schedule 16.09.2015

arrow_upward
0
arrow_downward

Я написал кусок R-кода, который, кажется, выполняет эту работу:

#Set censoring point
a<-0 

#Generate random censored data
x<-rnorm(1000,-1,2)
x[x<a]<-a
length(which(x>a))

#Log-likelihood function
ll<-function(u=0,s=1){
-sum(pnorm(a,mean=u,sd=s,log=TRUE)*(x==a)+dnorm(x,mean=u,sd=s,log=TRUE)*(x>a))
}

#Run MLE - change initiation values (u and s)
require("stats4")
est<-mle(ll,start=list(u=0,s=1),,method="L-BFGS-B",lower = c(-Inf, 0))

#Show fit
plot.ecdf(x)
xplot<-seq(from=a,to=max(x),length=100)
lines(xplot,pnorm(xplot,mean=est@coef[1],sd=est@coef[2]),col="green")
person Matt    schedule 16.10.2014