2010-09-01 38 views
4

在下面的代碼中,我使用引導來計算C.I.以及在零假設下施用於番茄植物的兩種不同肥料對植物產量沒有影響(並且「改良」肥料更好的替代方案)下的p值。第一個隨機樣本(x)來自使用標準肥料的植物,而第二個樣本(y)來自的植物則使用「改良」樣本。引導以比較兩個組

x <- c(11.4,25.3,29.9,16.5,21.1) 
y <- c(23.7,26.6,28.5,14.2,17.9,24.3) 
total <- c(x,y) 
library(boot) 
diff <- function(x,i) mean(x[i[6:11]]) - mean(x[i[1:5]]) 
b <- boot(total, diff, R = 10000) 

ci <- boot.ci(b) 
p.value <- sum(b$t>=b$t0)/b$R 

我不喜歡上面的代碼是什麼好像有11個值(分離第一5爲屬於X離開其餘樣本y採樣)的僅一個樣品重採樣完成。 您可以告訴我如何修改此代碼,以便從第一個樣本中取出大小爲5的重新取樣,並從第二個樣本中取出大小爲6的重新取樣,以便引導重新取樣可以模擬生成的「單獨樣本」設計原始數據?

回答

5

EDIT2: 刪除了Hack,因爲這是一個錯誤的解決方案。相反,人們必須使用開機功能的參數階層:

total <- c(x,y) 
id <- as.factor(c(rep("x",length(x)),rep("y",length(y)))) 
b <- boot(total, diff, strata=id, R = 10000) 
... 

要知道你不會得到甚至接近你的p.value一個正確的估計:

x <- c(1.4,2.3,2.9,1.5,1.1) 
y <- c(23.7,26.6,28.5,14.2,17.9,24.3) 

total <- c(x,y) 

b <- boot(total, diff, strata=id, R = 10000) 
ci <- boot.ci(b) 
p.value <- sum(b$t>=b$t0)/b$R 
> p.value 
[1] 0.5162 

會如何你解釋了兩個樣本的p值爲0.51,其中第二個值的所有值都高於第一個值的最高值?

上述代碼很好地獲得置信區間的偏置估計,但是關於差異的重要性測試應該通過對整個數據集進行置換來完成。

+0

謝謝!關於你的問題「我爲什麼要這麼做」,請查看第18頁底部的標題爲「用於比較兩個種羣的引導程序」(如果您想對此進行評論),請參閱http://bcs.whfreeman.com /ips5e/content/cat_080/pdf/moore14.pdf – 2010-09-01 08:55:04

+0

我的主要問題是如何定義diff.calc。我很驚訝地看不到裏面的第二個參數! – 2010-09-01 09:00:44

+0

@ gd047:我猜你已經從statexchange上的問題中得到了類似的東西。請注意,他們只說一個置信區間,並沒有提到那裏的p值。我的例子顯示你爲什麼。 – 2010-09-01 09:10:49

1

雖然在某些情況下實際土壤牀可能被認爲是分層變量,但這不是其中之一。你只有一個操作,在植物組之間。因此,你的零假設是他們確實來自完全相同的人口。處理這些項目就好像它們來自一組11個樣本,這是在這種情況下引導的正確方法。

如果你有兩個地塊,並且在每個地塊以平衡的方式嘗試過不同季節的不同肥料,那麼這些地塊將是分層樣本,你會希望像這樣對待它們。但這不是這種情況。

1

繼約翰,我覺得合適的方式來使用引導程序來測試這兩個不同羣體的總和是顯著不同的是如下:

x <- c(1.4,2.3,2.9,1.5,1.1) 
y <- c(23.7,26.6,28.5,14.2,17.9,24.3) 


b_x <- boot(x, sum, R = 10000) 
b_y <- boot(y, sum, R = 10000) 

z<-(b_x$t0-b_y$t0)/sqrt(var(b_x$t[,1])+var(b_y$t[,1])) 
pnorm(z) 

因此,我們可以清楚地拒絕它們是空相同的人口。我可能錯過了一定程度的自由度調整,但我不確定bootstrapping在這方面的工作原理,但這種調整不會顯着改變您的結果。