在R中繪製3D對數似然圖

Question

我想嘗試使用R爲具有兩個參數的模型繪製3D對數似然圖。

我的功能如下：

negL=function(theta,data){ 
    mu=theta[1] 
    tau=theta[2] 

    n11=data$ev.trt 
    n12=data$n.trt-data$ev.trt 
    n21=data$ev.ctrl 
    n22=data$n.ctrl-data$ev.ctrl 
    odds=(n11*n22)/(n12*n21) 
    logodds=log(odds) 
    varodds=(1/n11)+(1/n12)+(1/n21)+(1/n22) 

    x=logodds 
    sigma=varodds 

    L=-(1/2)*sum(log((2*pi)*(sigma+tau)))-(1/2)*sum((((x-mu)^2))/(sigma+tau)) 
    return(-L) 
    } 

這裏L是我的似然函數。

我想嘗試和借鑑基於此數據對數似然情節：

  Author ev.ctrl n.ctrl ev.trt n.trt duration 
1   RAPID1  15 199 146 393 6.150 
2   RAPID2  4 127  80 246 5.850 
3  Kim2007  9  63  28 65 6.850 
4   DE019  19 200  81 207 10.950 
5   ARMADA  5  62  37 67 11.650 
6 Weinblatt1999  1  30  23 59 13.000 
7   START  33 363 110 360 8.100 
8   ATTEST  22 110  61 165 7.850 
9  Abe2006  1  47  15 49 8.300 
10 Strand2006  5  40  5 40 11.250 
11  CHARISMA  14  49  26 50 0.915 
12  OPTION  22 204  90 205 7.650 

我已經計算theta的最大似然估計是（1.6760569 0.2112315），使用R的Optim（）

基本上我想讓mu出現在x軸上，tau在y軸上，然後在z軸上出現對數似然-L。當我嘗試和運行它

Error in (function (theta) : 
    unused argument (c(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, 
-0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, -0.96, 

：

我試着看到它的另一網頁後添加該到我的功能

persp3d(negL, c(0,6), c(-1,3), col="red", xlab="mu", ylab="tau^2", zlab="Likelihood", 
     main="Perspective plot of Likelihood") 

從RGL包，但我得到一個錯誤說。

任何想法如何做到這一點和功能使用將不勝感激。

非常感謝

Answer 1

這是你想要實現的嗎？

negL<- function(mu,tau,data){ 
    n11=data$ev.trt 
    n12=data$n.trt-data$ev.trt 
    n21=data$ev.ctrl 
    n22=data$n.ctrl-data$ev.ctrl 
    odds=(n11*n22)/(n12*n21) 
    logodds=log(odds) 
    varodds=(1/n11)+(1/n12)+(1/n21)+(1/n22) 
    x=logodds 
    sigma=varodds 
    L=-(1/2)*sum(log((2*pi)*(sigma+tau)))-(1/2)*sum((((x-mu)^2))/(sigma+tau)) 
    return(-L) 
} 
f <- Vectorize(negL,vectorize.args=c("mu","tau")) 
library(rgl) 
mu <- seq(0,6,0.1) 
tau <- seq(0,3,0.1) 
z <- outer(mu,tau,f,data=data) 
zlim <- range(z[!is.na(z)]) 
palette <- rev(rainbow(20)) 
colors <- palette[19*(z-zlim[1])/diff(zlim) + 1] 

persp3d(mu,tau,z,col=colors) 

所以有幾個問題。

首先，persp3d(...)接受3個參數：x是定義x軸的向量，y是定義y軸的矢量，並且z是矩陣與尺寸長度（X）×長度（y），其中對於x和y的組合，每個元素的值爲z。您必須設置矢量x和y以及矩陣z（如上面的代碼中所示）。

設置z矩陣的一種方法是使用outer(x,y,FUN)函數（請閱讀文檔）。這需要兩個向量，x和y，並返回在每個組合x和y上評估函數FUN的結果。但它需要FUN爲矢量化，您的negL(...)不是。所以我們創建一個新功能，f(...)這是一個矢量化版本negL(...)。

三，功能negL(...)必須採取兩個參數，mu和tau，而不是與兩個元件（theta）的載體。所以我改變了函數的簽名來做到這一點。除此之外，該功能與您的問題完全相同。

最後，在你的問題中，你（似乎......）想要使用負值tau。但negL(...)未定義爲tau < 0。