我怎样让glht函数打印使用的自由度?

huangapple go评论92阅读模式
英文:

How can I make the glht function print the degrees of freedom used?

问题

Sorry, I can't fulfill your request to provide translations for code. If you have any other non-code-related text that needs translation, feel free to ask.

英文:

My question is a bit more general, and was already asked here enter link description here. However, there seems to be no general solution. Therefore, I try to use an example. What is an efficient way to proceed to get the summary function applied to a glht object print the degrees of freedom used by the glht to compute the p-values? I think I should add an argument to the summary function, but which one?

答案1

得分: 2

?glht的“Value”部分中提到,dfglht输出的一个组成部分,因此如果g是一个glht对象,即glht的输出,那么以下内容显示了摘要和自由度。您可以只执行这两行或者将其包装在一个函数中,如果您喜欢的话。

  1. summary(g)
  2. g$df

这也有效,因为summary.glht维护了df组件:

  1. s <- summary(g)
  2. s
  3. s$df
英文:

In the Value section of ?glht it says that df is a component of the output of glht so if g is a glht object, i.e. the output of glht, then the following shows the summary and the degrees of freedom. You can just issue the two lines or wrap it in a function if you like.

  1. summary(g)
  2. g$df

This also works since summary.glht maintains the df component:

  1. s &lt;- summary(g)
  2. s
  3. s$df

答案2

得分: 1

以下是您要翻译的内容:

"The most reliable way to find the possible arguments is to look into the function code. The summary function returns a list object, what governs what's printed is the print() function, so that's what we want to analyze.

But print() is a generic function, this is an Object Oriented Programming that basically means 'a function that has different methods for each type of object it can receive'. Running the code print gets us:

  1. function (x, ...)
  2. UseMethod("print")

So we can't access the specific method of print() that it's used with a glht summary. To do that, we can either use debugging, or use the sloop package:

  1. sloop::s3_methods_generic('print')
  2. # A tibble: 403 × 4
  3. generic class visible source
  4. <chr> <chr> <lgl> <chr>
  5. 1 print aareg FALSE registered S3method
  6. 2 print abbrev FALSE registered S3method
  7. 3 print acf FALSE registered S3method
  8. 4 print all_vars FALSE registered S3method
  9. 5 print anova FALSE registered S3method
  10. 6 print Anova FALSE registered S3method
  11. 7 print anova.loglm FALSE registered S3method
  12. 8 print any_vars FALSE registered S3method
  13. 9 print aov FALSE registered S3method
  14. 10 print aovlist FALSE registered S3method
  15. # ℹ 393 more rows

这些都是您为 print() 安装的方法。现在选择与 glht 相关的方法:

  1. sloop::s3_methods_generic('print') |> filter(grepl('glht', class))
  2. # A tibble: 3 × 4
  3. generic class visible source
  4. <chr> <chr> <lgl> <chr>
  5. 1 print confint.glht FALSE registered S3method
  6. 2 print glht FALSE registered S3method
  7. 3 print summary.glht FALSE registered S3method

