Improving_MultComp/Tutorial on Categories.Rmd at main · Leitemfa/Improving_MultComp · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
---
title: "Categories of Model Results"
output: html_notebook
editor_options:
  chunk_output_type: console
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook to explain how to use the code to generate interpretation about the significant and non-significant effects based on the paper **the paper: I know what I shouldn’t say, but what should I say? An approach to present results without statistical significance** (download [here](https://doi.org/10.3389/fams.2024.1487750)). I am not an author of the paper, but I got really interested in their view

When you execute code within the notebook, the results appear in your console.

I want to make it in a way that you can extract their output from the simple lm model and also for multiple comparisons. For that last one, I made a code from the `multcomp` package.

# Loading libraries
```{r}
library(devtools)
library(multcomp)
```

# Running the simple lm model

For this example I will use the classic data from the `cars` library

```{r}
#loading the source code
source("Improving_MultComp/Categorize_lm.R")

# Example usage
lm_model <- lm(mpg ~ wt+hp+qsec+vs+am, data = mtcars)
summary(lm_model)
View(categorize_association(lm_model))

lm_model <- lm(mpg ~ -1 +wt+hp+qsec+vs+am, data = mtcars)
summary(lm_model)
View(categorize_association(lm_model))

# Output: Detailed summary with classifications for each coefficient

```

As you can see the output now extracted the p-value and the confidence interval for the estimates of every regression coefficient. Then, we can look at its range and create a classification for each estimate from the linear model.

We can also make a plot out of this output

```{r}
source("plot_association.R")
plot_association(lm_model)

```


Next step: made the code to also include mixed models

# Multiple comparisons

As I said I also want to make it so we can extract the same type of classification from multiple comparisons. For that I have choosen to start with the `multcomp` package

```{r}
source("Categories_glht.R")

### multiple linear model, swiss data
lmod <- lm(Fertility ~ ., data = swiss)

### test of H_0: all regression coefficients are zero
### (ignore intercept)

### define coefficients of linear function directly
K <- diag(length(coef(lmod)))[-1,]
rownames(K) <- names(coef(lmod))[-1]
K

### set up general linear hypothesis
lmod.glht <- glht(lmod, linfct = K)
str(lmod.glht)
summary(lmod.glht)
View(extract_glht_ci(lmod.glht))
plot_glht_ci(lmod.glht)
### alternatively, use a symbolic description
### instead of a matrix
lmod.glht2 <- glht(lmod, linfct = c("Agriculture = 0",
                      "Examination = 0",
                      "Education = 0",
                      "Catholic = 0",
                      "Infant.Mortality = 0"))
extract_glht_ci(lmod.glht2)
plot_glht_ci(lmod.glht2)

### multiple comparison procedures
### set up a one-way ANOVA
amod <- aov(breaks ~ tension, data = warpbreaks)

### set up all-pair comparisons for factor `tension'
### using a symbolic description (`type' argument
### to `contrMat()')
lmod.glht3 <- glht(amod, linfct = mcp(tension = "Tukey"))
extract_glht_ci(lmod.glht3)
plot_glht_ci(lmod.glht3)

### alternatively, describe differences symbolically
lmod.glht4 <- glht(amod, linfct = mcp(tension = c("M - L = 0",
                                    "H - L = 0",
                                    "H - M = 0")))
extract_glht_ci(lmod.glht3)
plot_glht_ci(lmod.glht3)

### mix of one- and two-sided alternatives
warpbreaks.aov <- aov(breaks ~ wool + tension,
                      data = warpbreaks)

### contrasts for `tension'
K <- rbind("L - M" = c( 1, -1,  0),
           "M - L" = c(-1,  1,  0),
           "L - H" = c( 1,  0, -1),
           "M - H" = c( 0,  1, -1))

warpbreaks.mc <- glht(warpbreaks.aov,
                      linfct = mcp(tension = K),
                      alternative = "less")
confint(warpbreaks.mc, level = 0.95)$confint
extract_glht_ci(warpbreaks.mc)
plot_glht_ci(warpbreaks.mc)

```