MvA_HW2_LPC
LPC
2023-01-11
Luka Pascal Cavazza, 11. January 2023
IMB; Multivariate Analysis;
HOMEWORK 2
Post-Coma Recovery of IQ
For my second homework I had to choose a different dataset, since the first one gave me a lot of headaches. I figured it might be interesting to check the built in datasets that R and its packages provide. To see all of what is available I used “data(package = .packages(all.available = TRUE))”. After goining through it I came across this in the library “carData” and I found it very interesting as it might show us if there can be significant recovery of brain functions (IQ to be specific) after patients woke up from comas. I don’t know what it is with these “medical” datasets that intrigue me, but this one looks much better.
So my unit of observation are the people that got injured and were consequently in a coma (and woke up), the sample size is 331 observations.
Description of variables:
- id: patients ID number
- days: number of days post coma at which IQs were measured
- duration: duration of the coma in days
- sex: gender
- age: in years at the time of injury
- piq: performance (i.e., mathematical) IQ
- viq: verbal IQ
Source of the dataset is as mentioned before the carData library.
The basic idea is to see whether duration of coma affects IQ and if it can recover over time (if there’ll be more than one measurement per patient).
So let’s see what we are dealing with.
library(carData)
data(Wong)
head(Wong)## id days duration sex age piq viq
## 1 3358 30 4 Male 20.67077 87 89
## 2 3535 16 17 Male 55.28816 95 77
## 3 3547 40 1 Male 55.91513 95 116
## 4 3592 13 10 Male 61.66461 59 73
## 5 3728 19 6 Male 30.12731 67 73
## 6 3790 13 3 Male 57.06229 76 69The first thing that I notice is how age is written. It’s obvious that it was calculated automatically but I’ll change it so it makes more sense.
Wong$age <- round(Wong$age, 0) #rounding the age column to the nearest number
data <- WongI just got an idea. Let’s see if there are any duplicates within ID and then we can see if there is any progress between 1st and 2nd test score.
#install.packages("misty")
library(misty)## |-------------------------------------|
## | misty 0.4.7 (2023-01-06) |
## | Miscellaneous Functions T. Yanagida |
## |-------------------------------------|data2M <- df.duplicated(data, id, first = TRUE, keep.all = TRUE, #extracting the duplicate values in ID in a new data frame
from.last = FALSE, keep.row.names = TRUE,
check = TRUE)
dataU <- df.unique(data, id, keep.all = TRUE, #extracting unique data in a new data frame
from.last = FALSE, keep.row.names = TRUE,
check = TRUE)
sum(duplicated(data2M$id)) #checking if it works, should be 224 observations## [1] 131This is a bit weird, but I believe the problem is some people were tested more than twice and some just once. I hope this will not be too big of a problem.
It would make sense to make a categorical variable with 3 levels, if people were in a coma up to 10 days I will assign them a value of 1, for 10 to 50 days a value of 2 and above 50 a value of 3
createDurationFactor <- function(duration) { #creating a duration factor so I have 3 levels
return (ifelse (duration < 10, 1, ifelse(duration > 50, 3, 2)))
}
data$DurationF <- createDurationFactor(data$duration)Now I can start with the hypothesis testing.
Hypothesis about the arithmetic mean
Let’s see if the variables are normally distributed.
ggplot(data2M, aes(x = duration)) +
geom_histogram(binwidth = 5, colour = "black") +
ylab("Frequency") +
xlab("Duration")
It’s asymetrical to the right to we should use non-parametric tests.
mean(data2M$duration)## [1] 13.79018sd(data2M$duration)## [1] 22.23046nrow(data2M)## [1] 224This will be too tricky so I’ll move on.
Wilcoxon signed rank test
Lets compare performance and verbal IQ.
dataU$Difference <- dataU$piq - dataU$viq #calculating differences
library(ggplot2)
ggplot(dataU, aes(x = Difference)) +
geom_histogram(binwidth = 4, color = "black") +
xlab("Differences")
It’s not normally distributed so I’ll use non parametric tests, I should
use the Wilcoxon Signed Rank Test. H0: distributions locations are the
same H1: distribution locations are different
library(effectsize)
wilcox.test(dataU$piq, dataU$viq,
paired = TRUE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided")##
## Wilcoxon signed rank test
##
## data: dataU$piq and dataU$viq
## V = 2712, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0effectsize(wilcox.test(data$piq, data$viq,
paired = TRUE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided"))## r (rank biserial) | 95% CI
## ----------------------------------
## -0.63 | [-0.70, -0.55]Based on the p-value (alfa=5%) we can reject the null hypothesis.
Let’s see the effect size:
interpret_rank_biserial(0.63)## [1] "very large"
## (Rules: funder2019)The effect size is very large.
Wilcoxon Rank Sum Test
data$IQ <- (data$piq + data$viq)/2 #I just made a new column with the average IQ score
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:effectsize':
##
## phi## The following object is masked from 'package:car':
##
## logit## The following objects are masked from 'package:ggplot2':
##
## %+%, alphadescribeBy(data$IQ, g = data$sex)##
## Descriptive statistics by group
## group: Female
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 71 91.77 14.76 91 91.25 15.57 62 125 63 0.27 -0.6 1.75
## ------------------------------------------------------------
## group: Male
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 260 91.12 12.85 90 90.6 12.6 61 126.5 65.5 0.34 -0.31 0.8I’m comparing 2 arithmetic means and they look the same. We can also see that it’s mostly men that suffer these injuries.
ggplot(data, aes(x = IQ)) +
geom_histogram(binwidth = 5, colour="gray") +
facet_wrap(~sex, ncol = 1) +
ylab("Frequency")
I would say that the Male is asymetrical so I will use non-parametrical
tests.
H0: distribution locations are the same H1: distribution locations are different
wilcox.test(data$IQ ~ data$sex,
paired = FALSE, #for independent variables
correct = FALSE,
exact = FALSE,
alternative = "two.sided")##
## Wilcoxon rank sum test
##
## data: data$IQ by data$sex
## W = 9357.5, p-value = 0.8584
## alternative hypothesis: true location shift is not equal to 0effectsize(wilcox.test(data$IQ ~ data$sex,
paired = FALSE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided"))## r (rank biserial) | 95% CI
## ---------------------------------
## 0.01 | [-0.14, 0.16]Based on the p-value we can not reject null hypothesis, meaning Male & Female IQ scores behave the same (are in the same range).
Let’s see the effect size:
interpret_rank_biserial(0.01)## [1] "tiny"
## (Rules: funder2019)The effect size is tiny.
Kruskal-Wallis Rank Sum test
data$ComaFactor <- factor(data$DurationF,
levels = c(1, 2, 3),
labels = c("Less than 10", "10 to 50", "50 or more"))psych::describeBy(x = data$IQ, group = data$DurationF)##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 200 93 13.87 92 92.6 13.34 62 126.5 64.5 0.24 -0.54 0.98
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 114 88.73 11.99 86.75 88.25 11.49 61 123 62 0.41 -0.05 1.12
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 17 87.74 10.93 92 88 11.12 70 101.5 31.5 -0.36 -1.53 2.65I’m comparing means, let’s check weather variances are the same in the popoulation (sd).
Levene test: homogenity of variances (if this one is not met (sig. p value) then you should use heteroskedastic anova)
H0: variances of happiness are the same in all groups H1: not the same
library(car)
leveneTest(data$IQ, group = data$DurationF)## Warning in leveneTest.default(data$IQ, group = data$DurationF): data$DurationF
## coerced to factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.0823 0.1263
## 328Because of the p-value we shouldn’t reject H0 We can assume homogenity.
Let’s check for normality (it must hold for each subgroup).
library(dplyr)
library(rstatix)##
## Attaching package: 'rstatix'## The following objects are masked from 'package:effectsize':
##
## cohens_d, eta_squared## The following object is masked from 'package:stats':
##
## filterdata %>%
group_by(DurationF) %>% #checking by groups
shapiro_test(IQ)## # A tibble: 3 × 4
## DurationF variable statistic p
## <dbl> <chr> <dbl> <dbl>
## 1 1 IQ 0.987 0.0553
## 2 2 IQ 0.982 0.140
## 3 3 IQ 0.893 0.0516H0: IQ is normally distributed H1: not normally
Done the same for all of the subgroups
All p-values are above 5% so everything is OK.
ANOVA_results <- aov(IQ ~ DurationF, #analysis of variance
data = data)
summary(ANOVA_results)## Df Sum Sq Mean Sq F value Pr(>F)
## DurationF 1 1440 1440 8.375 0.00406 **
## Residuals 329 56584 172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1let’s see the effect size:
library(effectsize)
eta_squared(ANOVA_results)## DurationF
## 0.02482393interpret_eta_squared(0.02, rules = "cohen1992")## [1] "small"
## (Rules: cohen1992)There is a small difference between groups that have been in a coma.
This has been a tricky homework, I think everything has to be done a couple of times to really understand it. To be honest I’m a bit confused and probably overcomplicated some things again. I hope I didn’t do any mistakes because of this.
Again I wanted to do too much and didn’t have time to check the things I really wanted to know.