Skip to content

Commit

Permalink
emmean chapter changes
Browse files Browse the repository at this point in the history
  • Loading branch information
Sidhuharp97 committed Dec 3, 2024
1 parent 9f511d7 commit b7218cb
Show file tree
Hide file tree
Showing 2 changed files with 81 additions and 47 deletions.
98 changes: 56 additions & 42 deletions chapters/means-and-contrasts.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -61,12 +61,12 @@ Now that we have a good model, let’s use the emmeans() function to obtain EMMs
Now that we have fitted a linear mixed model (model1) and it meets the model assumption. Let's use the `emmeans()` function to obtain estimated marginal means for main (variety and nitrogen) and interaction (variety x nitrogen) effects.

#### Main effects
```{r}
```{r, warning=FALSE}
m1 <- emmeans(model1, ~V, level = 0.95)
m1
```

```{r}
```{r, warning=FALSE}
m2 <- emmeans(model1, ~N)
m2
```
Expand Down Expand Up @@ -97,36 +97,20 @@ pairs(m2)
Here if we look at the results from code chunk above, it's easy to interpret results from `pairs()` function in case of varietey comparison becuase there were only 3 groups. It's bit confusing in case of Nitrogen treatments where we had 4 groups. We can further simplify it by using custom contrasts.


::: {.callout-caution}
## `pairs()`
Remember!!
The `pairs()` function can be used to calculate pairwise comparison when treatment groups are less than equal to 3.
:::
### Custom contrasts

1. Treatment

Firstly, let's run emmean object 'm1' for variety compariosn.
Firstly, let's run emmean object 'm2' for nitrogen treatment comparison.

```{r}
m2
```
```{r, echo=FALSE}
#emmeans(model1, specs = trt.vs.ctrlk ~ N)
```



```{r}
coef(pairs(m1))
```

Now, lets a create a vector for each variety in the same order as presented in output from m1.

Now, lets a create a vector for each nitrogen treatment in the same order as presented in output from m2.
```{r}
A1 = c(1, 0, 0)
A2 = c(0, 1, 0)
A3 = c(0, 0, 1)
```


```{r, echo =FALSE}
A1 = c(1, 0, 0, 0)
A2 = c(0, 1, 0, 0)
A3 = c(0, 0, 1, 0)
Expand All @@ -138,7 +122,6 @@ Let's create custom contrasts for comparing '0.0cwt' (A1) treatment to '0.2cwt'
This can be evaluated as shown below:
Here the output shows the difference in mean yield between these two varieties
```{r}
pairs(m2, adjust = "tukey")
contrast(m2, method = list(A1 - A2) )
contrast(m2, method = list(A1 - A3) )
contrast(m2, method = list(A1 - A4) )
Expand All @@ -149,13 +132,6 @@ Using custom contrasts is strongly recommended instead of `pairs()` when you are
:::


::: {.callout-caution}
## `pairs()`
Remember!!
The `pairs()` function can be used to calculate pairwise comparison when treatment groups are less than equal to 3.
:::


## Compact letter displays
Compact letter displays (CLDs) are a popular way to display multiple comparisons when there are more than few group means to compare. However, they are problematic as they are more prone to misinterpretation.
The R package `multcompView` (Graves et al., 2019) provides an implementation of CLDs creating a display where any two means associated with same symbol are not statistically different.
Expand All @@ -179,9 +155,49 @@ cld(m3, alpha=0.05, Letters=letters)
```
Interpretation of these letters is: Here we have a significant difference in grain yield with varieties "victory", with N treatments of 0.0cwt, 0.2cwt, 0.4cwt, and 0.6wt. Grain yield for Golden.rain variety was significantly lower with 0.0cwt N treatment compared to the 0.2cwt, 0.4cwt, and 0.6wt treatments.


In the data set we used for demonstration here, we had equal number of observations in each group. However, this might not be a case every time as it is common to have missing values in the data set. In such cases, readers usually struggle to interpret significant differences among groups. For example, estimated means of two groups are substantially different but they are no statistically different. This normally happens when SE of one group is large due to its small sample size, so it's hard for it to be statistically different from other groups. In such cases, we can use alternatives to CLDs as shown below.

## Alternatives to CLD

1. Equivalence test

Let's assume based on subject matter considerations, if mean yield of two groups differ by less than 30 can be considered equivalent.
Let's try equivalence test on clds of nitrogen treatment emmeans (m2)
```{r}
cld(m2, delta = 30, adjust = "none")
```
Here, two treatment groups '0.0cwt' and '0.2cwt', '0.4cwt' and '0.6cwt' can be considered equivalent.

2. Significance Sets

Another alternative is to simply reverse all the boolean flags we used in constructing CLDs for m3 first time.


```{r}
cld(m2, signif = TRUE)
```
## Export emmeans to excel sheet

The emmean objects can be directly exported

### start from here
```{r}
result_n <- as.data.frame(summary(m1))
```


```{r}
library(writexl)
write_xlsx(result_n)
```


```{r}
```



```{r, eval=FALSE, echo = FALSE}
One-way estimated marginal means and plot
Expand All @@ -198,14 +214,6 @@ CLD = cld(marginal,
CLD
Location emmean SE df lower.CL upper.CL .group
Olympia 8.333333 0.6718548 16 6.449596 10.21707 a
Northampton 11.833333 0.6718548 16 9.949596 13.71707 b
Ventura 13.333333 0.6718548 16 11.449596 15.21707 b
Burlington 21.833333 0.6718548 16 19.949596 23.71707 c
### Order the levels for printing
CLD$Location = factor(CLD$Location,
Expand Down Expand Up @@ -258,6 +266,12 @@ ggplot(CLD,
```




## creating plots with emmeans



### Interactions using Emmeans

```{r, eval=FALSE, echo=FALSE}
Expand Down
Loading

0 comments on commit b7218cb

Please sign in to comment.