「下標越界」錯誤對我的自定義函數？

Question

大家好。因此，我決定嘗試寫一個簡單的自定義函數，對某些產生的迴歸估計量進行t檢驗（例如H_0：Beta_j =「某些常量」與H_1：Beta_j <「某些常量」）。 這是我第一次創建自己的代碼，但我一直與R一起工作了幾個月，我想我對它有一個體面的理解，所以我不明白爲什麼我不斷收到「下標越界」運行它。

我的代碼：

custom_test<-function(data,coeff,alt,alternative=c("two.sided","greater","less"),clevel=.95){ 
    dof<-data$df.residual 
    top<-data$coefficients["coeff"]-alt 
    bottom=coef(summary(data))["coeff","Std. Error"] 
    stat<-abs(top/bottom) 
    if (alternative=="two.sided") { 
    tstat<-qt(clevel/2,dof) 
    pstat<-2*pt(tstat,dof) 
    return(pstat) 
    } else if (alternative=="greater") { 
     tstat<-qt(clevel/2,dof) 
     pstat<-pt(tstat,dof) 
     return(pstat) 
    } else if (alternative=="less") { 
     tstat<-qt(clevel/2,dof) 
     pstat<-pt(tstat,dof) 
     return(pstat) 
    } else { 
     return("Error") 
    } 

} 

我嘗試用標準lm()結果運行此，hrsemp是一個變種，並且得到錯誤：

custom_test(fit9,hrsemp,0,alternative="less") 
Error in coef(summary(data))["coeff", "Std. Error"] : 
    subscript out of bounds 

但每次我運行有問題的代碼手動自己，我確實得到一個答案：

> coef(fit9) 
(Intercept)  hrsemp log(sales) log(employ) 
12.45837237 -0.02926893 -0.96202698 0.76147045 
> coef(summary(fit9))["hrsemp", "Std. Error"] 
[1] 0.02280484 

其他關於這個錯誤的堆棧交換問題似乎有微妙的不同，而且迄今爲止我還沒有能夠將他們的經驗總結爲我的代碼。

有人請解釋我要去哪裏嗎？

Answer 1

Frank是對的;基於與其他所有人相同的原因，您會得到此錯誤：您嘗試訪問不存在的對象元素。更具體地說，在你的情況下，你試圖訪問"coeff"行和"Std. Error"列中的元素coef(summary(data))。這是一個問題，因爲可能不存在名爲"coeff"的係數。你要做到以下幾點：

custom_test<-function(data,coeff,alt,alternative=c("two.sided","greater","less"),clevel=.95){ 
    dof<-data$df.residual 
    top<-data$coefficients[coeff]-alt 
    bottom=coef(summary(data))[coeff,"Std. Error"] 
    stat<-abs(top/bottom) 
    if (alternative=="two.sided") { 
     tstat<-qt(clevel/2,dof) 
     pstat<-2*pt(tstat,dof) 
     return(pstat) 
    } else if (alternative=="greater") { 
     tstat<-qt(clevel/2,dof) 
     pstat<-pt(tstat,dof) 
     return(pstat) 
    } else if (alternative=="less") { 
     tstat<-qt(clevel/2,dof) 
     pstat<-pt(tstat,dof) 
     return(pstat) 
    } else { 
     return("Error") 
    } 

} 

和傳遞變量的字符串名字：

set.seed(42) 
hrsemp <- rnorm(10) 
Y <- 1 + 5 * hrsemp + rnorm(10) 
fit9 <- lm(Y ~ hrsemp) 
custom_test(fit9, 'hrsemp', 0, alternative="less") 
[1] 0.475 

（請注意，你可以交替喂函數的實際變量對象，並使用deparse(substitute(coeff)) - 例如，見this SO question）。

但是，您可能會注意到這會給您錯誤的答案。這是因爲你寫錯了你的功能。你可能想要更多的東西是這樣的：

custom_test <- function(data, coeff, alt, 
         alternative = c("two.sided", "greater", "less"), 
         clevel = .95){ 
    dof <- data$df.residual 
    top <- data$coefficients[coeff] - alt 
    bottom <- coef(summary(data))[coeff, "Std. Error"] 
    stat <- abs(top/bottom) 
    if (alternative == "two.sided") { 
     return(2 * (1 - pt(stat, dof))) 
    } else if (alternative == "greater") { 
     return(1 - pt(stat, dof)) 
    } else if (alternative == "less") { 
     return(1 - pt(stat, dof)) 
    } else { 
     stop("Provide a valid alternative hypothesis.", call.=FALSE) 
    } 
} 


custom_test(fit9, 'hrsemp', 0, alternative="less") 
     hrsemp 
7.858176e-05 
custom_test(fit9, 'hrsemp', 0, alternative="two.sided") 
     hrsemp 
0.0001571635 
coef(summary(fit9))['hrsemp', 'Pr(>|t|)'] 
[1] 0.0001571635 

一個爲什麼這是正確的計算可以發現here很多很好的解釋。