https://web.stanford.edu/class/psych252/tutorials/index.html
Anmerkung: Das ist ziemlich viel Spielerei: Die Aufgabe braucht nur eigentlich nur ~ 8 Zeilen Code. Daten eingeben – Jede Spalte abtippen und in einen Vektor speichern, Vektoren dann zu einem Dataframe zusammenfassen. Wenns mehr Daten wären würde ichs extern speichern (.csv):
library(gdata)
library(tidyr)
library(dplyr)
Attacken <- c(0,1,2,3,4,5,6)
Kognitiv <- c(29,29,16,3,0,2,0)
Interozeptiv <- c(12,23,19,8,4,1,1)
df <- data.frame(Attacken,Kognitiv,Interozeptiv)
attacks <- rep(Attacken,Kognitiv)
cog <- data.frame(cbind(therapy = "Kognitiv", attacks))
attacks <- rep(Attacken,Interozeptiv)
int <- data.frame(cbind(therapy = "Interozeptiv", attacks))
output <- bind_rows(cog,int)
library(ggplot2)
output$therapy <- as.factor(output$therapy)
p <- ggplot(output, aes(x=therapy, y=attacks))
p + #geom_jitter(shape=16, position=position_jitter(0.2))
geom_dotplot(binaxis='y', stackdir='center', dotsize=.4)
Success (0) und Failure(1) codieren, Wilcox-U-Test durchführen:
output$attacks[output$attacks<2]<-0
output$attacks[output$attacks>=2]<-1
wilcox.test(as.numeric(attacks) ~ therapy, data=output)
##
## Wilcoxon rank sum test with continuity correction
##
## data: as.numeric(attacks) by therapy
## W = 3275.5, p-value = 0.006136
## alternative hypothesis: true location shift is not equal to 0
\(\chi^2\)-Test durchführen:
mytable <- table(output$attacks, output$therapy)
print(mytable)
##
## Interozeptiv Kognitiv
## 0 35 58
## 1 33 21
chisq.test(output$attacks, output$therapy, correct=FALSE)
##
## Pearson's Chi-squared test
##
## data: output$attacks and output$therapy
## X-squared = 7.5741, df = 1, p-value = 0.005921
Und weils so schön war, \(\chi^2\) jetzt nochmal händisch. Über jede Zeile und jede Spalte des Dataframes iterieren und den Wert ausrechnen:
df <-as.data.frame.matrix(mytable)
df$Attacken <- NULL
count <- 0
for (row in 1:nrow(df)) {
for (column in names(df)){
#print(df[row,column])
expected <- rowSums(df)[row]*colSums(df)[column]/sum(rowSums(df))
count <- count + ((df[row,column]-expected)^2/expected)
}
}
print(as.numeric(count))
## [1] 7.574121
degrees <- ((dim(df)[1]-1)*(dim(df)[2]-1))
print(degrees)
## [1] 1
Und jetzt bestimmen wir noch den p-Wert, indem wir die \(\chi^2\)-Funktion integrieren:
x <- seq(3, 20, length=1000)
hx <- dchisq(x,df=degrees)
integrand <- function(x) {dchisq(x,df=degrees)}
integrate(integrand, lower = count, upper = Inf)
## 0.005921195 with absolute error < 4.9e-05
plot(x, hx, type="l", xlab="x", ylab="y", main="chi^2 -Verteilung",col="#8441b4", lwd=2)
x <- seq(count, 20, length=1000)
hx <- dchisq(x,df=degrees)
polygon(c(count,x,20),c(0,hx,0),col=rgb(((100*132)/255*0.01), ((100*65)/255*0.01), ((100*180)/255*0.01),0.5),border=rgb(((100*132)/255*0.01), ((100*65)/255*0.01), ((100*180)/255*0.01),0.5))