私のカスタム関数で "subscript out of bounds"エラーが発生しましたか？

Question

こんにちは、すべて。そこで私は、一部のlm()で生成された回帰の推定値（たとえば、H_0：Beta_j = "some constant" vs H_1：Beta_j < "一定の定数"）でt検定を実行する簡単なカスタム関数を書くことにしました。 これは自分のコードを初めて作成するときですが、数ヶ月間Rを使って作業していましたが、私はそれをよく理解していると思いますので、どうして "subscript out of bounds"それを実行しています。

私のコード：。

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") 
    } 

} 

そして、私はhrsempはVARされ、標準lm()結果でこれを実行しようとすると、エラーが発生します。

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見つけることができるである理由の多くの良い説明の一つ：あなたは、おそらくより多くのこのような何かをしたいです。