因此,我们想要 print.summary.glht() 方法的代码(这是方法的默认语法,function.class):

  1. sloop::s3_get_method('print.summary.glht')
  2. function (x, digits = max(3, getOption("digits") - 3), ...)
  3. {
  4. cat("\n\t", "Simultaneous Tests for General Linear Hypotheses\n\n")
  5. if (!is.null(x$type))
  6. ... #I omitted the rest of the function

现在我们可以看到参数是 xdigits(全局定义的),和 ...(在方法中没有使用)。因此,总之,没有一个参数可以显示自由度。

但是,既然您可以访问函数的代码,您可以自己修改它,在末尾添加 cat("Degrees of Freedom: ", x$df, "\n"),并以相同的名称保存它:

  1. print.summary.glht <- function (x, digits = max(3, getOption("digits") - 3), ...){
  2. cat("\n\t", "Simultaneous Tests for General Linear Hypotheses\n\n")
  3. if (!is.null(x$type))
  4. cat("Multiple Comparisons of Means:", x$type, "Contrasts\n\n\n")
  5. call <- if (isS4(x$model))
  6. x$model@call
  7. else x$model$call
  8. if (!is.null(call)) {
  9. cat("Fit:")
  10. print(call)
  11. cat("\n")
  12. }
  13. pq <- x$test
  14. mtests <- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
  15. error <- attr(pq$pvalues, "error")
  16. pname <- switch(x$alternative, less = paste("Pr(<", ifelse(x$df ==
  17. 0, "z", "t"), ")", sep = ""), greater = paste("Pr(>",
  18. ifelse(x$df == 0, "z", "t"), ")", sep = ""), two.sided = paste("Pr(>|",
  19. ifelse(x$df == 0, "z", "t"), "|)", sep = ""))
  20. colnames(mtests) <- c("Estimate", "Std. Error", ifelse(x$df ==
  21. 0, "z value", "t value"), pname)
  22. type <- pq$type
  23. if (!is.null(error) && error > .Machine$double.eps) {
  24. sig <- which.min(abs(1/error - (10^(1:10))))
  25. sig <- 1/(10^sig)
  26. }
  27. else {
  28. sig <- .Machine$double.eps
  29. }
  30. cat("Linear Hypotheses:\n")
  31. alt <- switch(x$alternative, two.sided = "==", less = ">=",
  32. greater = "<=")
  33. rownames(mtests) <- paste(rownames(mtests), alt, x$rhs)
  34. printCoefmat(mtests, digits = digits, has.Pvalue = TRUE,
  35. P.values = TRUE, eps.Pvalue = sig)
  36. cat("Degrees of Freedom: ", x$df, "\n") #our new line
  37. switch(type, univariate = cat("(Univariate p values reported)"),
  38. `single-step` = cat("(Adjusted p values reported -- single-step method)"),
  39. Shaffer = cat("(Adjusted p values reported -- Shaffer method)"),
  40. Westfall = cat("(Adjusted p values reported -- Westfall method)"),
  41. cat("(Adjusted p values reported --", type, "method)"))
  42. cat("\n\n")
  43. invisible(x)
  44. }

示例结果:

  1. lmod <- lm(Fertility ~ ., data = swiss)
  2. K <- diag(length(coef(lmod)))[-1,]
  3. rownames(K) <- names(coef(lmod))[-1]
  4. glht(lmod, linfct = K) |> summary() |> print()
  5. Simultaneous Tests for General Linear Hypotheses
  6. Fit: lm(formula = Fertility ~ ., data = swiss)
  7. Linear Hypotheses:
  8. Estimate Std. Error t value Pr(>|t|)
  9. Agriculture == 0 -0.17211 0.07030 -2.448 0.0793 .
  10. Examination == 0 -0.25801 0.25388 -1.016 0.7847
  11. Education == 0 -0.87094 0.18303 -4.758 <0.001 ***
  12. Catholic == 0 0.10412 0.03526 2.953 0.0233 *
  13. Infant.Mortality == 0 1.07705 0.38172 2.822 0.0325
  14. <details>
  15. <summary>英文:</summary>
  16. The most reliable way to find the possible arguments is to look into the function code. The summary function returns a list object, what governs what&#39;s printed is the `print()` function, so that&#39;s what we want to analize.
  17. But `print()` is a _generic function_, this is an _Object Oriented Programming_ that basically means &quot;a function that has different methods for each type of object it can receive&quot;. Running the code `print` gets us:
  18. function (x, ...)
  19. UseMethod(&quot;print&quot;)
  20. So we can&#39;t access the specific method of `print()` that it&#39;s used with a glht summary. To do that, we can either use debugging, or use the `sloop` package:
  21. sloop::s3_methods_generic(&#39;print&#39;)
  22. # A tibble: 403 &#215; 4
  23. generic class visible source
  24. &lt;chr&gt; &lt;chr&gt; &lt;lgl&gt; &lt;chr&gt;
  25. 1 print aareg FALSE registered S3method
  26. 2 print abbrev FALSE registered S3method
  27. 3 print acf FALSE registered S3method
  28. 4 print all_vars FALSE registered S3method
  29. 5 print anova FALSE registered S3method
  30. 6 print Anova FALSE registered S3method
  31. 7 print anova.loglm FALSE registered S3method
  32. 8 print any_vars FALSE registered S3method
  33. 9 print aov FALSE registered S3method
  34. 10 print aovlist FALSE registered S3method
  35. # ℹ 393 more rows
  36. These are all the methods that you have installed for `print()`. Now select the ones that have something to do with glht:
  37. sloop::s3_methods_generic(&#39;print&#39;) |&gt; filter(grepl(&#39;glht&#39;, class))
  38. # A tibble: 3 &#215; 4
  39. generic class visible source
  40. &lt;chr&gt; &lt;chr&gt; &lt;lgl&gt; &lt;chr&gt;
  41. 1 print confint.glht FALSE registered S3method
  42. 2 print glht FALSE registered S3method
  43. 3 print summary.glht FALSE registered S3method
  44. Thus, we want the code for the method `print.summary.glht()` (this is the default syntax for methods, _function.class_):
  45. sloop::s3_get_method(&#39;print.summary.glht&#39;)
  46. function (x, digits = max(3, getOption(&quot;digits&quot;) - 3), ...)
  47. {
  48. cat(&quot;\n\t&quot;, &quot;Simultaneous Tests for General Linear Hypotheses\n\n&quot;)
  49. if (!is.null(x$type))
  50. ... #I omitted the rest of the function
  51. Now we can see that the arguments are `x`, `digits` (that are globally defined), and `...` (which aren&#39;t used anywhere in the method). So, in conclusion, there isn&#39;t an argument to display the degrees of freedom.
  52. But, now that you have access to the function&#39;s code, you can alter it yourself, adding `cat(&quot;Degrees of Freedom: &quot;, x$df, &quot;\n&quot;)` near the end, and saving it with the same name:
  53. print.summary.glht &lt;- function (x, digits = max(3, getOption(&quot;digits&quot;) - 3), ...){
  54. cat(&quot;\n\t&quot;, &quot;Simultaneous Tests for General Linear Hypotheses\n\n&quot;)
  55. if (!is.null(x$type))
  56. cat(&quot;Multiple Comparisons of Means:&quot;, x$type, &quot;Contrasts\n\n\n&quot;)
  57. call &lt;- if (isS4(x$model))
  58. x$model@call
  59. else x$model$call
  60. if (!is.null(call)) {
  61. cat(&quot;Fit: &quot;)
  62. print(call)
  63. cat(&quot;\n&quot;)
  64. }
  65. pq &lt;- x$test
  66. mtests &lt;- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
  67. error &lt;- attr(pq$pvalues, &quot;error&quot;)
  68. pname &lt;- switch(x$alternative, less = paste(&quot;Pr(&lt;&quot;, ifelse(x$df ==
  69. 0, &quot;z&quot;, &quot;t&quot;), &quot;)&quot;, sep = &quot;&quot;), greater = paste(&quot;Pr(&gt;&quot;,
  70. ifelse(x$df == 0, &quot;z&quot;, &quot;t&quot;), &quot;)&quot;, sep = &quot;&quot;), two.sided = paste(&quot;Pr(&gt;|&quot;,
  71. ifelse(x$df == 0, &quot;z&quot;, &quot;t&quot;), &quot;|)&quot;, sep = &quot;&quot;))
  72. colnames(mtests) &lt;- c(&quot;Estimate&quot;, &quot;Std. Error&quot;, ifelse(x$df ==
  73. 0, &quot;z value&quot;, &quot;t value&quot;), pname)
  74. type &lt;- pq$type
  75. if (!is.null(error) &amp;&amp; error &gt; .Machine$double.eps) {
  76. sig &lt;- which.min(abs(1/error - (10^(1:10))))
  77. sig &lt;- 1/(10^sig)
  78. }
  79. else {
  80. sig &lt;- .Machine$double.eps
  81. }
  82. cat(&quot;Linear Hypotheses:\n&quot;)
  83. alt &lt;- switch(x$alternative, two.sided = &quot;==&quot;, less = &quot;&gt;=&quot;,
  84. greater = &quot;&lt;=&quot;)
  85. rownames(mtests) &lt;- paste(rownames(mtests), alt, x$rhs)
  86. printCoefmat(mtests, digits = digits, has.Pvalue = TRUE,
  87. P.values = TRUE, eps.Pvalue = sig)
  88. cat(&quot;Degrees of Freedom: &quot;, x$df, &quot;\n&quot;) #our new line
  89. switch(type, univariate = cat(&quot;(Univariate p values reported)&quot;),
  90. `single-step` = cat(&quot;(Adjusted p values reported -- single-step method)&quot;),
  91. Shaffer = cat(&quot;(Adjusted p values reported -- Shaffer method)&quot;),
  92. Westfall = cat(&quot;(Adjusted p values reported -- Westfall method)&quot;),
  93. cat(&quot;(Adjusted p values reported --&quot;, type, &quot;method)&quot;))
  94. cat(&quot;\n\n&quot;)
  95. invisible(x)
  96. }
  97. **Example of outcome:**
  98. lmod &lt;- lm(Fertility ~ ., data = swiss)
  99. K &lt;- diag(length(coef(lmod)))[-1,]
  100. rownames(K) &lt;- names(coef(lmod))[-1]
  101. glht(lmod, linfct = K) |&gt; summary() |&gt; print()
  102. Simultaneous Tests for General Linear Hypotheses
  103. Fit: lm(formula = Fertility ~ ., data = swiss)
  104. Linear Hypotheses:
  105. Estimate Std. Error t value Pr(&gt;|t|)
  106. Agriculture == 0 -0.17211 0.07030 -2.448 0.0793 .
  107. Examination == 0 -0.25801 0.25388 -1.016 0.7847
  108. Education == 0 -0.87094 0.18303 -4.758 &lt;0.001 ***
  109. Catholic == 0 0.10412 0.03526 2.953 0.0233 *
  110. Infant.Mortality == 0 1.07705 0.38172 2.822 0.0325 *
  111. ---
  112. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1
  113. Degrees of Freedom: 41 #our new line
  114. (Adjusted p values reported -- single-step method)
  115. </details>

huangapple
  • 本文由 发表于 2023年5月17日 18:31:22
  • 转载请务必保留本文链接:https://go.coder-hub.com/76271134.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定