Load libraries and download data
library(readxl)
library(tidyverse)
library(dplyr)
library(ordinal)
library(MASS)
library(brant)
library(bannerCommenter)
data<- read_excel('Emil.xlsx')
Descriptive statistics
library(tidyverse)
pre_columns <- c("Q1a pre", "Q1b Pre", "Q2 a&b pre", "Q2 C pre", "Q2 D pre", "Q3 Pre", "Q4 A&B pre")
post_columns <- c("Q1a Post", "Q1b Post", "Q2 a&b Post", "Q2 C Post", "Level Q2 D Post", "Level Q3 Post", "Level Q4 A&B post")
teacher_column <- "Teachers"
# Reshape data into long format
long_data <- data %>%
dplyr::select(all_of(c(pre_columns, post_columns, teacher_column))) %>%
pivot_longer(cols = c(pre_columns, post_columns),
names_to = "Test",
values_to = "Score") %>%
mutate(Test = dplyr::recode(Test,
`Q1a pre` = "Q1a Pre", `Q1a Post` = "Q1a Post",
`Q1b Pre` = "Q1b Pre", `Q1b Post` = "Q1b Post",
`Q2 a&b pre` = "Q2 a&b Pre", `Q2 a&b Post` = "Q2 a&b Post",
`Q2 C pre` = "Q2 C Pre", `Q2 C Post` = "Q2 C Post",
`Q2 D pre` = "Q2 D Pre", `Level Q2 D Post` = "Q2 D Post",
`Q3 Pre` = "Q3 Pre", `Level Q3 Post` = "Q3 Post",
`Q4 A&B pre` = "Q4 a&b Pre", `Level Q4 A&B post` = "Q4 a&b Post"))
# Create the summary statistics table
summary_table <- long_data %>%
group_by(Teachers, Test) %>%
summarise(
Mean = mean(Score, na.rm = TRUE),
SD = sd(Score, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_wider(names_from = Test, values_from = c(Mean, SD)) %>%
mutate(
`Q1a Change` = `Mean_Q1a Post` - `Mean_Q1a Pre`,
`Q1b Change` = `Mean_Q1b Post` - `Mean_Q1b Pre`,
`Q2ab Change` = `Mean_Q2 a&b Post` - `Mean_Q2 a&b Pre`,
`Q2C Change` = `Mean_Q2 C Post` - `Mean_Q2 C Pre`,
`Q2D Change` = `Mean_Q2 D Post` - `Mean_Q2 D Pre`,
`Q3 Change` = `Mean_Q3 Post` - `Mean_Q3 Pre`,
`Q4 Change` = `Mean_Q4 a&b Post` - `Mean_Q4 a&b Pre`
) %>%
dplyr::select(
Teachers,
`Mean_Q1a Pre`, `Mean_Q1a Post`, `Q1a Change`,
`Mean_Q1b Pre`, `Mean_Q1b Post`, `Q1b Change`,
`Mean_Q2 a&b Pre`, `Mean_Q2 a&b Post`, `Q2ab Change`,
`Mean_Q2 C Pre`, `Mean_Q2 C Post`, `Q2C Change`,
`Mean_Q2 D Pre`, `Mean_Q2 D Post`, `Q2D Change`,
`Mean_Q3 Pre`, `Mean_Q3 Post`, `Q3 Change`,
`Mean_Q4 a&b Pre`, `Mean_Q4 a&b Post`, `Q4 Change`
)
# Create the frequency table
frequencies <- long_data %>%
dplyr::group_by(Teachers, Test, Score) %>%
dplyr::count(name = "Count") %>%
dplyr::group_by(Teachers, Test) %>%
dplyr::mutate(Percent = Count / sum(Count)) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = Test, values_from = c(Count, Percent),
names_glue = "{Test}_{.value}") %>%
dplyr::arrange(Teachers, Score)
# Export to CSV
write.csv(summary_table, "descriptive_statistics_by_teacher.csv", row.names = FALSE)
write.csv(frequencies, "frequencies_by_teacher_test_etc.csv", row.names = FALSE)
Question 1a
long_dataq1a <- data%>%
dplyr::select("Q1a pre", "Q1a Post", "Teachers")%>%
pivot_longer(cols = c("Q1a pre", "Q1a Post"),
names_to = "Time",
values_to = "Q1aScore")
long_dataq1a$Q1aScore<- as.factor(long_dataq1a$Q1aScore)
Q1amodel <- polr(Q1aScore ~ Time + Teachers, data = long_dataq1a, method = "logistic")
model_summary <- summary(Q1amodel)
z_values <- model_summary$coefficients[, "Value"] / model_summary$coefficients[, "Std. Error"]
p_values <- 2 * (1 - pnorm(abs(z_values)))
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
model_summary
## Call:
## polr(formula = Q1aScore ~ Time + Teachers, data = long_dataq1a,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ1a pre -0.9065 0.7165 -1.265
## TeachersM -0.9506 0.6659 -1.428
##
## Intercepts:
## Value Std. Error t value
## 0|1 -4.4590 0.7928 -5.6247
## 1|2 -4.2251 0.7591 -5.5661
## 2|3 -3.4756 0.6815 -5.1000
##
## Residual Deviance: 87.24779
## AIC: 97.24779
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
p_values
## TimeQ1a pre TeachersM 0|1 1|2 2|3
## 2.058129e-01 1.534235e-01 1.858308e-08 2.604626e-08 3.396253e-07
banner("Proportional Odds")
##
## #################################################################
## ## Proportional Odds ##
## #################################################################
brant(Q1amodel)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 3.05 4 0.55
## TimeQ1a pre 1.68 2 0.43
## TeachersM 1.3 2 0.52
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
#Check balance between scores
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_dataq1a$Q1aScore)
##
## 0 1 2 3
## 4 1 5 130
#Check residuals
banner("predicted probability for residuals")
##
## #################################################################
## ## predicted probability for residuals ##
## #################################################################
predicted_probs <- predict(Q1amodel, type = "probs")
head(predicted_probs)
## 0 1 2 3
## 1 0.06901345 0.016633802 0.07975920 0.8345935
## 2 0.02907241 0.007384079 0.03766359 0.9258799
## 3 0.06901345 0.016633802 0.07975920 0.8345935
## 4 0.02907241 0.007384079 0.03766359 0.9258799
## 5 0.06901345 0.016633802 0.07975920 0.8345935
## 6 0.02907241 0.007384079 0.03766359 0.9258799
Question 1b
long_dataq1b <- data%>%
dplyr::select("Q1b Pre", "Q1b Post", "Teachers")%>%
pivot_longer(cols = c("Q1b Pre", "Q1b Post"),
names_to = "Time",
values_to = "Q1bScore")
long_dataq1b$Q1bScore<- as.factor(long_dataq1b$Q1bScore)
Q1bmodel <- polr(Q1bScore ~ Time + Teachers, data = long_dataq1b, method = "logistic")
model_summary <- summary(Q1bmodel)
z_values <- model_summary$coefficients[, "Value"] / model_summary$coefficients[, "Std. Error"]
p_values <- 2 * (1 - pnorm(abs(z_values)))
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
model_summary
## Call:
## polr(formula = Q1bScore ~ Time + Teachers, data = long_dataq1b,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ1b Pre -3.2298 0.4490 -7.194
## TeachersM -0.7333 0.3518 -2.085
##
## Intercepts:
## Value Std. Error t value
## 0|1 -5.9198 0.6060 -9.7679
## 1|2 -4.0797 0.4720 -8.6438
## 2|3 -1.2380 0.3026 -4.0913
## 3|4 -0.1546 0.2612 -0.5921
##
## Residual Deviance: 320.4833
## AIC: 332.4833
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
p_values
## TimeQ1b Pre TeachersM 0|1 1|2 2|3 3|4
## 6.292744e-13 3.710584e-02 0.000000e+00 0.000000e+00 4.290177e-05 5.538130e-01
banner("Proportional Odds")
##
## #################################################################
## ## Proportional Odds ##
## #################################################################
brant(Q1bmodel)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 18.8 6 0
## TimeQ1b Pre 18.01 3 0
## TeachersM 0.26 3 0.97
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
#Check balance between scores
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_dataq1b$Q1bScore)
##
## 0 1 2 3 4
## 6 22 57 19 36
#Check residuals
banner("predicted probability for residuals")
##
## #################################################################
## ## predicted probability for residuals ##
## #################################################################
predicted_probs <- predict(Q1bmodel, type = "probs")
head(predicted_probs)
## 0 1 2 3 4
## 1 0.123831231 0.34706156 0.4676028 0.03980437 0.02170002
## 2 0.005560597 0.02845281 0.3424316 0.26431596 0.35923902
## 3 0.123831231 0.34706156 0.4676028 0.03980437 0.02170002
## 4 0.005560597 0.02845281 0.3424316 0.26431596 0.35923902
## 5 0.123831231 0.34706156 0.4676028 0.03980437 0.02170002
## 6 0.005560597 0.02845281 0.3424316 0.26431596 0.35923902
Question 2 ab
banner("Running Q2a&b Model: Pre vs Post")
##
## ##################################################################
## ## Running Q2a&b Model: Pre vs Post ##
## ##################################################################
long_data_q2ab <- data %>%
dplyr::select("Q2 a&b pre", "Q2 a&b Post", Teachers) %>%
pivot_longer(cols = c("Q2 a&b pre", "Q2 a&b Post"), names_to = "Time", values_to = "Score") %>%
mutate(Score = as.factor(Score))
Q2ab_model <- polr(Score ~ Time + Teachers, data = long_data_q2ab, method = "logistic")
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
summary(Q2ab_model)
##
## Re-fitting to get Hessian
## Call:
## polr(formula = Score ~ Time + Teachers, data = long_data_q2ab,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ2 a&b pre -1.186 0.3733 -3.176
## TeachersM -1.218 0.3963 -3.074
##
## Intercepts:
## Value Std. Error t value
## 0|1 -4.0332 0.4952 -8.1448
## 1|2 -3.7800 0.4662 -8.1086
## 2|3 -2.6652 0.3794 -7.0254
## 3|4 0.7719 0.2741 2.8163
##
## Residual Deviance: 275.7718
## AIC: 287.7718
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
2 * (1 - pnorm(abs(coef(summary(Q2ab_model))[, "t value"])))
##
## Re-fitting to get Hessian
## TimeQ2 a&b pre TeachersM 0|1 1|2 2|3
## 1.491714e-03 2.110622e-03 4.440892e-16 4.440892e-16 2.134959e-12
## 3|4
## 4.858074e-03
banner("Proportional Odds Assumption")
##
## ##################################################################
## ## Proportional Odds Assumption ##
## ##################################################################
brant(Q2ab_model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 4.06 6 0.67
## TimeQ2 a&b pre 1.49 3 0.68
## TeachersM 2.3 3 0.51
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
## Warning in brant(Q2ab_model): 2 combinations in table(dv,ivs) do not occur.
## Because of that, the test results might be invalid.
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_data_q2ab$Score)
##
## 0 1 2 3 4
## 8 2 15 90 25
banner("Predicted Probabilities")
##
## #################################################################
## ## Predicted Probabilities ##
## #################################################################
head(predict(Q2ab_model, type = "probs"))
## 0 1 2 3 4
## 1 0.16390936 0.03771160 0.2334063 0.5248839 0.0400888
## 2 0.05651572 0.01511993 0.1188271 0.6892917 0.1202456
## 3 0.16390936 0.03771160 0.2334063 0.5248839 0.0400888
## 4 0.05651572 0.01511993 0.1188271 0.6892917 0.1202456
## 5 0.16390936 0.03771160 0.2334063 0.5248839 0.0400888
## 6 0.05651572 0.01511993 0.1188271 0.6892917 0.1202456
Question 2c
banner("Running Q2C Model: Pre vs Post")
##
## ##################################################################
## ## Running Q2C Model: Pre vs Post ##
## ##################################################################
long_data_q2c <- data %>%
dplyr::select("Q2 C pre", "Q2 C Post", Teachers) %>%
pivot_longer(cols = c("Q2 C pre", "Q2 C Post"), names_to = "Time", values_to = "Score") %>%
mutate(Score = as.factor(Score))
Q2c_model <- polr(Score ~ Time + Teachers, data = long_data_q2c, method = "logistic")
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
summary(Q2c_model)
## Call:
## polr(formula = Score ~ Time + Teachers, data = long_data_q2c,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ2 C pre -0.5688 0.3461 -1.6435
## TeachersM -0.1770 0.3758 -0.4709
##
## Intercepts:
## Value Std. Error t value
## 0|1 -3.0480 0.4196 -7.2646
## 1|2 -2.1574 0.3370 -6.4026
## 2|3 -0.8773 0.2852 -3.0761
##
## Residual Deviance: 278.6206
## AIC: 288.6206
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
2 * (1 - pnorm(abs(coef(summary(Q2c_model))[, "t value"])))
## TimeQ2 C pre TeachersM 0|1 1|2 2|3
## 1.002765e-01 6.377241e-01 3.741452e-13 1.527534e-10 2.097310e-03
banner("Proportional Odds Assumption")
##
## ##################################################################
## ## Proportional Odds Assumption ##
## ##################################################################
brant(Q2c_model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 8.42 4 0.08
## TimeQ2 C pre 1.62 2 0.44
## TeachersM 6.75 2 0.03
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_data_q2c$Score)
##
## 0 1 2 3
## 9 11 32 88
banner("Predicted Probabilities")
##
## #################################################################
## ## Predicted Probabilities ##
## #################################################################
head(predict(Q2c_model, type = "probs"))
## 0 1 2 3
## 1 0.09093684 0.10503873 0.2711709 0.5328536
## 2 0.05360445 0.06766975 0.2104531 0.6682727
## 3 0.09093684 0.10503873 0.2711709 0.5328536
## 4 0.05360445 0.06766975 0.2104531 0.6682727
## 5 0.09093684 0.10503873 0.2711709 0.5328536
## 6 0.05360445 0.06766975 0.2104531 0.6682727
Question 2d
banner("Running Q2D Model: Pre vs Post")
##
## ##################################################################
## ## Running Q2D Model: Pre vs Post ##
## ##################################################################
long_data_q2d <- data %>%
dplyr::select("Q2 D pre", "Level Q2 D Post", Teachers) %>%
pivot_longer(cols = c("Q2 D pre", "Level Q2 D Post"), names_to = "Time", values_to = "Score") %>%
mutate(Score = as.factor(Score))
Q2d_model <- polr(Score ~ Time + Teachers, data = long_data_q2d, method = "logistic")
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
summary(Q2d_model)
## Call:
## polr(formula = Score ~ Time + Teachers, data = long_data_q2d,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ2 D pre -1.3039 0.3387 -3.8493
## TeachersM -0.3559 0.3757 -0.9475
##
## Intercepts:
## Value Std. Error t value
## 0|1 -3.8931 0.4762 -8.1750
## 1|2 -0.6560 0.2687 -2.4414
## 2|3 1.2768 0.2968 4.3017
## 3|4 1.7606 0.3371 5.2227
##
## Residual Deviance: 321.5268
## AIC: 333.5268
## (1 observation deleted due to missingness)
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
2 * (1 - pnorm(abs(coef(summary(Q2d_model))[, "t value"])))
## TimeQ2 D pre TeachersM 0|1 1|2 2|3 3|4
## 1.184352e-04 3.433902e-01 2.220446e-16 1.463081e-02 1.694992e-05 1.763216e-07
banner("Proportional Odds Assumption")
##
## ##################################################################
## ## Proportional Odds Assumption ##
## ##################################################################
brant(Q2d_model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 7.58 6 0.27
## TimeQ2 D pre 0.01 3 1
## TeachersM 7.57 3 0.06
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_data_q2d$Score)
##
## 0 1 2 3 4
## 7 67 47 6 12
banner("Predicted Probabilities")
##
## #################################################################
## ## Predicted Probabilities ##
## #################################################################
head(predict(Q2d_model, type = "probs"))
## 0 1 2 3 4
## 1 0.09680244 0.6349993 0.2178257 0.01870785 0.03166469
## 2 0.02827410 0.3972651 0.4110089 0.05595019 0.10750170
## 3 0.09680244 0.6349993 0.2178257 0.01870785 0.03166469
## 4 0.02827410 0.3972651 0.4110089 0.05595019 0.10750170
## 5 0.09680244 0.6349993 0.2178257 0.01870785 0.03166469
## 6 0.02827410 0.3972651 0.4110089 0.05595019 0.10750170
Q3
banner("Running Q3 Model: Pre vs Post")
##
## #################################################################
## ## Running Q3 Model: Pre vs Post ##
## #################################################################
long_data_q3 <- data %>%
dplyr::select("Q3 Pre", "Level Q3 Post", Teachers) %>%
pivot_longer(cols = c("Q3 Pre", "Level Q3 Post"), names_to = "Time", values_to = "Score") %>%
mutate(Score = as.factor(Score))
Q3_model <- polr(Score ~ Time + Teachers, data = long_data_q3, method = "logistic")
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
summary(Q3_model)
## Call:
## polr(formula = Score ~ Time + Teachers, data = long_data_q3,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ3 Pre -1.6653 0.3753 -4.437
## TeachersM -0.6446 0.3866 -1.667
##
## Intercepts:
## Value Std. Error t value
## 0|1 -3.6284 0.4439 -8.1733
## 1|2 -0.1140 0.2640 -0.4317
## 2|3 0.8614 0.2858 3.0140
## 3|4 2.3769 0.4433 5.3620
##
## Residual Deviance: 310.2686
## AIC: 322.2686
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
2 * (1 - pnorm(abs(coef(summary(Q3_model))[, "t value"])))
## TimeQ3 Pre TeachersM 0|1 1|2 2|3 3|4
## 9.123358e-06 9.547156e-02 2.220446e-16 6.659281e-01 2.577894e-03 8.230123e-08
banner("Proportional Odds Assumption")
##
## ##################################################################
## ## Proportional Odds Assumption ##
## ##################################################################
brant(Q3_model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 16.89 6 0.01
## TimeQ3 Pre 6.33 3 0.1
## TeachersM 9.53 3 0.02
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_data_q3$Score)
##
## 0 1 2 3 4
## 13 84 21 16 6
banner("Predicted Probabilities")
##
## #################################################################
## ## Predicted Probabilities ##
## #################################################################
head(predict(Q3_model, type = "probs"))
## 0 1 2 3 4
## 1 0.2110650 0.6888209 0.05985454 0.03112747 0.009132025
## 2 0.0481603 0.5814660 0.18883730 0.13507015 0.046466193
## 3 0.2110650 0.6888209 0.05985454 0.03112747 0.009132025
## 4 0.0481603 0.5814660 0.18883730 0.13507015 0.046466193
## 5 0.2110650 0.6888209 0.05985454 0.03112747 0.009132025
## 6 0.0481603 0.5814660 0.18883730 0.13507015 0.046466193
Question 4a&b
# Prepare long-form data for Q4 A&B
long_dataq4 <- data %>%
dplyr::select("Q4 A&B pre", "Level Q4 A&B post", "Teachers") %>%
pivot_longer(cols = c("Q4 A&B pre", "Level Q4 A&B post"),
names_to = "Time",
values_to = "Q4Score")
# Convert scores to factor
long_dataq4$Q4Score <- as.factor(long_dataq4$Q4Score)
# Fit ordinal logistic regression model
Q4model <- polr(Q4Score ~ Time + Teachers, data = long_dataq4, method = "logistic")
# Model summary
model_summary <- summary(Q4model)
# Compute p-values
z_values <- model_summary$coefficients[, "Value"] / model_summary$coefficients[, "Std. Error"]
p_values <- 2 * (1 - pnorm(abs(z_values)))
banner("Model Summary")
##
## #################################################################
## ## Model Summary ##
## #################################################################
model_summary
## Call:
## polr(formula = Q4Score ~ Time + Teachers, data = long_dataq4,
## method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## TimeQ4 A&B pre -1.9997 0.3590 -5.571
## TeachersM -0.4402 0.3563 -1.235
##
## Intercepts:
## Value Std. Error t value
## 0|1 -3.9427 0.4415 -8.9298
## 1|2 -1.2081 0.2951 -4.0940
## 2|3 0.5574 0.2695 2.0684
## 3|4 2.3853 0.4394 5.4286
##
## Residual Deviance: 343.5313
## AIC: 355.5313
banner("P Values")
##
## ##################################################################
## ## P Values ##
## ##################################################################
p_values
## TimeQ4 A&B pre TeachersM 0|1 1|2 2|3
## 2.537040e-08 2.166947e-01 0.000000e+00 4.240488e-05 3.859844e-02
## 3|4
## 5.680497e-08
banner("Proportional Odds Test")
##
## ##################################################################
## ## Proportional Odds Test ##
## ##################################################################
brant(Q4model)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 19.04 6 0
## TimeQ4 A&B pre 14.46 3 0
## TeachersM 3.61 3 0.31
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
# Check score balance
banner("Balance Between Scores")
##
## ##################################################################
## ## Balance Between Scores ##
## ##################################################################
table(long_dataq4$Q4Score)
##
## 0 1 2 3 4
## 12 59 42 21 6
# Check predicted probabilities
banner("Predicted Probabilities for Residuals")
##
## #################################################################
## ## Predicted Probabilities for Residuals ##
## #################################################################
predicted_probs <- predict(Q4model, type = "probs")
head(predicted_probs)
## 0 1 2 3 4
## 1 0.1820009 0.5921220 0.1783277 0.03958788 0.007961545
## 2 0.0292406 0.2876862 0.4136592 0.21345065 0.055963357
## 3 0.1820009 0.5921220 0.1783277 0.03958788 0.007961545
## 4 0.0292406 0.2876862 0.4136592 0.21345065 0.055963357
## 5 0.1820009 0.5921220 0.1783277 0.03958788 0.007961545
## 6 0.0292406 0.2876862 0.4136592 0.21345065 0.055963357