#.Set Working Directory----------------------------------------------------------------------------------------------------
setwd("~/Cow/S21-MAW636 Analysis/") #Production Data Analysis/
setwd("~/Cow/S21-MAW636 Analysis/Production Data Analysis/")
# Load packages----------------------------------------------------------------------------------------------------.
library(readxl)
library(tidyverse)
library(lubridate)
library(lmerTest)
library(lme4)
library(broom.mixed)
library(emmeans)
library(ggplot2)
library(gtsummary)
library(afex)
library(broom.mixed)
# Upload data, by cow*period, filter production square, n = 32--------------------------------------------------------------------------------
# treatment schedule
trtsched_upload <- read_excel("MAW636 Treatment Schedule.xlsx") %>%
pivot_longer(cols = -Cow, names_to = "period", values_to = "trt") %>%
rename(cow = "Cow") %>%
mutate(cow = as.character(cow)) %>%
na_if("NA") %>%
drop_na() %>%
separate(trt, c("osc_lvl", "prot_lvl"), sep = cumsum(c(1, 2))) %>%
mutate(trt = paste0(osc_lvl, prot_lvl)) # separate into factorial trt
milk_comp_wt_summary <- read.csv("milk_comp_wt_summary.csv") %>%
mutate(cow = as.character(cow)) %>%
filter(period != "P2") %>%
filter(cow == "164" | cow == "8359" | cow == "8156" | cow == "148" |
cow == "7961" | cow == "9140" | cow == "145" | cow == "8254") %>%
mutate(cow = as.character(cow)) %>%
select(-X) %>%
left_join(trtsched_upload) %>%
mutate(osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP"))
ghgdat <- read_excel("~/2021-Fall/F21 GGAA/MAW636 Methane and CO2 data.xlsx") %>%
mutate(
period = paste0("P", period),
cow = as.character(id)
) %>%
select(-id) %>%
filter(cow == "164" | cow == "8359" | cow == "8156" | cow == "148" |
cow == "7961" | cow == "9140" | cow == "145" | cow == "8254") %>%
mutate(cow = as.character(cow)) %>%
mutate(osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP"))
DMI_dat_period <- read.csv("DMI_dat_period.csv") %>%
filter(cow == "164" | cow == "8359" | cow == "8156" | cow == "148" |
cow == "7961" | cow == "9140" | cow == "145" | cow == "8254") %>%
filter(period != "P2") %>%
select(-X) %>%
mutate(cow = as.character(cow))
# DMI - we will want treatment means instead.
DMI_dat_periodGG <- DMI_dat_period %>%
left_join(trtsched_upload) %>%
mutate(osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP"))
cp_avg <- read_excel("~/Cow/S21-MAW636 Analysis/Results from USDA/MAW636 CP Avg from Mary.xlsx")
# Join datasets -----------------------------------------------------------------------------
ghgdat_complete <- ghgdat %>%
left_join(milk_comp_wt_summary) %>%
left_join(DMI_dat_period, by = c("period", "cow"))
##%######################################################%##
# #
#### Calculate Emissions ####
# #
##%######################################################%##
# Production = g/d
# Intensity = g/kg FPCM milk
# Yield = g/kg DMI
ghgdat_calcs <- ghgdat_complete %>%
mutate(
ch_prod_g = ch4gd,
ch_intensity_g_kg = ch4gd / fpcm,
ch_yield_g_kg = ch4gd / DMI_kg,
co_prod_g = co2gd,
co_intensity_g_kg = co2gd / fpcm,
co_yield_g_kg = co2gd / DMI_kg,
o_prod_g = o2gd,
ch_co_ratio = (ch_prod_g/1000/0.657)/(co_prod_g/1000/1.977), # convert to liter
rq = (co2gd/1.977)/(o2gd/1.429) # divide by density (kg/m3) of each molecule at STP
) %>%
mutate(period = factor(period, levels = c("P1", "P2", "P3", "P4", "P5"))) %>%
mutate(osc_lvl = factor(osc_lvl, levels = c("S", "O"))) %>%
mutate(period_num = as.numeric(period)) %>% # create numeric period to use with CORCAR1
select(cow, period, period_num, osc_lvl, prot_lvl, starts_with("ch_"), starts_with("co_"), o_prod_g, rq)
write.csv(ghgdat_calcs, "MAW636_ghgdat_calcs.csv")Where \(w\) is the mass of milk, \(c\) is the proportion of component, \(i\) is cow, \(p\) is period, \(t\) is timepoint, and the \(\cdot\) is used to indicate indeces averaged over. This is attempting to express the procedure used to calculate the mean component production for each cow in each period. Missing observations were replaced with the mean of existing observations for that period, cow, and timepoint. A similar approach was used to calculate dry matter intake.
\[ x_{i p \space \cdot}\left(kg\right)= \frac{\sum_{t=1}^{n_{t} = 8} w_{i p t}\left(kg\right) c_{i p t}(\%)} {n_{t} = 8} \]
\[ x_{i p \space \cdot} (\%)=\sum_{t=1}^{n_{t}=8} \left( \frac{w_{i p t}\left(kg\right)}{w_{ i p \space \cdot}\left(kg\right)} \times\left(\operatorname{c_{ipt}}\left({\%} \right)\right) \right) \]
Using values from Mary Becker corrected for 105 DM. The mean and sd are of N = 80 observations representing N = 40 samples, each with 2 Leco replicates.
cp_avg %>%
group_by(Diet) %>%
summarise(mean = mean(Cpcorr),
sd = sd(Cpcorr))## # A tibble: 6 x 3
## Diet mean sd
## <chr> <dbl> <dbl>
## 1 A 12.2 0.332
## 2 B 13.9 0.607
## 3 C 15.5 0.758
## 4 Corn Silage 6.31 0.557
## 5 D 16.9 1.10
## 6 Haylage 20.5 1.03Set orthogonal (sum) contrasts and define convenience functions to make output of mixed() more readable/workable. The mixed() function is a wrapper for lmer(). Whereas lmer() by default outputs dummy (treatment) contrasts, mixed() outputs orthogonal (sum) contrasts with F tests based on Type III sums of squares.
# Set contrast options ----------------------------------------------------
afex::set_sum_contrasts()
# Define convenience functions to make output more readable/workable --------------------------------------------
a_mix <- function(mod) {
data.frame(mod$anova_table) %>%
rownames_to_column() %>%
select(rowname, Pr..F.) %>%
mutate(P = round(Pr..F., 3)) %>%
select(-Pr..F.)
}
e_mix <- function(mod) {
data.frame(emmeans(mod, ~ prot_lvl * osc_lvl)) %>%
mutate(across(where(is.numeric), round, 3))
} # Not used in 2022.01.17 due to computing overall mu when F test not significantIn the ANOVA-style notation, a general model can be expressed for the values of each dependent variable: \[ value_{ijkl} = \mu + period_{j} + protlevel_k + feedpattern_l + (protlevel_k \times feedpattern_l) + cow_i + \epsilon_{ijkl}\] where
\(\mu\) is the grand mean of the observations
\(period\) is the fixed effect of period (\(j\) = P1, P3, P4, P5)
\(protlevel\) is the fixed main effect of protein level (\(k\) = low, high)
\(feedpattern\) is the fixed main effect of feeding pattern (\(l\) = O, S)
\(protlevel_k \times feedpattern_l\) is the interaction of feeding pattern and protein level.
\(cow\) is the random effect of cow (\(i\) = 1,…,8)
For production models, data frames contain 32 observations, representing 8 cows in 4 periods. When no significant treatment effects, get overall mean of treatment means. This is different from the grand mean.
dmi_mod <- mixed(DMI_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = DMI_dat_periodGG)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]dmi_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: DMI_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: DMI_dat_periodGG
## Effect df F p.value
## 1 period 3, 18 13.73 *** <.001
## 2 prot_lvl 1, 18 0.74 .401
## 3 osc_lvl 1, 18 0.01 .935
## 4 prot_lvl:osc_lvl 1, 18 0.62 .442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(dmi_mod, "1") # emmean overall, on intercept. ## 1 emmean SE df lower.CL upper.CL
## overall 26.3 0.819 7 24.4 28.3
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95########################################################
dmi_mod$anova_table## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: DMI_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: DMI_dat_periodGG
## num Df den Df F Pr(>F)
## period 3 18 13.7270 6.688e-05 ***
## prot_lvl 1 18 0.7395 0.4011
## osc_lvl 1 18 0.0069 0.9345
## prot_lvl:osc_lvl 1 18 0.6187 0.4418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1data.frame(dmi_mod$anova_table)## num.Df den.Df F Pr..F.
## period 3 18 13.727041849 6.688112e-05
## prot_lvl 1 18 0.739523635 4.011225e-01
## osc_lvl 1 18 0.006937554 9.345387e-01
## prot_lvl:osc_lvl 1 18 0.618704426 4.417626e-01anmod <- tidy(dmi_mod$anova_table)
anmod## # A tibble: 4 x 5
## term num.Df den.Df statistic p.value
## <chr> <int> <dbl> <dbl> <dbl>
## 1 period 3 18.0 13.7 0.0000669
## 2 prot_lvl 1 18.0 0.740 0.401
## 3 osc_lvl 1 18.0 0.00694 0.935
## 4 prot_lvl:osc_lvl 1 18.0 0.619 0.442anmod %>% bind_rows## # A tibble: 4 x 5
## term num.Df den.Df statistic p.value
## <chr> <int> <dbl> <dbl> <dbl>
## 1 period 3 18.0 13.7 0.0000669
## 2 prot_lvl 1 18.0 0.740 0.401
## 3 osc_lvl 1 18.0 0.00694 0.935
## 4 prot_lvl:osc_lvl 1 18.0 0.619 0.442milk_mod <- mixed(milk_kgperd ~ period + prot_lvl * osc_lvl + (1 | cow), data = milk_comp_wt_summary)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]milk_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: milk_kgperd ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: milk_comp_wt_summary
## Effect df F p.value
## 1 period 3, 18.00 22.62 *** <.001
## 2 prot_lvl 1, 18.00 0.11 .745
## 3 osc_lvl 1, 18.00 0.98 .335
## 4 prot_lvl:osc_lvl 1, 18.00 0.22 .648
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1a_mix(milk_mod)## rowname P
## 1 period 0.000
## 2 prot_lvl 0.745
## 3 osc_lvl 0.335
## 4 prot_lvl:osc_lvl 0.648emmeans(milk_mod, "1")## 1 emmean SE df lower.CL upper.CL
## overall 36.7 1.28 7 33.6 39.7
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95fpcm_mod <- mixed(fpcm ~ period + prot_lvl * osc_lvl + (1 | cow), data = milk_comp_wt_summary)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]fpcm_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: fpcm ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: milk_comp_wt_summary
## Effect df F p.value
## 1 period 3, 18 26.94 *** <.001
## 2 prot_lvl 1, 18 0.01 .907
## 3 osc_lvl 1, 18 0.91 .353
## 4 prot_lvl:osc_lvl 1, 18 0.01 .910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1a_mix(fpcm_mod)## rowname P
## 1 period 0.000
## 2 prot_lvl 0.907
## 3 osc_lvl 0.353
## 4 prot_lvl:osc_lvl 0.910emmeans(fpcm_mod, "1")## 1 emmean SE df lower.CL upper.CL
## overall 36.6 1.33 7 33.5 39.7
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95For GHG models, data frames contain 32 observations, representing 8 cows in 4 periods.
# Methane Production g/d -----------------------------------------------------------------------------
ch_prod_mod01 <- mixed(ch_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]ch_prod_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: ch_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 9.75 *** <.001
## 2 prot_lvl 1, 18.00 0.02 .892
## 3 osc_lvl 1, 18.00 0.53 .477
## 4 prot_lvl:osc_lvl 1, 18.00 0.57 .461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(ch_prod_mod01, "1") ## 1 emmean SE df lower.CL upper.CL
## overall 474 23.9 7 417 530
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Methane Intensity g/kg FPCM -----------------------------------------------------------------------------
ch_intensity_mod01 <- mixed(ch_intensity_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]ch_intensity_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: ch_intensity_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 2.47 + .095
## 2 prot_lvl 1, 18.00 0.27 .612
## 3 osc_lvl 1, 18.00 0.01 .943
## 4 prot_lvl:osc_lvl 1, 18.00 0.47 .500
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(ch_intensity_mod01, "1")## 1 emmean SE df lower.CL upper.CL
## overall 13.1 0.557 7 11.7 14.4
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Methane Yield g/kg DMI -----------------------------------------------------------------------------
ch_yield_mod01 <- mixed(ch_yield_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]ch_yield_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: ch_yield_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 2.48 + .094
## 2 prot_lvl 1, 18.00 1.24 .280
## 3 osc_lvl 1, 18.00 0.49 .491
## 4 prot_lvl:osc_lvl 1, 18.00 0.01 .917
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(ch_yield_mod01, "1")## 1 emmean SE df lower.CL upper.CL
## overall 18 0.65 7 16.5 19.6
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Carbon dioxide Production g/d -----------------------------------------------------------------------------
co_prod_mod01 <- mixed(co_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]co_prod_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: co_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 20.29 *** <.001
## 2 prot_lvl 1, 18.00 3.70 + .070
## 3 osc_lvl 1, 18.00 2.62 .123
## 4 prot_lvl:osc_lvl 1, 18.00 3.34 + .084
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(co_prod_mod01, ~prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP S 14890 410 10.4 13980 15799
## LP S 14870 410 10.4 13961 15779
## HP O 15597 410 10.4 14688 16506
## LP O 14827 410 10.4 13918 15736
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Carbon dioxide Intensity g/kg FPCM -----------------------------------------------------------------------------
co_intensity_mod01 <- mixed(co_intensity_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]co_intensity_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: co_intensity_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 11.89 *** <.001
## 2 prot_lvl 1, 18.00 1.57 .226
## 3 osc_lvl 1, 18.00 0.03 .870
## 4 prot_lvl:osc_lvl 1, 18.00 1.34 .262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(co_intensity_mod01,"1")## 1 emmean SE df lower.CL upper.CL
## overall 418 15.6 7 381 455
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# Carbon dioxide Yield g/kg DMI -----------------------------------------------------------------------------
co_yield_mod01 <- mixed(co_yield_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]co_yield_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: co_yield_g_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18 11.85 *** <.001
## 2 prot_lvl 1, 18 6.34 * .022
## 3 osc_lvl 1, 18 0.54 .474
## 4 prot_lvl:osc_lvl 1, 18 0.01 .927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(co_yield_mod01, ~prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP S 586 16.3 16.3 551 620
## LP S 557 16.3 16.3 522 592
## HP O 595 16.3 16.3 561 630
## LP O 564 16.3 16.3 530 599
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# O2 Consumption g/d -----------------------------------------------------------------------------
o2_prod_mod01 <- mixed(o_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]o2_prod_mod01## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: o_prod_g ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18 16.34 *** <.001
## 2 prot_lvl 1, 18 2.98 .102
## 3 osc_lvl 1, 18 0.22 .646
## 4 prot_lvl:osc_lvl 1, 18 0.34 .565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(o2_prod_mod01,"1")## 1 emmean SE df lower.CL upper.CL
## overall 10362 296 7 9663 11062
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(o2_prod_mod01, ~prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP S 10421 332 10.8 9689 11152
## LP S 10223 332 10.8 9491 10955
## HP O 10603 332 10.8 9872 11335
## LP O 10202 332 10.8 9471 10934
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95\[ \mathrm{RQ}=\frac{\text { Volume of } \mathrm{CO}_{2} \text { emitted }}{\text { Volume of } \mathrm{O}_{2} \text { consumed }} \\ \\ \text{where the density of }\mathrm{CO}_{2} \text{ is 1.977 g/L} \text{ and} \text{ the density of }\mathrm{O}_{2} \text{ is 1.429 g/L} \]
# RQ -----------------------------------------------------------------------------
rq_mod <- mixed(rq ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]rq_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: rq ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18 3.43 * .039
## 2 prot_lvl 1, 18 0.03 .865
## 3 osc_lvl 1, 18 1.31 .268
## 4 prot_lvl:osc_lvl 1, 18 1.33 .263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(rq_mod, "1")## 1 emmean SE df lower.CL upper.CL
## overall 1.05 0.0121 7 1.02 1.08
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95# RQ -----------------------------------------------------------------------------
ch_co_ratio_mod <- mixed(ch_co_ratio ~ period + prot_lvl * osc_lvl + (1 | cow), data = ghgdat_calcs)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]ch_co_ratio_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: ch_co_ratio ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: ghgdat_calcs
## Effect df F p.value
## 1 period 3, 18.00 3.25 * .046
## 2 prot_lvl 1, 18.00 1.10 .309
## 3 osc_lvl 1, 18.00 0.00 .950
## 4 prot_lvl:osc_lvl 1, 18.00 0.48 .497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(ch_co_ratio_mod, "1")## 1 emmean SE df lower.CL upper.CL
## overall 0.0947 0.00398 7 0.0853 0.104
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95anova_proc <- function(model) {
result = tidy(model$anova_table)
depvar = as.character(terms(model$full_model)[[2]])
depvarc = rep(depvar, each = 4)
cbind(depvarc, result)
}
em_proc <- function(model) {
result <- data.frame(emmeans(model, ~ prot_lvl * osc_lvl)) %>%
mutate(across(where(is.numeric), round, 3))
depvar = as.character(terms(model$full_model)[[2]])
depvarc = rep(depvar, each = 4)
cbind(depvarc, result)
}
depvar = terms(dmi_mod$full_model)[[2]]
deparse(dmi_mod$full_model@call[[2]])
terms(dmi_mod$full_model)[[2]]
dmi_mod$anova_table
########################################################