當我們只有關於當前素食者的調查數據時，如何計算平均堅持素食主義的時間？
對隨機人口樣本進行了調查。他們被問到是否吃素食。如果他們的回答是肯定的，他們還被要求說明他們連續吃素食多久了。我想用這個數據來計算堅持素食的平均時間。換句話說，當有人成為素食主義者時，我想知道他們平均保持素食多久。讓我們假設：
- 所有受訪者都給出了正確和準確的回答
- 世界是穩定的：素食主義的流行度沒有改變，平均堅持時間也沒有改變。
到目前為止我的推理
我發現分析一個世界的玩具模型很有幫助，每年年初都有兩個人成為素食者。每次，他們中的一個人吃素1年，另一個人吃素3年。顯然，這個世界的平均堅持時間是（1 + 3）/ 2 = 2 年。這是一個說明該示例的圖表。每個矩形代表一個素食主義時期：
假設我們在第 4 年年中進行了一項調查（紅線）。我們得到以下數據：
如果我們在任何一年（從第 3 年開始）進行調查，我們都會得到相同的數據。如果我們只是平均我們得到的回答：
(2* 0.5 + 1.5 + 2.5)/4 = 1.25
我們低估了，因為我們假設每個人在調查後都不再吃素，這顯然是不正確的。為了獲得更接近這些參與者保持素食的真實平均時間的估計值，我們可以假設他們平均報告了大約一半的素食時間，並將報告的持續時間乘以 2。在一項隨機抽取的大型調查中從人口（就像我正在分析的人一樣）來看，我認為這是一個現實的假設。至少它會給出正確的期望值。但是，如果我們唯一要做的事情是加倍，我們得到的平均值是 2.5，這是一個高估。這是因為一個人保持素食的時間越長，他就越有可能成為當前素食者的樣本。
然後我認為某人在當前素食者樣本中的概率與他們的素食時間成正比。為了解釋這種偏差，我試圖將當前素食者的數量除以他們預計的堅持時間：
但是，這也給出了不正確的平均值：
(2*1 + ⅓ * 3 + ⅕ * 5)/(2 + ⅓ + ⅕) = 4 / 2.533333 = 1.579 年
如果素食者的數量除以他們正確的堅持時間，它將給出正確的估計：
(1 + ⅓ * (1 + 3 + 5))/(1 + ⅓ * 3) = 2 年
但是，如果我使用預測的依從性長度並且它們是我在現實中所擁有的一切，那麼它就不起作用了。我不知道還能嘗試什麼。我讀了一些關於生存分析的文章，但我不確定如何在這種情況下應用它。理想情況下，我還希望能夠計算 90% 的置信區間。任何提示將非常感謝。
編輯：上面的問題可能沒有答案。但也有另一項研究詢問隨機樣本的人是否/曾經是素食者，以及他們過去吃過多少次素食。我也知道每個人在學習和其他方面的年齡。也許這些信息可以與對當前素食者的調查結合使用，以某種方式獲得平均值。實際上，我談到的研究只是其中的一部分，但卻是一個非常重要的問題，我想從中得到更多。
讓 fX(x) 表示依從長度的pdf X 人口中的素食主義。我們的目標是估計 EX=∫∞0xfX(x);dx .
假設被納入調查的概率（事件 S ) 與 X , 粘附長度的 pdf X 參與調查的有 fX|S(x)=xfX(x)∫xfX(x)dx=xfX(x)EX.
在被納入調查時，只有一次 Z 已通過。有條件的 X 和 S ，據報導吃素的時間與pdf統一 fZ|X=x,S(z)=1x,for 0≤z≤x.其他地方的密度為零。因此，利用總概率定律，時間的整體分佈 Z 被納入調查的人成為素食者 fZ|S(z)=∫∞0fZ|X=x,S(z)fX|S(x)dx&=∫∞z1xxfX(x)EXdx&=1−FX(z)EX,在哪裡 FX(z) 是的 cdf X . 自從 X 是一個正變量 FX(0)=P(X≤0)=0 所以 fZ(0)=1/EX .這表明估計 EX 通過也許首先估計 fZ|S(z) 從觀察到的數據非參數 z1,z2,…,zn . 一種選擇是核密度估計，使用Silverman 的反射 方法 z=0 由於域 fZ(z) 有一個下界 z=0 . 該方法應用於模擬數據，如下圖紅色曲線所示。已獲得估價 ˆfZ(0) 的 fZ(z) 在 z=0 ，估計 EX 然後由 ^EX=1/ˆfZ(0) .
然而，這種非參數方法並不理想，因為它沒有利用以下事實： fZ|S(z) 是非增函數。另外，如果 fX(0)=F′X(0)>0 , fZ|S(0) 可能會被嚴重低估和 EX 高估了。找到一個估計 EX 在這種情況下，不做更多假設似乎很困難，主要是因為在這種情況下存在的短遵守時間幾乎不會出現在觀察到的數據中，這是有偏抽樣的結果。
或者，可以做出一些關於分佈的假設 fX(x) 並通過最大化似然性來擬合參數模型 L(θ)=n∏i=11−FX(zi;θ)EX(θ)
數字（上圖中的藍色曲線）。R代碼模擬數據並實現這兩種方法：
# Simulate lognormal duration length in population set.seed(1) n <- 1e+4 x <- rlnorm(n,mean=2,sd=.2) # Biased sampling y <- sample(x, size=n/10, prob=x, replace=TRUE) # Duration at time of sampling z <- runif(length(y),min=0, max=y) hist(z,prob=TRUE,main="") # Compute kernel density estimate with reflection aroudn z=0 to <- max(x) + 3 fhat <- density(z,from = -to, to=to) m <- length(fhat$y) fhat$y <- fhat$y[(m/2+1):m] + fhat$y[(m/2):1] fhat$x <- fhat$x[(m/2+1):m] lines(fhat,col="red") # Estimate of EX 1/fhat$y[1] #> [1] 7.032448 # True value (mean of above lognormal) exp(2+.2^2/2) #> [1] 7.538325 # Maximum likelihood nll <- function(theta, z) { - sum(plnorm(z, theta[1], theta[2], log.p=TRUE, lower.tail = FALSE)) + length(z)*(theta[1] + theta[2]^2/2) } fit <- optim(c(0,1),nll,z=z) fit$par #> [1] 1.9860464 0.2019859 EXhat <- exp(fit$par[1]+fit$par[2]^2/2) # MLE of EX EXhat #> [1] 7.436836 curve(plnorm(z, fit$par[1], fit$par[2], lower.tail=FALSE)/EXhat, xname="z", col="blue",add=TRUE)
