363_PSET_4
2025-06-16
In this analysis I will be focusing on players from top 5 leagues in the world according to UEFA Country Coefficients. I will be loading players only from the top 5 leagues and removing all players that are goalkeepers (GK), as the attributes that determine their ratings are singular, as they only affect the GK postion. I will be examining the relationship between attributes such as Acceleration, Sprint Speed, Positioning, Finishing, Shot Power, Long Shots, Volleys, Penalties, Vision, Crossing, Free Kick Accuracy, Short Passing, Long Passing, Curve, Dribbling, Agility, Balance, Reactions, Ball Control, Composure, Interceptions, Heading Accuracy, Defensive Awareness, Standing Tackle, Sliding Tackle, Jumping, Stamina, Strength and Aggression to determine the players overalls.
df <- read.csv("male_players.csv")
df <- df %>%
filter(Position != "GK")
df <- df %>%
filter(League %in% c("Premier League", "LALIGA EA SPORTS",
"Serie A Enilive", "Bundesliga", "Ligue 1 McDonald's"))ggplot(df, aes(x = League, y = OVR, fill = League)) + geom_violin(alpha = 0.7,
trim = FALSE) + geom_jitter(aes(color = League), width = 0.2,
size = 1.2, alpha = 0.6, show.legend = FALSE) + scale_fill_brewer(palette = "Paired") +
coord_flip() + labs(title = "Violin + Jitter: OVR Distribution by League",
x = "League", y = "Overall Rating (OVR)") + theme_minimal() +
theme(legend.position = "none", axis.text.y = element_text(size = 10),
axis.title = element_text(size = 12), plot.title = element_text(hjust = 0.5,
size = 14, face = "bold"), panel.grid.major.y = element_blank())
library(tidyverse)
## Assumes your raw data frame is called df
## ------------------------------- If not, change the
## object name in the next few lines.
# 1 – Create a broader Position grouping (helps keep cell
# sizes healthy)
df <- df |>
mutate(PositionGroup = case_when(Position %in% c("ST", "CF",
"LW", "RW") ~ "Attack", Position %in% c("CAM", "CM",
"CDM", "LM", "RM") ~ "Midfield", Position %in% c("CB",
"LB", "RB", "LWB", "RWB") ~ "Defense", TRUE ~ "Other") |>
factor(levels = c("Attack", "Midfield", "Defense")))
# 2 – Create Age groups for the age-interaction example
df <- df |>
mutate(AgeGroup = cut(Age, breaks = c(-Inf, 22, 29, Inf),
labels = c("Young", "Prime", "Veteran")))
# 3 – Choose the response variables you care about
resp <- c("PAC", "SHO", "PAS", "DRI", "DEF", "PHY", "OVR")# Plot mean + SE bar charts for each attribute
for (v in resp) {
ggplot(df, aes(x = PositionGroup, y = .data[[v]], fill = PositionGroup)) +
stat_summary(fun = mean, geom = "col", colour = "black",
width = 0.7) + stat_summary(fun.data = mean_se, geom = "errorbar",
width = 0.2) + labs(title = paste(v, "by Position"),
x = NULL, y = v) + theme_minimal(base_size = 13) + theme(legend.position = "none") ->
p
print(p)
}
# Optional – table of group means (helps your discussion
# section)
df |>
group_by(PositionGroup) |>
summarise(across(all_of(resp), ~round(mean(.x, na.rm = TRUE),
1)))## # A tibble: 3 × 8
## PositionGroup PAC SHO PAS DRI DEF PHY OVR
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Attack 75 72.7 63.6 73.4 35.5 68.2 73.9
## 2 Midfield 70.5 65.3 69.8 73.9 57.7 66.2 73.4
## 3 Defense 68.6 46.2 61.7 64.9 72 72.2 73.5# Classic interaction plots (base R style)
for (v in resp) {
interaction.plot(df$PositionGroup, # x-axis factor
df$League, # trace factor
df[[v]], # response
fun = mean,
type = "b", pch = 19,
col = rainbow(length(unique(df$League))),
xlab = "Position group",
ylab = paste("Mean", v),
main = paste("Position × League:", v))
}
# Tip: If you prefer ggplot aesthetics ----------------------------
# library(ggh4x) ## (for nicer facet strips; optional)
for (v in resp) {
ggplot(df,
aes(x = PositionGroup, y = .data[[v]],
colour = League, group = League)) +
stat_summary(fun = mean, geom = "line") +
stat_summary(fun = mean, geom = "point") +
labs(title = paste("Position × League:", v),
x = "Position group", y = v) +
theme_minimal(base_size = 13) -> p
print(p)
}
Per these visualisations I have decided to go with the first option of conducting One-Way MANOVA by Primary Position, as we can see in the interaction plots made for the other two options, that there are limited examples of significant interaction effects.
# Check multivariate normality via chi‐square Q–Q plots for
# each PositionGroup
# vector of your multivariate responses
resp <- c("PAC", "SHO", "PAS", "DRI", "DEF", "PHY", "OVR")
# Attackers
dat_attack <- df[df$PositionGroup == "Attack", resp, drop = FALSE]
cqplot(dat_attack, main = "Attack")
# Midfielders
dat_midfield <- df[df$PositionGroup == "Midfield", resp, drop = FALSE]
cqplot(dat_midfield, main = "Midfield")
# Defenders
dat_defense <- df[df$PositionGroup == "Defense", resp, drop = FALSE]
cqplot(dat_defense, main = "Defense")
We can see from the plots above that there is multivariate normality for the attack and defense groups, but the midfield group shows some deviations from normality, this is not suprising as the midfield group both includes midfielders that play on the wings, which prioritise pace,dribbliung and shooting, and midfielders that play centrally, which prioritise passing, defending and physicality.
# –– 0. Setup libraries and data prep ––#
library(car) # for Anova(), linearHypothesis(), qqPlot()
library(candisc) # for cqplot()## Warning: package 'candisc' was built under R version 4.3.3##
## Attaching package: 'candisc'## The following object is masked from 'package:vegan':
##
## scores## The following object is masked from 'package:stats':
##
## cancorlibrary(vegan) # for mrpp()
# assume your data is already in df and Position is a
# factor of interest and responses are the core FIFA
# attributes:
resp <- c("PAC", "SHO", "PAS", "DRI", "DEF", "PHY", "OVR")
mv_form <- as.formula(paste0("cbind(", paste(resp, collapse = ","),
") ~ PositionGroup"))
manova_mod <- lm(mv_form, data = df)
# Fit the MANOVA
manova_mod <- lm(mv_form, data = df)
# multivariate tests:
pillai <- anova(manova_mod)
wilks <- anova(manova_mod, test = "Wilks")
roy <- anova(manova_mod, test = "Roy")
print(pillai)## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99416 55983 7 2304 < 2.2e-16 ***
## PositionGroup 2 1.09654 400 14 4610 < 2.2e-16 ***
## Residuals 2310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(wilks)## Analysis of Variance Table
##
## Df Wilks approx F num Df den Df Pr(>F)
## (Intercept) 1 0.005845 55983 7 2304 < 2.2e-16 ***
## PositionGroup 2 0.160517 492 14 4608 < 2.2e-16 ***
## Residuals 2310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(roy)## Analysis of Variance Table
##
## Df Roy approx F num Df den Df Pr(>F)
## (Intercept) 1 170.089 55983 7 2304 < 2.2e-16 ***
## PositionGroup 2 3.114 1025 7 2305 < 2.2e-16 ***
## Residuals 2310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 2d) Univariate ANOVAs
uni_results <- summary.aov(manova_mod)
print(uni_results)## Response PAC :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 12305 6152.6 49.574 < 2.2e-16 ***
## Residuals 2310 286696 124.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response SHO :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 266665 133333 1334.9 < 2.2e-16 ***
## Residuals 2310 230735 100
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PAS :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 32404 16202.0 223.55 < 2.2e-16 ***
## Residuals 2310 167420 72.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DRI :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 42876 21437.9 356.04 < 2.2e-16 ***
## Residuals 2310 139091 60.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DEF :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 395314 197657 1439.3 < 2.2e-16 ***
## Residuals 2310 317222 137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PHY :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 17013 8506.7 101.2 < 2.2e-16 ***
## Residuals 2310 194167 84.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response OVR :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 89 44.528 1.1707 0.3103
## Residuals 2310 87863 38.036The multivariate analysis shows that a player’s primary position is decisively linked to the overall mix of on-field skills. Pillai’s trace is 1.10 (approx. F = 400, p < 0.001), Wilks’ lambda is 0.16 (F = 492, p < 0.001), and Roy’s greatest root is 3.11 (F = 1 025, p < 0.001); each statistic rejects the null hypothesis that attackers, midfielders, and defenders share the same joint mean vector of attributes. Follow-up univariate ANOVAs clarify how that global difference manifests. Defensive Awareness exhibits the greatest separation (F 1 439), with defenders far surpassing the other groups, while Shooting is next (F 1 335), reflecting the scoring edge of attackers. Dribbling (F 356), Passing (F 224), Pace (F 50), and Physicality (F 101) are likewise highly position-dependent, all with p-values well below 0.001. In contrast, overall rating (OVR) shows no significant variation across positions (F 1.17, p = 0.31), suggesting that although EA Sports engineers markedly different skill profiles for each role, it balances the composite rating so that elite attackers, midfielders, and defenders are comparably valued in absolute terms.
colnames(model.matrix(manova_mod))## [1] "(Intercept)" "PositionGroupMidfield" "PositionGroupDefense"print(linearHypothesis(manova_mod, "PositionGroupMidfield = 0"))##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 6201.3146 10201.9012 -8452.8487 -575.47279 -30637.561 2781.5564
## SHO 10201.9012 16783.3427 -13905.9429 -946.72128 -50402.437 4575.9917
## PAS -8452.8487 -13905.9429 11521.8555 784.41180 41761.253 -3791.4664
## DRI -575.4728 -946.7213 784.4118 53.40302 2843.120 -258.1243
## DEF -30637.5610 -50402.4372 41761.2527 2843.12015 151364.703 -13742.2644
## PHY 2781.5564 4575.9917 -3791.4664 -258.12431 -13742.264 1247.6478
## OVR 726.0512 1194.4407 -989.6612 -67.37647 -3587.052 325.6652
## OVR
## PAC 726.05118
## SHO 1194.44067
## PAS -989.66125
## DRI -67.37647
## DEF -3587.05188
## PHY 325.66520
## OVR 85.00622
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 286695.69 82239.04 55609.34 100456.54 -77493.15 -20376.20 35114.36
## SHO 82239.04 230735.37 145487.57 130905.23 28043.71 49085.19 86173.22
## PAS 55609.34 145487.57 167420.25 130274.19 93438.93 46488.48 92781.73
## DRI 100456.54 130905.23 130274.19 139090.63 39432.30 23116.80 80366.35
## DEF -77493.15 28043.71 93438.93 39432.30 317221.99 152202.25 83155.05
## PHY -20376.20 49085.19 46488.48 23116.80 152202.25 194166.82 72189.68
## OVR 35114.36 86173.22 92781.73 80366.35 83155.05 72189.68 87862.77
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.5276222 367.636 7 2304 < 2.22e-16 ***
## Wilks 1 0.4723778 367.636 7 2304 < 2.22e-16 ***
## Hotelling-Lawley 1 1.1169495 367.636 7 2304 < 2.22e-16 ***
## Roy 1 1.1169495 367.636 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# (b) Attack vs Defense
print(linearHypothesis(manova_mod, "PositionGroupDefense = 0"))##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 12280.0597 50478.17 3725.2112 16322.801 -69468.358 -7585.8886
## SHO 50478.1741 207494.60 15312.7807 67096.187 -285555.277 -31182.4058
## PAS 3725.2112 15312.78 1130.0595 4951.595 -21073.538 -2301.2133
## DRI 16322.8007 67096.19 4951.5949 21696.460 -92338.163 -10083.2529
## DEF -69468.3576 -285555.28 -21073.5381 -92338.163 392982.838 42913.4087
## PHY -7585.8886 -31182.41 -2301.2133 -10083.253 42913.409 4686.1096
## OVR 856.0422 3518.83 259.6843 1137.861 -4842.635 -528.8118
## OVR
## PAC 856.04222
## SHO 3518.83049
## PAS 259.68425
## DRI 1137.86146
## DEF -4842.63504
## PHY -528.81184
## OVR 59.67465
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 286695.69 82239.04 55609.34 100456.54 -77493.15 -20376.20 35114.36
## SHO 82239.04 230735.37 145487.57 130905.23 28043.71 49085.19 86173.22
## PAS 55609.34 145487.57 167420.25 130274.19 93438.93 46488.48 92781.73
## DRI 100456.54 130905.23 130274.19 139090.63 39432.30 23116.80 80366.35
## DEF -77493.15 28043.71 93438.93 39432.30 317221.99 152202.25 83155.05
## PHY -20376.20 49085.19 46488.48 23116.80 152202.25 194166.82 72189.68
## OVR 35114.36 86173.22 92781.73 80366.35 83155.05 72189.68 87862.77
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.7478097 975.9942 7 2304 < 2.22e-16 ***
## Wilks 1 0.2521903 975.9942 7 2304 < 2.22e-16 ***
## Hotelling-Lawley 1 2.9652601 975.9942 7 2304 < 2.22e-16 ***
## Roy 1 2.9652601 975.9942 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# (c) Defense vs Midfield
print(linearHypothesis(manova_mod, "PositionGroupDefense - PositionGroupMidfield = 0"))##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 1774.65247 17377.2780 7353.1538 8174.1552 -12969.3640 -5456.8065
## SHO 17377.27796 170157.1402 72001.5887 80040.7797 -126995.1425 -53432.6835
## PAS 7353.15380 72001.5887 30467.3008 33869.0653 -53737.6922 -22609.9128
## DRI 8174.15523 80040.7797 33869.0653 37650.6469 -59737.6649 -25134.3766
## DEF -12969.36401 -126995.1425 -53737.6922 -59737.6649 94781.6012 39878.9686
## PHY -5456.80653 -53432.6835 -22609.9128 -25134.3766 39878.9686 16778.9119
## OVR -70.00827 -685.5163 -290.0746 -322.4623 511.6285 215.2656
## OVR
## PAC -70.008271
## SHO -685.516289
## PAS -290.074586
## DRI -322.462273
## DEF 511.628479
## PHY 215.265577
## OVR 2.761756
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 286695.69 82239.04 55609.34 100456.54 -77493.15 -20376.20 35114.36
## SHO 82239.04 230735.37 145487.57 130905.23 28043.71 49085.19 86173.22
## PAS 55609.34 145487.57 167420.25 130274.19 93438.93 46488.48 92781.73
## DRI 100456.54 130905.23 130274.19 139090.63 39432.30 23116.80 80366.35
## DEF -77493.15 28043.71 93438.93 39432.30 317221.99 152202.25 83155.05
## PHY -20376.20 49085.19 46488.48 23116.80 152202.25 194166.82 72189.68
## OVR 35114.36 86173.22 92781.73 80366.35 83155.05 72189.68 87862.77
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.6052851 504.7321 7 2304 < 2.22e-16 ***
## Wilks 1 0.3947149 504.7321 7 2304 < 2.22e-16 ***
## Hotelling-Lawley 1 1.5334743 504.7321 7 2304 < 2.22e-16 ***
## Roy 1 1.5334743 504.7321 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Step 3b: Univariate contrasts for each attribute
for (v in resp) {
cat("\n=== Univariate contrasts for", v, "===\n")
uni_mod <- lm(as.formula(paste(v, "~ PositionGroup")), data = df)
# Attack vs Midfield & Attack vs Defense
print(linearHypothesis(uni_mod, c("PositionGroupMidfield = 0",
"PositionGroupDefense = 0")))
}##
## === Univariate contrasts for PAC ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: PAC ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 299001
## 2 2310 286696 2 12305 49.574 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for SHO ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: SHO ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 497401
## 2 2310 230735 2 266665 1334.9 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for PAS ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: PAS ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 199824
## 2 2310 167420 2 32404 223.55 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for DRI ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: DRI ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 181966
## 2 2310 139091 2 42876 356.04 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for DEF ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: DEF ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 712536
## 2 2310 317222 2 395314 1439.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for PHY ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: PHY ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 211180
## 2 2310 194167 2 17013 101.2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## === Univariate contrasts for OVR ===
## Linear hypothesis test
##
## Hypothesis:
## PositionGroupMidfield = 0
## PositionGroupDefense = 0
##
## Model 1: restricted model
## Model 2: OVR ~ PositionGroup
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2312 87952
## 2 2310 87863 2 89.055 1.1707 0.3103# –– 4. Add a continuous covariate (e.g. Age) ––# 4a)
# Scatterplots to check linearity
pairs(df[c(resp, "Age")], pch = 20, cex = 0.5, main = "Attribute vs. Age scatterplots")
options(contrasts = c("contr.sum", "contr.poly"))
# construct: cbind(PAC,SHO,...,OVR) ~ PositionGroup + Age
mv_age_form <- as.formula(paste0("cbind(", paste(resp, collapse = ","),
") ~ PositionGroup + Age"))
# fit the model
manova_cov <- lm(mv_age_form, data = df)
# 4b.1) Sequential (Type I) multivariate tests (Pillai's
# trace):
print(anova(manova_cov))## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99504 65993 7 2303 < 2.2e-16 ***
## PositionGroup 2 1.10808 409 14 4608 < 2.2e-16 ***
## Age 1 0.41200 231 7 2303 < 2.2e-16 ***
## Residuals 2309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # DISCUSS: Pillai's trace for PositionGroup (first) and
# Age (second)
# 4b.2) Sequential (Type I) univariate ANOVAs:
print(summary.aov(manova_cov))## Response PAC :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 12305 6152.6 51.202 < 2.2e-16 ***
## Age 1 9238 9238.0 76.878 < 2.2e-16 ***
## Residuals 2309 277458 120.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response SHO :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 266665 133333 1640.93 < 2.2e-16 ***
## Age 1 43119 43119 530.66 < 2.2e-16 ***
## Residuals 2309 187617 81
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PAS :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 32404 16202 282.28 < 2.2e-16 ***
## Age 1 34891 34891 607.89 < 2.2e-16 ***
## Residuals 2309 132529 57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DRI :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 42876 21437.9 393.07 < 2.2e-16 ***
## Age 1 13160 13159.9 241.29 < 2.2e-16 ***
## Residuals 2309 125931 54.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DEF :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 395314 197657 1684.18 < 2.2e-16 ***
## Age 1 46235 46235 393.95 < 2.2e-16 ***
## Residuals 2309 270987 117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PHY :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 17013 8506.7 120.72 < 2.2e-16 ***
## Age 1 31459 31458.7 446.43 < 2.2e-16 ***
## Residuals 2309 162708 70.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response OVR :
## Df Sum Sq Mean Sq F value Pr(>F)
## PositionGroup 2 89 44.5 1.5541 0.2116
## Age 1 21707 21706.8 757.6204 <2e-16 ***
## Residuals 2309 66156 28.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# # DISCUSS: Which attrs show a significant Age effect in
# the sequential ANOVAs?
# 4b.3) Type II sums of squares (multivariate & univariate)
library(car)
print(summary(Anova(manova_cov, type = "II"), multivariate = TRUE))##
## Type II MANOVA Tests:
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 277457.710 102197.25 73562.65 111482.472 -56826.36 -3328.765 49275.12
## SHO 102197.249 187616.60 106700.29 107084.222 -16605.93 12255.005 55579.56
## PAS 73562.646 106700.29 132529.35 108846.121 53274.57 13358.073 65261.36
## DRI 111482.472 107084.22 108846.12 125930.695 14765.57 2769.934 63464.85
## DEF -56826.356 -16605.93 53274.57 14765.571 270987.14 114064.470 51475.21
## PHY -3328.765 12255.00 13358.07 2769.934 114064.47 162708.079 46057.90
## OVR 49275.125 55579.56 65261.36 63464.852 51475.21 46057.900 66155.93
##
## ------------------------------------------
##
## Term: PositionGroup
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 11429.099 45356.737 1634.576 13851.459 -64772.204 -4903.987
## SHO 45356.737 284288.231 70522.222 111575.189 -292534.744 -51799.148
## PAS 1634.576 70522.222 39552.804 36737.805 -31051.692 -20557.268
## DRI 13851.459 111575.189 36737.805 47511.136 -97760.410 -23495.368
## DEF -64772.204 -292534.744 -31051.692 -97760.410 379157.433 38795.247
## PHY -4903.987 -51799.148 -20557.268 -23495.368 38795.247 12131.328
## OVR 1463.869 9388.982 2407.293 3717.034 -9514.143 -1738.060
## OVR
## PAC 1463.8686
## SHO 9388.9822
## PAS 2407.2931
## DRI 3717.0338
## DEF -9514.1425
## PHY -1738.0598
## OVR 310.3604
##
## Multivariate Tests: PositionGroup
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 1.122422 420.9736 14 4608 < 2.22e-16 ***
## Wilks 2 0.152122 514.5289 14 4606 < 2.22e-16 ***
## Hotelling-Lawley 2 3.768912 619.7168 14 4604 < 2.22e-16 ***
## Roy 2 3.205975 1055.2238 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Age
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 9237.977 -19958.21 -17953.31 -11025.93 -20666.80 -17047.44 -14160.76
## SHO -19958.213 43118.78 38787.28 23821.01 44649.64 36830.18 30593.66
## PAS -17953.310 38787.28 34890.90 21428.07 40164.36 33130.40 27520.38
## DRI -11025.932 23821.01 21428.07 13159.94 24666.73 20346.87 16901.50
## DEF -20666.798 44649.64 40164.36 24666.73 46234.86 38137.78 31679.84
## PHY -17047.437 36830.18 33130.40 20346.87 38137.78 31458.74 26131.77
## OVR -14160.764 30593.66 27520.38 16901.50 31679.84 26131.77 21706.84
##
## Multivariate Tests: Age
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4119986 230.5225 7 2303 < 2.22e-16 ***
## Wilks 1 0.5880014 230.5225 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## Roy 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(summary(Anova(manova_cov, type = "II"), univariate = TRUE))##
## Type II MANOVA Tests:
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 277457.710 102197.25 73562.65 111482.472 -56826.36 -3328.765 49275.12
## SHO 102197.249 187616.60 106700.29 107084.222 -16605.93 12255.005 55579.56
## PAS 73562.646 106700.29 132529.35 108846.121 53274.57 13358.073 65261.36
## DRI 111482.472 107084.22 108846.12 125930.695 14765.57 2769.934 63464.85
## DEF -56826.356 -16605.93 53274.57 14765.571 270987.14 114064.470 51475.21
## PHY -3328.765 12255.00 13358.07 2769.934 114064.47 162708.079 46057.90
## OVR 49275.125 55579.56 65261.36 63464.852 51475.21 46057.900 66155.93
##
## ------------------------------------------
##
## Term: PositionGroup
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 11429.099 45356.737 1634.576 13851.459 -64772.204 -4903.987
## SHO 45356.737 284288.231 70522.222 111575.189 -292534.744 -51799.148
## PAS 1634.576 70522.222 39552.804 36737.805 -31051.692 -20557.268
## DRI 13851.459 111575.189 36737.805 47511.136 -97760.410 -23495.368
## DEF -64772.204 -292534.744 -31051.692 -97760.410 379157.433 38795.247
## PHY -4903.987 -51799.148 -20557.268 -23495.368 38795.247 12131.328
## OVR 1463.869 9388.982 2407.293 3717.034 -9514.143 -1738.060
## OVR
## PAC 1463.8686
## SHO 9388.9822
## PAS 2407.2931
## DRI 3717.0338
## DEF -9514.1425
## PHY -1738.0598
## OVR 310.3604
##
## Multivariate Tests: PositionGroup
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 1.122422 420.9736 14 4608 < 2.22e-16 ***
## Wilks 2 0.152122 514.5289 14 4606 < 2.22e-16 ***
## Hotelling-Lawley 2 3.768912 619.7168 14 4604 < 2.22e-16 ***
## Roy 2 3.205975 1055.2238 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Age
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 9237.977 -19958.21 -17953.31 -11025.93 -20666.80 -17047.44 -14160.76
## SHO -19958.213 43118.78 38787.28 23821.01 44649.64 36830.18 30593.66
## PAS -17953.310 38787.28 34890.90 21428.07 40164.36 33130.40 27520.38
## DRI -11025.932 23821.01 21428.07 13159.94 24666.73 20346.87 16901.50
## DEF -20666.798 44649.64 40164.36 24666.73 46234.86 38137.78 31679.84
## PHY -17047.437 36830.18 33130.40 20346.87 38137.78 31458.74 26131.77
## OVR -14160.764 30593.66 27520.38 16901.50 31679.84 26131.77 21706.84
##
## Multivariate Tests: Age
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4119986 230.5225 7 2303 < 2.22e-16 ***
## Wilks 1 0.5880014 230.5225 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## Roy 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type II Sums of Squares
## df PAC SHO PAS DRI DEF PHY OVR
## PositionGroup 2 11429 284288 39553 47511 379157 12131 310.36
## Age 1 9238 43119 34891 13160 46235 31459 21706.84
## residuals 2309 277458 187617 132529 125931 270987 162708 66155.93
##
## F-tests
## PAC SHO PAS DRI DEF PHY OVR
## PositionGroup 47.56 3498.74 344.56 871.14 1615.34 172.16 5.42
## Age 76.88 265.33 607.89 120.65 393.95 223.22 757.62
##
## p-values
## PAC SHO PAS DRI DEF PHY
## PositionGroup < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16
## Age < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16
## OVR
## PositionGroup 0.0045008
## Age < 2.22e-16# # DISCUSS: Do PositionGroup and Age both enter
# significantly under Type II?
# 4b.4) (If you want) Type III sums of squares Make sure
# you’ve set contrasts correctly before this:
# options(contrasts=c('contr.sum','contr.poly'))
print(summary(Anova(manova_cov, type = "III"), multivariate = TRUE))##
## Type III MANOVA Tests:
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 277457.710 102197.25 73562.65 111482.472 -56826.36 -3328.765 49275.12
## SHO 102197.249 187616.60 106700.29 107084.222 -16605.93 12255.005 55579.56
## PAS 73562.646 106700.29 132529.35 108846.121 53274.57 13358.073 65261.36
## DRI 111482.472 107084.22 108846.12 125930.695 14765.57 2769.934 63464.85
## DEF -56826.356 -16605.93 53274.57 14765.571 270987.14 114064.470 51475.21
## PHY -3328.765 12255.00 13358.07 2769.934 114064.47 162708.079 46057.90
## OVR 49275.125 55579.56 65261.36 63464.852 51475.21 46057.900 66155.93
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 431548.5 183304.80 215532.44 291683.63 145835.82 241593.67 285416.23
## SHO 183304.8 77860.65 91549.69 123895.71 61945.31 102619.47 121233.56
## PAS 215532.4 91549.69 107645.45 145678.37 72836.19 120661.46 142548.18
## DRI 291683.6 123895.71 145678.37 197148.96 98570.43 163293.16 192912.82
## DEF 145835.8 61945.31 72836.19 98570.43 49283.19 81643.22 96452.45
## PHY 241593.7 102619.47 120661.46 163293.16 81643.22 135251.31 159784.48
## OVR 285416.2 121233.56 142548.18 192912.82 96452.45 159784.48 188767.71
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.789763 1235.9 7 2303 < 2.22e-16 ***
## Wilks 1 0.210237 1235.9 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 3.756534 1235.9 7 2303 < 2.22e-16 ***
## Roy 1 3.756534 1235.9 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: PositionGroup
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 11429.099 45356.737 1634.576 13851.459 -64772.204 -4903.987
## SHO 45356.737 284288.231 70522.222 111575.189 -292534.744 -51799.148
## PAS 1634.576 70522.222 39552.804 36737.805 -31051.692 -20557.268
## DRI 13851.459 111575.189 36737.805 47511.136 -97760.410 -23495.368
## DEF -64772.204 -292534.744 -31051.692 -97760.410 379157.433 38795.247
## PHY -4903.987 -51799.148 -20557.268 -23495.368 38795.247 12131.328
## OVR 1463.869 9388.982 2407.293 3717.034 -9514.143 -1738.060
## OVR
## PAC 1463.8686
## SHO 9388.9822
## PAS 2407.2931
## DRI 3717.0338
## DEF -9514.1425
## PHY -1738.0598
## OVR 310.3604
##
## Multivariate Tests: PositionGroup
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 1.122422 420.9736 14 4608 < 2.22e-16 ***
## Wilks 2 0.152122 514.5289 14 4606 < 2.22e-16 ***
## Hotelling-Lawley 2 3.768912 619.7168 14 4604 < 2.22e-16 ***
## Roy 2 3.205975 1055.2238 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Age
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 9237.977 -19958.21 -17953.31 -11025.93 -20666.80 -17047.44 -14160.76
## SHO -19958.213 43118.78 38787.28 23821.01 44649.64 36830.18 30593.66
## PAS -17953.310 38787.28 34890.90 21428.07 40164.36 33130.40 27520.38
## DRI -11025.932 23821.01 21428.07 13159.94 24666.73 20346.87 16901.50
## DEF -20666.798 44649.64 40164.36 24666.73 46234.86 38137.78 31679.84
## PHY -17047.437 36830.18 33130.40 20346.87 38137.78 31458.74 26131.77
## OVR -14160.764 30593.66 27520.38 16901.50 31679.84 26131.77 21706.84
##
## Multivariate Tests: Age
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4119986 230.5225 7 2303 < 2.22e-16 ***
## Wilks 1 0.5880014 230.5225 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## Roy 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(summary(Anova(manova_cov, type = "III"), univariate = TRUE))##
## Type III MANOVA Tests:
##
## Sum of squares and products for error:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 277457.710 102197.25 73562.65 111482.472 -56826.36 -3328.765 49275.12
## SHO 102197.249 187616.60 106700.29 107084.222 -16605.93 12255.005 55579.56
## PAS 73562.646 106700.29 132529.35 108846.121 53274.57 13358.073 65261.36
## DRI 111482.472 107084.22 108846.12 125930.695 14765.57 2769.934 63464.85
## DEF -56826.356 -16605.93 53274.57 14765.571 270987.14 114064.470 51475.21
## PHY -3328.765 12255.00 13358.07 2769.934 114064.47 162708.079 46057.90
## OVR 49275.125 55579.56 65261.36 63464.852 51475.21 46057.900 66155.93
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 431548.5 183304.80 215532.44 291683.63 145835.82 241593.67 285416.23
## SHO 183304.8 77860.65 91549.69 123895.71 61945.31 102619.47 121233.56
## PAS 215532.4 91549.69 107645.45 145678.37 72836.19 120661.46 142548.18
## DRI 291683.6 123895.71 145678.37 197148.96 98570.43 163293.16 192912.82
## DEF 145835.8 61945.31 72836.19 98570.43 49283.19 81643.22 96452.45
## PHY 241593.7 102619.47 120661.46 163293.16 81643.22 135251.31 159784.48
## OVR 285416.2 121233.56 142548.18 192912.82 96452.45 159784.48 188767.71
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.789763 1235.9 7 2303 < 2.22e-16 ***
## Wilks 1 0.210237 1235.9 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 3.756534 1235.9 7 2303 < 2.22e-16 ***
## Roy 1 3.756534 1235.9 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: PositionGroup
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY
## PAC 11429.099 45356.737 1634.576 13851.459 -64772.204 -4903.987
## SHO 45356.737 284288.231 70522.222 111575.189 -292534.744 -51799.148
## PAS 1634.576 70522.222 39552.804 36737.805 -31051.692 -20557.268
## DRI 13851.459 111575.189 36737.805 47511.136 -97760.410 -23495.368
## DEF -64772.204 -292534.744 -31051.692 -97760.410 379157.433 38795.247
## PHY -4903.987 -51799.148 -20557.268 -23495.368 38795.247 12131.328
## OVR 1463.869 9388.982 2407.293 3717.034 -9514.143 -1738.060
## OVR
## PAC 1463.8686
## SHO 9388.9822
## PAS 2407.2931
## DRI 3717.0338
## DEF -9514.1425
## PHY -1738.0598
## OVR 310.3604
##
## Multivariate Tests: PositionGroup
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 1.122422 420.9736 14 4608 < 2.22e-16 ***
## Wilks 2 0.152122 514.5289 14 4606 < 2.22e-16 ***
## Hotelling-Lawley 2 3.768912 619.7168 14 4604 < 2.22e-16 ***
## Roy 2 3.205975 1055.2238 7 2304 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Age
##
## Sum of squares and products for the hypothesis:
## PAC SHO PAS DRI DEF PHY OVR
## PAC 9237.977 -19958.21 -17953.31 -11025.93 -20666.80 -17047.44 -14160.76
## SHO -19958.213 43118.78 38787.28 23821.01 44649.64 36830.18 30593.66
## PAS -17953.310 38787.28 34890.90 21428.07 40164.36 33130.40 27520.38
## DRI -11025.932 23821.01 21428.07 13159.94 24666.73 20346.87 16901.50
## DEF -20666.798 44649.64 40164.36 24666.73 46234.86 38137.78 31679.84
## PHY -17047.437 36830.18 33130.40 20346.87 38137.78 31458.74 26131.77
## OVR -14160.764 30593.66 27520.38 16901.50 31679.84 26131.77 21706.84
##
## Multivariate Tests: Age
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4119986 230.5225 7 2303 < 2.22e-16 ***
## Wilks 1 0.5880014 230.5225 7 2303 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## Roy 1 0.7006762 230.5225 7 2303 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type III Sums of Squares
## df PAC SHO PAS DRI DEF PHY OVR
## (Intercept) 1 431549 77861 107645 197149 49283 135251 188767.71
## PositionGroup 2 11429 284288 39553 47511 379157 12131 310.36
## Age 1 9238 43119 34891 13160 46235 31459 21706.84
## residuals 2309 277458 187617 132529 125931 270987 162708 66155.93
##
## F-tests
## PAC SHO PAS DRI DEF PHY OVR
## (Intercept) 3591.34 479.12 1875.46 3614.82 209.96 1919.36 6588.44
## PositionGroup 47.56 3498.74 689.11 435.57 3230.69 172.16 5.42
## Age 76.88 530.66 303.94 241.29 393.95 223.22 757.62
##
## p-values
## PAC SHO PAS DRI DEF PHY
## (Intercept) < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16
## PositionGroup < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16
## Age < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16 < 2.22e-16
## OVR
## (Intercept) < 2.22e-16
## PositionGroup 0.0045008
## Age < 2.22e-16# –– 5. Check model assumptions ––# Chi‐square Q–Q of
# residuals (overall multivariate residuals)
cqplot(manova_cov$residuals, label = "Residuals (Position MANOVA)")
we will attempt box cox transformation
## bocxcox
pt <- powerTransform(as.matrix(df[resp]))
summary(pt)## bcPower Transformations to Multinormality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## PAC 2.2493 2.25 2.0611 2.4374
## SHO 2.1594 2.16 2.0334 2.2853
## PAS 2.0349 2.00 1.8781 2.1916
## DRI 2.8890 2.89 2.7275 3.0505
## DEF 2.1504 2.15 2.0253 2.2755
## PHY 2.7714 2.77 2.5568 2.9861
## OVR 1.8479 2.00 1.5858 2.1100
##
## Likelihood ratio test that transformation parameters are equal to 0
## (all log transformations)
## LRT df pval
## LR test, lambda = (0 0 0 0 0 0 0) 4609.838 7 < 2.22e-16
##
## Likelihood ratio test that no transformations are needed
## LRT df pval
## LR test, lambda = (1 1 1 1 1 1 1) 1457.332 7 < 2.22e-16# you’ll get one λ per response (or a common λ if you
# specify)
lambda_vec <- coef(pt, round = 2)
print(lambda_vec)## PAC SHO PAS DRI DEF PHY OVR
## 2.249272 2.159376 2.034865 2.888965 2.150368 2.771437 1.847903# to apply the transformation in one step:
df_bc <- as.data.frame(bcPower(df[resp], lambda = lambda_vec))
names(df_bc) <- paste0(resp, "_bc")
df2 <- bind_cols(df, df_bc)
# now fit your MANOVA on the '_bc' columns
mv_bc_form <- as.formula(paste0("cbind(", paste(paste0(resp,
"_bc"), collapse = ","), ") ~ PositionGroup + Age"))
manova_bc <- lm(mv_bc_form, data = df2)
cqplot(manova_bc$residuals, label = "Residuals (Position MANOVA, Box-Cox transformed)")
# If severe non-normality appears, consider
# transformations: e.g. df$PAC_t <- log(df$PAC + 1) and
# refit the model
# –– 6. MRPP on the same grouping ––# Multi‐response
# permutation procedure
resp_mat <- as.matrix(df[resp])
mrpp_res <- mrpp(resp_mat, df$Position, permutations = 999)
print(mrpp_res)##
## Call:
## mrpp(dat = resp_mat, grouping = df$Position, permutations = 999)
##
## Dissimilarity index: euclidean
## Weights for groups: n
##
## Class means and counts:
##
## CAM CB CDM CM LB LM LW RB RM RW ST
## delta 31.06 29.42 25.72 26.31 26.26 28.46 29.66 26.23 29.58 28.92 27.69
## n 153 478 210 312 192 155 52 210 159 53 339
##
## Chance corrected within-group agreement A: 0.3105
## Based on observed delta 27.91 and expected delta 40.48
##
## Significance of delta: 0.001
## Permutation: free
## Number of permutations: 999# End of script