stp4/stp25stat
Functions for statistical computations

Package index
Search the stp4/stp25stat package
Vignettes
Functions
572
Source code
93
Man pages
29
Home
/
GitHub
/
stp4/stp25stat
/

In stp4/stp25stat: Functions for statistical computations 

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Einleitung

Die von mir verwendeten statistischen Methoden basieren vor allem auf den Empfehlungen von [@Bortz2006] sowie [@Sachs2006]. Die Tabellen und Grafiken werden mit der Software R [@R-base] im APA-Style Format [@APA] erstellt. Das Literaturverzeichnis erstelle ich mit dem Literaturgenerator und dem Wikipedia BibTeX Generator sowie DOI zo BibTeX Generator.

Link: - Grafiken - Grafiken2 - Beispiele - Kano-Analyse - Regressionsanalyse 

 #getwd()
library(stpvers)
require(broom)
#require(nlme)

set.seed(1)
Projekt("", "Introduction" )
#set_my_options(output="html")

hkarz$Lai<- factor(hkarz$lai, 0:1, c("negativ", "positiv"))
hkarz<- upData2(hkarz, c(gruppe="Gruppe"
                      ,tzell="T-Zelltypisierung"
                      ,Lai="LAI-Test"))

Methoden

Deskriptive Auswertungen

Je nach Kontext wird der Mittelwert, der Median oder die Häufigkeiten berechnet, dabei wird die Schreibweise M(SD) oder Median[unteres Quartil, oberes Quartil], verwendet. Bei Faktoren (Nominal- Skala) wird der Anteil in Prozent angegeben, dabei wird die Schreibweise Prozent(Anzahl) verwendet. Wenn ein Signifikanz Test angegeben ist handelt es sich um einen nicht-parametrischen Test der abhängig vom Skalenniveau entweder ein Chi-Quadrat-test bei Nominal-Skala oder Wilcoxon bzw. Kruskal-Wallis-Test bei ordinalen und metrischen-Skalen ist. (Es werden auf Empfehlungen des APA-Styles die F-Werte mit ausgegeben) Die nicht Parametrischen Tests sind mit Hilfe von Hmisc::spearman2 bzw. stats::chisq.test erstellt.

Hmisc::spearman2[14] Wilcoxon, or Kruskal-Wallis tests are approximated by using the t or F distributions. spearman2 computes the square of Spearman's rho rank correlation and a generalization of it in which x can relate non-monotonically to y. This is done by computing the Spearman multiple rho-squared between (rank(x), rank(x)^2) and y. When x is categorical, a different kind of Spearman correlation used in the Kruskal-Wallis test is computed (and spearman2 can do the Kruskal-Wallis test). This is done by computing the ordinary multiple R^2 between k-1 dummy variables and rank(y), where x has k categories. x can also be a formula, in which case each predictor is correlated separately with y, using non-missing observations for that predictor. biVar is used to do the looping and bookkeeping. By default the plot shows the adjusted rho^2, using the same formula used for the ordinary adjusted R^2. The F test uses the unadjusted R2.

Anmerkung zu einfachen Auswertungen

Wilcoxon-Test vs. Mann-Whitney-U-Test

Beide Tests sind die gleichen Testverfahren, die genaue Bezeichnung ist Wilcoxon-Mann-Whitney-Test und auch der H-Test ist mathematisch dasselbe Verfahren [@Bortz12010] (Seite 140 und 203). (Die Verwechslung rührt wahrscheinlich daher, dass noch der Wilcoxon-Vorzeichenrangtest als Verfahren existiert.)

Testen auf Normalverteilung

Die Normalverteilung ist keine zwingende Voraussetzung für den T-Test, die ANOVA oder die Regressionsanalyse. Bortz schreibt hier - ab ca. 30 Untersuchungsobjekten erübrigt sich die Forderung nach einer Normalverteilung [@Bortz2006] (Seite 218).
Normal verteilt sollen die Residuen sein! Bei den Daten erwarten wir in der Regel keine Normalverteilung also die Abweichungen (Messfehler) sollen normalverteilt sein. Der KS-Test (Kolmogorow-Smirnow-Test) ist kein Nachweis der Normalverteilung sondern testet ob die Daten einer angenommenen Wahrscheinlichkeitsverteilung folgen.

Häufigkeiten

Häufigkeiten werden primär mit der Funktion APA2() und Tabelle() erstellt die Ausgabe der Konfidenzintervalle erfolgt die Funktion durch APA_CI .

Tabelle()

APA2()

smokers  <- c( 83, 90, 129, 70 )
patients <- c( 86, 93, 136, 82 )

n<-999
g =gl(3, n/3, labels = c("Control", "Treat A", "Treat B"))
g2<- g[sample.int(n)]
levels(g2)<- c("male", "female", "female")
data<- data.frame(g=g, g2=g2,
                  x=rnorm(n) )[sample.int(n)[1:78],]

 APA_CI( data, g )  
#  library(BayesianFirstAid)
#  x<-as.data.frame(table(data$g))$Freq
#  bayes.prop.test(x, rep(nrow(data), nlevels(data$g)), p=1/3)

Mittelwerte

Mittelwert ist ein Überbegriff von Maße der zentralen Tendenz. Folgende Mittelwerte können berechnet werden: Modus (Ausprägung mit höchster Häufigkeit), Median (sortierten Werte), Arithmetisches Mittel (was meist unter Mittelwert verstanden wird) Geometrisches Mittel (für logarithmische Verteilungen)

hkarz %>% Tabelle2(tzell[median], Lai, gruppe) 

set_my_options(median=list(style="IQR"))

APA2(tzell[1]+Lai~gruppe, hkarz, 
     type=c("auto", "median"),
     caption="Einfache Auswertung",
     test=TRUE, include.n = TRUE)

T-Test

Student's t-Test ist der ueberbegriff von Two Sample t-test, Welch Two Sample t-test und Paired t-test Dabei ist zu beruecksichtigen dass der Welch Test von ungleichen Varianzen ausgeht und bei mir die Standart -Funktion für den T-Test ist. Interessant ist auch dass Anova und T-Test identisch sind.

Zum T-Test gehört auch Cohen's d als Effektstärke.

# T-Test
 APA( t.test(m1 ~ geschl, varana, var.equal=TRUE))

 APA( t.test(m1 ~ geschl, varana))

 APA(aov( m1 ~ geschl, varana ))

APA2(m1 + m2 ~ geschl, varana, test = "t.test")

APA_Ttest(m1 + m2 ~ geschl, varana)

#varanax<-Melt2(m1+m2~nr,varana , key="time", value="m")
# broom::tidy(with(varana, t.test(m1,m2 ) ))
# broom::tidy(  t.test(m~time, varanax, var.equal=FALSE)) )

#APA_Ttest(m~time, varanax, paired = TRUE, include.mean=FALSE)
#APA_Ttest(m~time, varanax, paired = TRUE, type="U-Test", include.mean=FALSE)
#APA_Ttest(m~time, varanax, var.equal=TRUE, include.mean=FALSE)



#https://www.r-bloggers.com/is-it-time-to-ditch-the-comparison-of-means-t-test/

#exakt das selbe wie t.test(m1 ~ geschl, varana, var.equal=TRUE))
#t.test(m1 ~ geschl, varana, var.equal=TRUE))

res <-
summary(nlme::gls(m1 ~ geschl, varana,  weights = nlme::varIdent(form = ~ 1 |
geschl)))
res$tTable %>% fix_format() %>% Output(caption = "Generalized Least Squares")

Kreuztabellen

Kontingenztafeln, Kontingenztabellen oder Kreuztabellen sind Tabellen, die Kombinationen von Häufigkeiten darstellen. Hier gibt es wieder die Funktionen APA2() und Tabelle() sowie die R-Funktion xtabs(). Bei der Verwendung von xtabs kann die Ergebnisstabelle mit APA2 formatiert Tabellen erstellen. 2x2 Tabellen werden mit Häufigkeit und wahlweise mit Prozent (das Verhalten wird über margin = 2 gesteuert) ausgebeben. Berechnet werden mittels fisher.test() die Odds-Ratio mit 95%-CI (vcd:: oddsratio). Weiteres lässt sich mit der Option type = 2 eine Sensitivitätsanalyse erstellen. Bei NxM-Tabellen wird als Test-Statistik Pearson und der Kontingentkoeffizient berechnet, alternativ steht auch der Phi-Coefficient zur Verfügung (auch hier mit type = 2). Die Berechnung erfolgt hier mit der Funktion vcd ::assocstat.

APA_Xtabs()

2x2 Tabelle

Mit Tabelle Fischer Test und Diagnostiktest

hkarz$LAI<- factor(hkarz$lai, 0:1, c("pos", "neg"))
hkarz$Tzell<- cut(hkarz$tzell, 3, c("low", "med", "hig"))

xtab <- xtabs(~ LAI+ gruppe, hkarz)
APA2(xtab, caption="Harnblasenkarzinom")
APA_Xtabs(~ LAI+ gruppe, hkarz, include.percent=FALSE, include.test=TRUE)

NxMxO Tabelle

data(infert, package = "datasets")
infert$case  <- factor(infert$case ,1:0, c("case", "control") )

infert$spontaneous <- factor(infert$spontaneous)
infert$induced2    <- factor(infert$induced==0)

tab_1<- xtabs(~  case, infert)
tab_2x2<- xtabs(~ induced2 + case, infert)
tab_3x2<- xtabs(~ induced + case, infert)
tab_3x3<- xtabs(~ induced + education, infert)
tab_3x3x2<- xtabs(~ induced + education+case, infert)

#APA2(summary(tab_3x3x2))

(APA2(tab_1, include.test=TRUE, output=FALSE))
(APA2(tab_2x2, include.test=TRUE, output=FALSE))
(APA2(tab_3x2, include.test=TRUE, output=FALSE))
(APA2(tab_3x3, include.test=TRUE, output=FALSE))
(APA2(tab_3x3x2, include.test=TRUE, output=FALSE))

Sensitivity und Specificity

Sind spezielle Masszahlen bei 2x2 Tabellen Klassifikation (richtig/falsch) für xtabs und glm-Objekte können mit Klassifikation() berechnet werden, nicht zu vergessen ist die Kappa-Statistik die über MethComp() angefordert wird.

Sensitivität: richtig positive Rate eines Tests

Spezifität: richtig-negative Rate eines Tests

$$\begin{array} {rrr} A & B \ C & D \ \end{array}$$

$$Sensitivity = \frac{A}{A+C}$$ $$Specificity = \frac{D}{B+D}$$ Es gibt zwei Verfahren die anwendbar sind epiR epi.tests und epiR epi.tests Sensitivity, specificity and predictive value of a diagnostic test

require(epiR)
DF<- GetData("
Diagnose RT.qPCR Anzahl
krank positiv 111
krank negativ 12
gesund positiv 2
gesund negativ 62 ", 
             Tabel_Expand =TRUE, id.vars=1:2, output=FALSE)
DF$Diagnose<- factor( DF$Diagnose, rev(levels( DF$Diagnose)))
DF$RT.qPCR<- factor( DF$RT.qPCR, rev(levels( DF$RT.qPCR)))


APA2(epi.tests(xtabs(~ RT.qPCR + Diagnose, DF)))

APA_Xtabs(~ RT.qPCR + Diagnose, DF, include.diagnostic=TRUE, 
          caption="caret confusionMatrix")

fit1<- glm(Diagnose~RT.qPCR, DF, family = binomial)


  Klassifikation(fit1)$statistic[c(1,7,8),]

ROC Analyse

Dient der Auffindung des cut off value (Trennwert). Receiver Operating Characteristic (ROC), ist eine statistische Methode die ursprünglich in der Rundfunktechnik entwickelt wurde. Man versteht darunter eine Grafik die die Werte aus Sensitivität und Spezifität aufträgt, also Richtig-Positiv-Rate gegen die Falsch-Positiv-Rate.

Die optimale Klassifikations-Schwelle finden man, in dem man jenen ROC-Wert sucht, der den größten Normalabstand zur Diagonale des Diagramms aufweist. Weitere Info zu ROC findet sich unter https://www.hranalytics101.com/how-to-assess-model-accuracy-the-basics/

J = max(Sensitivity + Specificity - 1)

[@Sachs2006] Seite 159 Blutzucker - Beispieldaten (simuliert)

require(pROC)

blz.diab <- round(rnorm(100, mean = 115, sd = 20), 1)
blz.cntr <- round(rnorm(100, mean = 85, sd = 20), 1)
data <- data.frame(Gruppe = factor(c(
  rep("Diabetiker", length(blz.diab)),
  rep("Kontrollen", length(blz.cntr))
)),
Blutzucker = c(blz.diab, blz.cntr))


roc_curve <-  roc(data$Gruppe, data$Blutzucker)
plot(roc_curve, print.thres = "best", print.auc = TRUE)

Das selbe aber mit Regression daher sind die cut off -Werte nur indirekt interpretierbar

fit1  <- glm(Gruppe ~ Blutzucker, data, family = binomial)
x1 <- Klassifikation(fit1)
roc_curve   <- roc(x1$response, x1$predictor)

plot(roc_curve, print.thres = "best", print.auc = TRUE)
abline(v = 1, lty = 2)
abline(h = 1, lty = 2)

fit1 <- glm(gruppe ~ lai, hkarz, family = binomial)
fit2 <- glm(gruppe ~ lai + tzell, hkarz, family = binomial)
#thkarz <- as.data.frame(xtabs(~gruppe+lai, hkarz))
#fit2<- glm(Freq ~ gruppe*lai, thkarz, family = poisson())
x1 <- Klassifikation(fit1)
x2 <- Klassifikation(fit2)
#require(pROC)
roc1   <- roc(x1$response, x1$predictor)
roc2   <- roc(x2$response, x2$predictor)

plot(roc1, print.auc = TRUE, print.auc.y = 0.6)
plot(
  roc2,  lty = 2,  col = "blue",
  print.auc.y = 0.7,  print.auc = TRUE,
  add = TRUE
)
legend(
  "bottomright",
  legend = c("Lai", "Lai + T-Zell"),
  col = c(par("fg"), "blue"),
  lty = 1:2,  lwd = 2
)
roc.test(roc1, roc2)

Effektstärke

Mehr zum Thema Unstandardisierte Effektgrößen ist unter Regressionsanalyse zu finden.

Effektstärke ist eine dimensionslose Zahl zur Verdeutlichung der praktischen Relevanz von statistisch signifikanten Ergebnisse. Pearson-Korrelation r Cohens d

$$r={\frac {d}{{\sqrt {d^{2}+{\frac {(n_{1}+n_{2})^{2}}{n_{1}n_{2}}}}}}}$$

Cohens f2 $$f^{2}={\frac {R_{{included}}^{2}-R_{{excluded}}^{2}}{1-R_{{included}}^{2}}}$$

partielle Eta-Quadrat

$$\eta ^{2}={\frac {QS_{{{\rm {Effekt}}}}}{QS_{{{\rm {Effekt}}}}+QS_{{{{\rm {Res}}}}}}}$$

# APA2(tzell~gruppe,hkarz, type="cohen.d")  # APA_Effsize ist das gleiche

1+1

ANOVA/Regression: f, $\eta$, $\eta_part$ und d, $\Delta$, Hedge's g (bei Vergleich der Stichproben) sowie Odds Ratio, Risk Ratio, r, R

Korrelation

Der Korrelationskoeffizient ist ein dimensionsloses Maß für den Grad des linearen Zusammenhangs. Korrelationstabelle Interkorrelationen

 n<- 2*8
 e<- rnorm(n)*10
 DF<- data.frame(a=rnorm(n) + e,
                    b=rnorm(n), c=rnorm(n), d=rnorm(n) + e,
                    g=gl(2, 8, labels = c("Control", "Treat")))

Kreuztabelle: Cohen's w, $\phi$, Cramer's V, C

   APA_Correlation(~a+b+c, DF, caption="Formula a+b+c", output= "text")
#   in der Tabelle gibt es seltsame sonerzeichen ist nicht in caption

Likert-Skalen

Analyse von Likertskalen

 set.seed(1)
n<-100
lvs<-c("--","-","o","+","++")
DF2<- data.frame(
  Magazines=gl(length(lvs),1,n,lvs),
  Comic.books=gl(length(lvs),2,n,lvs),
  Fiction=gl(length(lvs),3,n,lvs),
  Newspapers=gl(length(lvs),5,n,lvs))
DF2$Comic.books[sample.int(n/2)]<- lvs[length(lvs)]
DF2$Newspapers[sample.int(n/2)]<- lvs[1]
DF2$Magazines[sample.int(n/2)]<- lvs[2]

DF2<- transform(DF2, Geschlecht= cut( rnorm(n), 2, Cs(m, f)))

Res1 <- Likert(~., DF2 )
APA2(Res1)

Gruppenvergleich 

APA2(Res2 <- Likert(.~ Geschlecht, DF2 ))
APA2(Res2, ReferenceZero=3, na.exclude=TRUE, type="freq")

Likert-Plots

#require(HH)
# ?likertplot
#
#windows(7,3)
#likertplot( Item   ~ .| Geschlecht , data=Res2$results,
#            main='',ylab="", sub="" ,
#            # col=brewer_pal_likert(5, "RdBu", "Geschlechtay80") ,
#            positive.order=TRUE,
#            as.percent = TRUE,
#            auto.key=list(space="right", columns=1)
#)
#data<- Res2$names
#data$mean<- Res2$m
#  barchart( Item~ mean |Geschlecht, Mymean2, xlim=c(1,5))
#windows(3,6)
#dotplot(Item ~ mean, data,
#        groups=Geschlecht, xlim=c(.085, 5.15),
#        type=c("p", "l"))

Rangreihen

Rangordnungen von Objekten können durch eine Transformation der Rangreihen in Intervallskalierte Merkmale überführt werden. Die Grundidee dieser Methode geht auf Thurstone (1927) nach dem "Law of Categorical Judgement" zurück. Dabei werden die kumulierten Häufigkeiten in Normalverteilte z-Werte übergeführt und aus diesen die Intervallskalierten Merkmalsauspraegungen gebildet [@Bortz2006] Seite 155.

library(PlackettLuce)

nlv <- 5
n <- 2 * 3 * nlv * 1
set.seed(n)

DF <-
  data.frame(
    Geschlecht = gl(2, n / 2, labels = c("Maennlich", "Weiblich")),
    Alter = gl(4, n / 4,   labels = c("20-29", "30-39", "40-49", "50-59")),
    Landwirtschaft = gl(2, n / 2, labels = c("konventionell", "biologisch"))
  )

Attribute <-
  as.data.frame(t(apply(matrix(NA, ncol = n, nrow = 5), 2,
                        function(x)
                          sample.int(5))))

Attribute[1, ] <- c(5, 1, 4, 2, 3)
Attribute[2, ] <- c(5, 1, 4, 2, 3)
Attribute[3, ] <- c(5, 2, 4, 3, 1)
Attribute[4, ] <- c(5, 1, 4, 3, 2)
Attribute[5, ] <- c(5, 1, 4, 3, 2)

Attribute[21, ] <- c(1, 2, 5, 4, 3)
Attribute[22, ] <- c(1, 4, 5, 3, 2)
Attribute[23, ] <- c(2, 5, 1, 4, 3)
Attribute[24, ] <- c(1, 4, 2, 5, 3)
Attribute[25, ] <- c(1, 4, 3, 5, 2)

attribute  <- c("Verfuegbarkeit",
                "Vielfalt",
                "Qualitaet",
                "Geschmack",
                "Preis")

Attribute<- dapply2(Attribute, function(x) factor(x, 1:5, attribute))

DF <- cbind(DF, Attribute)

 res <-
  Rangreihe( ~ V1+V2+V3+V4+V5,
             DF, 
             include.percent=FALSE, 
             order=FALSE, 
             include.na=FALSE,
             caption="Produkte aus konventioneller und biologischer  Landwirtschaft")

names(res)


R <- as.rankings(res$items)

mod <- PlackettLuce( R )
coef(mod)

summary(mod)


summary(mod, ref = "Verfuegbarkeit")
round(coef(mod, log = TRUE, ref = "Verfuegbarkeit") ,2)
round(coef(mod, log = TRUE ) ,2)


res$res$pc <-  round(coef(mod, log = FALSE) ,2)
res$res$log.pc <- round(coef(mod, log = TRUE) ,2)
res$res[order(res$res$pc,decreasing=TRUE),]

Tests

Quelle: https://cran.r-project.org/doc/contrib/Ricci-refcard-regression.pdf

Methodenvergleich

Oft interessiert die Zuverlässigkeit und Reproduzierbarkeit einer Diagnose. Die Beurteilung kann dabei durch einen Bewerter (Messverfahren) in wiederholter Form erfolgen und wird dann als Intra-Rater bezeichnet oder die Beurteilung eines Merkmals erfolgt durch mehrere Bewerter (Messverfahren). und hier spricht man von Inter-Rater. Die Methode der Beurteilung der Übereinstimmung hängt von den jeweiligen Verteilungseigenschaften ab.

Bei Nominalen wird abgezählt und die Rate der Übereinstimmung bewertet (Cohen-Koeffizient) Bei Ordinalen-Daten werden die gewichteten Übereinstimmungen ausgezählt (gewichteter Cohen-Koeffizient). Bei metrischen(stetigen) Daten werden die Differenzen beurteilt (Bland-Altman-Methode oder auch Tukey Mean Difference).

Bland-Altman-Methode Bias (d) systematische Abweichung Messfehler (s) Standardabweichung der Differenz Limits of agreement (LOA) Intervall von 95 (entspricht d+-1.96 -> es wird eine Normalverteilung unterstellt).

MetComp

Die generische Funktion MetComp() kann sowohl Kappa als auch Tukey-means berechnen. Kappa kann aber auch über die xtab() und APA2 berechnet werden. Wobei hier nur 2x2-Tabellen untersucht werden können, hingegen sind bei Kappa() auch mehrere ordinale Kategorien erlaubt.

Ähnliche Methode ist ICC die aber diese zählt eher zur Reliabilitätsanalyse.

Funktionen: MetComp, BAP, und Kappa

Beispiel Altman and Bland [@Giavarina2015]

set.seed(0815)
Giavarina <- data.frame(
A=c(1,5,10,20,50,
    40,50,60,70,80,
    90,100,150,200,250,
    300,350,400,450,500,
    550,600,650,700,750,
    800,850,900,950,1000),
B=c(8,16,30,14,39,
    54,40,68,72,62,
    122,80,181,259,275,
    380,320,434,479,587,
    626,648,738,766,793,
    851,871,957,1001,980),
group= sample(gl(2, 15, labels = c("Control", "Treat")))
)

MetComp(~A+B, Giavarina, caption = "Giavarina")

Beispiel Botulinum A

Sachs Angewandte Statistik Seite 627

Botulinum <- data.frame(
  A= factor(c(rep(1, 14), rep(1, 3),
              rep(0, 5),rep(0, 18)),
            1:0, c("+", "-")),
  B= factor(c(rep(1, 14), rep(0, 3),
              rep(1, 5),rep(0, 18)),
            1:0, c("+", "-")))

# APA2(xtabs(~A+B, Botulinum), test=TRUE)

MetComp(~A+B, Botulinum)

 xt <-xtabs(~A+B, Botulinum)
 Klassifikation(xt)$statistic[c(1,3,4), ]

Giavarina <- transform(Giavarina, C = round( A + rnorm(30,0,20)),
                D = round( A + rnorm(30,0,10) + A/10 ),
                E = round( A + rnorm(30,5,10) + (100-A/10) ))


#ICC2(~A+E, Giavarina, caption="ICC (Korrelationen)")

# A - Goldstandart

x <- MetComp_BAP(~A+B+E, Giavarina)
# x %>% Output("BA-Analyse der Messwertreihe")
plot(x)

x <- MetComp_BAP(~A+E+B, Giavarina)
# x %>% Output("BA-Analyse der Messwertreihe")
plot(x)

Verschiedene Situationen

set.seed(0815)

n<-100
DF<- data.frame(
  A=rnorm(n, 100,50),
  B=rnorm(n, 100,50),
  C=NA,  D=NA,  E=NA,
  group= sample(gl(2, n/2, labels = c("Control", "Treat")))
)

cutA<-mean(DF$A)
DF <- transform(DF, C = round( A + rnorm(n, -5, 20)),
                D = round( A + rnorm(n,0,10) + A/10 ),
                E = A + ifelse(A<cutA, A/5, -A/5 )+ rnorm(n, 0, 10)
)

x<- MetComp_BAP(~A+C, DF)
plot(x)

x<- MetComp_BAP(~A+B, DF)
plot(x)

x<- MetComp_BAP(~A+D, DF)
plot(x)

x<- MetComp_BAP(~A+E, DF)
plot(x)

Regressionsanalyse

Regressionsanalysen sind eine Gruppe von statistische Verfahren, die Beziehungen zwischen Variablen modellieren.

Regression versus Kausalität: Eine Regression liefert Aussagen über statistische Zhä jdh ih üb di Kliä Zusammenhänge, jedoch nicht über die Kausalität der Variablen (ob eine kausale Beziehung besteht, kann nicht die Regeion beantoten

Darunter fallen alle Methoden wie ANOVA, lineare-Regression, logistische-Regression usw. ein anderer Überbegriff ist Generalized Linear Model (GLM).

T-Test und ANOVA unterscheiden sich nicht wenn $d = 1$ ist. In diesen Fall sind F-Wert und T-Wert dieselben Zahlen $F=T^2$. Die Regressionsanalyse ist die Verallgemeinerung der (Mittelwert)-Analyse und liefert exakt die gleichen Ergebnisse wie die ANOVA. Der Unterschied ist, dass die Regression mehr Informationen liefert. Aber die die statistische Hypothese ist verschieden.

ANOVA: $H_0 : \mu_1 = \dots = \mu_n$

Regression $H_0 : \beta_1 = \beta_2 = 0$ wobei $y_i=\beta_0+\beta_1x_i + \dots$

Bei SPSS wird die ANOVA mit vielen verschiedenen Bezeichnungen ausgegeben: als Test-der-Zwischen-Subjekt-Effekte, Test-der-Inner-Subjekt-Effekte, Multivariate-Test, usw.. Aber hinter den Begriffen steht immer eine ANOVA oder noch allgemeiner ein GLM (Generalisiertes-Lineares-Modell). Ich verwende immer nur die Begriffe ANOVA und Regressionsanalyse.

Effektstärke, Eta-Quadrat Das Eta-Quadrat Eta² ist nur eine Variante von Maßzahlen zur Effektstärke und nicht die Effektstärke und nicht immer ist das Eta² die passende Zahl. Eta² hat zudem den prinzipiellen Nachteil, dass es "unanschaulich" ist. Trotzdem verwende ich meist, aus pragmatischen Gründen, auch wenn es falsch ist das Eta-Quadrat. Anmerkung zur Effektstärke das APA-Manual schreibt:

include some measure of effect size ..., whenever possible, provide a confidence interval of each effect size. Effect size may be expressed in the original units (mean, ... regression slope). It can be valuable to report Effect size, also in some standardized or unites-free unit, … however, it, to provide the reader enough information to assess the magnitude of the observed effect(APA [1] Seite 34).

Frei interpretiert heißt das; unstandardisierte Effektgrößen (Regressionskoeffizienten mit 95%-CI) müssen angeben werden, andere können angegeben werden. Die "magnitude of the observed effect" ist am besten über die Effectplots darzustellen (Fox [6]). Es gibt zwei Arten von (standardisierten) Effektstärken die d-Familie, die Unterschiede zwischen Gruppen betrachtet, und die r-Familie, welche ein Maß für Zusammenhänge zwischen Daten ist [@Hemmerich]. Unstandardisierte Effektgrößen sind Differenzen von Gruppenmittelwerten (raw mean difference) und unstandardisierte Regressionskoeffizienten. Der p-Wert ist keine Effektstärke.

ANOVA und multiple Vergleiche

Die ANOVA kann als Erweiterung des t-Tests betrachtet werden und ist fast das gleiche wie die lineare Regression.

Quelle: http://www.uni-koeln.de/~a0032/statistik/texte/mult-comp.pdf

Kontraste

R bietet die o.a. Standard-Kontraste, die über die folgenden Funktionen erreichbar sind: contr.treatment(k,base=j) (j=Nummer der Vergleichsgruppe)contr.sum(k), contr.helmert(k), contr.poly(k)

$\alpha$-Adjustierungen

R bietet $\alpha$-Adjustierungen sowohl unabhängig von irgendwelchen Vergleichsverfahren an: p.adjust die Standard-Adjustierungen (Bonferroni, Holm, Hochberg und Bejamini), multcomp die Verfahren von Westfall und Shaffer,

#' breaks ~ wool + tension 
#' warpbreaks %>% Tabelle2(breaks, by= ~ wool + tension)
fm1 <- aov(breaks ~ wool + tension, data = warpbreaks)
APA2(fm1, caption="ANOVA")

SPSS vs R

Vergleich der Ergebnisse mit SPSS dort heit die ANOVA auch Tests of Between-Subjects Effects und wird über den Button Univariat aufgerufen.

Source|Type III Sum of Squares|df|Mean Square|F|Sig.|Partial Eta Squared ------|-----------------------|--|-----------|-|----|------------------- Corrected Model|12587|3|4196|124|0.000|0.43 Intercept|3827|1|3827|113|0.000|0.19 grade|3329|1|3329|98|0.000|0.17 treatment|241|1|241|7|0.008|0.01 stdTest|6971|1|6971|206|0.000|0.30 Error|16468|486|34||| Total|65279|490|||| Corrected Total|29056|489||||

#'  R Squared = .43 (Adjusted R Squared = .43)
#' Levene's Test of Equality of Error Variances F=1.8, df1=3,df2=486, p=.146
#' 
schools2 <-  transform(
  schools,
  classroom = factor(classroom),
  grade = factor(grade),
  treatment = factor(treatment),
  score10 = score > 10,
  score2 =  round((log((
  score + 20
  )) - 1.78) / .023)

  )

  fit <- aov(score ~ grade + treatment + stdTest, schools2)

  APA2(fit, include.eta = TRUE)

#APA_Validation(fit, include.levene = TRUE,
#               include.aic = FALSE, include.ftest = FALSE,include.deviance = FALSE,
#               include.heteroskedasticity = FALSE,
#               include.durbin = FALSE, include.normality = FALSE)

Tukey HSD (honestly significant differnce)

Das wohl bekannteste Verfahren ist das von Tukey (für gleiche ni) bzw. von Tukey & Kramer(für ungleiche ni). Geprüft wird wie bei der Varianzanalyse: alle k Mittelwerte sind gleich. 

TukeyHSD(fm1, "tension", ordered = TRUE) %>%
  APA2(caption="TukeyHSD" )

#' plot(TukeyHSD(fm1, "tension"))
#' levels(warpbreaks$tension)
#' Lm Split

fm1_split <-  summary(fm1,
                      split=list(tension=list( M=1,  H=3, L=2)),
                      expand.split=FALSE)
APA2(fm1_split)

das gleiche mit multcomp

require(multcomp)

fit_Tukey <- glht(fm1,
linfct = mcp(tension = "Tukey"),
alternative = "less")

APA2(fit_Tukey, caption = "APA2: multcomp mcp Tukey")

Kontraste

### contrasts for `tension'
K <- rbind(
  "L - M" = c(1,-1,  0),
  "M - L" = c(-1,  1,  0),
  "L - H" = c(1,  0,-1),
  "M - H" = c(0,  1,-1)
  )

  warpbreaks.mc <- glht(fm1,
  linfct = mcp(tension = K),
  alternative = "less")
  APA2(warpbreaks.mc, caption = "APA2: multcomp mcp mit Contrasten")
### correlation of first two tests is -1
#cov2cor(vcov(fm1))

### use smallest of the two one-sided
### p-value as two-sided p-value -> 0.0232
#summary(fm1)

Interaction

fm2 <- aov(breaks ~ wool * tension, data = warpbreaks)
APA2(fm2)
x <- TukeyHSD(fm2, "tension",
              ordered = TRUE)
APA2(x, caption="Interaction: TukeyHSD" )

warpbreaks$WW<-interaction(warpbreaks$wool,warpbreaks$tension )
mod2<-aov(breaks~WW, warpbreaks)
APA2(mod2, caption="ANOVA interaction haendich zu den Daten hinzugefuehgt")

Lineare Regressionsanalyse (lm)

Die Bezeichnung Regression stammte historisch gesehen von Francis Galton, er untersuchte den Zusammenhang der Körpergröße von Eltern und Kindern (Regression to the Mean). Ziel der Regressionsanalyse ist eine funktionale Beziehung zwischen zwei Größen zu finden.[1] Mathematisch lässt sich das folgend formulieren ($\hat{Y} = c + bX + e$), dabei ist $X$ die unabhängige und $\hat{Y}$ die abhängige Variable und $e$ der statistische Fehler. Gesucht wird, die “Formel” der Gerade, die in der graphischen Darstellung durch den Mittelwert verläuft. Die Regression ist quasi die Erweiterung der Korrelationsanalyse die ja die Stärke des Zusammenhangs ermittelt.

Linear model

fit1 <- lm(tzell ~ gruppe, hkarz)
fit2 <- lm(tzell ~ gruppe + Lai, hkarz)
fit3 <- lm(tzell ~ gruppe * Lai, hkarz)

APA_Table(fit1, fit2, fit3, type = "long", caption = "Regression")

Kennzahlen der (linearen) Regressionsanalyse)

~~Die library(psycho) kann ueber analyze() Betas recht einfach ausgeben. Das geht sowohl bei lm als auch bei lmer. Aber die Ergebnisse sind anders als bei PSPP und bei meiner Funktion?? - Gehoert also noch ueberprueft.~~

## Geht nicht mehr easystats
# 
# Note: Many functions of the 'psycho' package have been (improved and) moved to other packages of the new 'easystats' collection (https://github.com/easystats). If you don't find where a function is gone, please open an issue at: https://github.com/easystats/easystats/issues
# Error in analyze(fit) : could not find function "analyze"

 #library(psycho)
 # df <- psycho::affective  # Load a dataset from the psycho package
 #  #df <- standardize(df)  # Standardize all numeric variables
 # fit <- lm(Age ~ Salary, data=df)  # Fit a Bayesian linear model
 # results <- analyze(fit)  # Format the output
 # APA2(results)

#  library(lmerTest)
#  fit <- lmerTest::lmer(Sepal.Length ~ Sepal.Width + (1|Species), data=iris)
#  results <- analyze(fit)
#  APA2(results)

Logistische Regression Logit-Modell

Das Logit-Modell ist eine Gruppe von Regressionsverfahren die als Klassifikationsverfahren verwendet werden. Zähldaten.

Binomiale Regressionsanalyse

Die binomiale logistische Regression wird angewendet, wenn geprueft werden soll, ob ein Zusammenhang zwischen einer binaeren Ziel-Variable und einer oder mehreren Einfluss-Variablen besteht. Die folgenden Beispiele stammen von Backhaus Seite 456.

Anmerkung zu SPSS - warum sind die Ergebnisse von R und SPSS so verscheieden?

Bei SPSS gibt es zwei verschiedene Methoden klassisch und ueber GLM Klassisch kann SPSS aber keine Faktoren ins Model aufnehmen und beherrscht nur den Wald-Test. Bei GLM rechnet SPSS die Intercept anderst kann aber Faktoren berechnen und es kounnen verschiedene Tests gewaehlt werden. Bei SPSS sind die Ergebnisse der Parameter identisch aber das Intercept ist unterschiedlich. Bei SPSS wird quasi die ANOVA bei den Koeffizienten ausgegeben und die Inercepts werden je nach Einstellung unterschiedlich berechnet - was problematisch wird wenn man aus den Parametern die zugehourigen Wahrscheinlichkeiten errechnen will. Bei R wird der Parameter b mit dem LR-Test getestet und Faktoren werden automatisch erkannt und das Intercept wird immer gleich berechnet. Zusammenfassend: R und SPSS werden zwei verschiedene Hypothesen um die p-Werte zu bestimmen. Und bei SPSS werden die Intercepts quasi kreativ berechnet.

Ergebnisse von SPSS ueber Regression

Parameter|B|S.E.|Wald|df|Sig.|Exp(B) ---------|-|----|-----|--|---|----- tzell|0.2|0.094|4.57|1|0.032|1.222 lai|2.21|0.877|6.3|1|0.012|9.075 Constant|-14.64|6.329|5.37|1|0.021|0

Ergebnisse von SPSS ueber GML

Parameter|B|Std. Error|Lower|Upper|Wald Chi-Square|df|Sig. ---------|-|----------|-----|-----|---------------|--|---- Threshold [gruppe=1.00]|12.44|6.52|-0.35|25.2|3.63|1|0.057 tzell|0.2|0.09|0.01|0.34|4.5|1|0.032 [lai=.00]|-2.2|0.88|-3.9|-0.48|6.3|1|0.012

Der Wald-Test bei R kann mit car::Anova(fit, type = “III”, test.statistic = “Wald”) aufgerufen werden. Voraussetzung ist dass das Modell streng linear ist.

Beispiel Harnblasenkarzinom

fit1 <- glm(gruppe~tzell, hkarz, family = binomial)
APA_Table(fit1, type="long")

Signifikanz der Regressionskoeffizienten; bei SPSS wird der Wald-Test verwendet, ich verwende den LRT (ein Vergleich findet man unter http://thestatsgeek.com/2014/02/08/wald-vs-likelihood-ratio-test/)

Odds ratio* OR Exp(b) - Exponent des Regressionskoeffizienten Odds Ratios ist nur bei der logistischen Regression sinnvoll. Logarithmierte Odds (Logits, Effekt-Koeffizienten)

res <- broom::tidy(fit1) 
cbind(res[1:2], 
      confint(fit1),
      res[5]) %>% 
  fix_format %>% Output("broom::tidy(fit1): Odds Ratios")

Signifikanz des Regressionsmodells: Gütekriterien auf Basis der LogLikelihood-Funktion Devianz -2LL-Wert H0: Modell besitzt eine perfekte Anpassung Gute Werte sind -2LL nahe 0 und p-Wert nahe 1 deviance(fit1)

LR-Test (Likeihood Ratio-Test) H0: Alle Regressionskoeffizienten sind gleich Null. Chi-Quadret möglichst hoch p<0.05

Log-likelihood values cannot be used alone as an index of fit because they are a function of sample size but can be used to compare the fit of different coefficients. Because you want to maximize the log-likelihood, the higher value is better. For example, a log-likelihood value of -3 is better than -7

Der Wald Test sind eigentlich zwei Methoden die erste ist der Test ob irgendein Koeffizient signifikant ist, die zweite Methode welcher Koeffizient signifikant ist. Der erstere Wald Test (Model Coefficients) ist in der Regel wenig informativ. Der zweite Wald-Test ist die ANOVA (Type II) andere Bezeichnung ist Wald Chi-square Test oder Analysis of Deviance Table (Die Anova testet die Hypothese ob der zugehörige Koeffizient Null ist.)

Likelihood ratio test = Wald-Test Analysis of Deviance Table = Wald-Test (Type II), Anova Log Likelihood die Zahl kann für sich nicht interpretiert werden AIC, BIC (Informationskriterium) kleiner Wert steht für eine höhere Informativität der Models. Akaikes Informationskriterium (AIC), Bayesian Information Criterion (BIC)

#--Omnibus Test  Fuer F-Test :  lmtest::waldtest(fit1)
 lmtest::lrtest(fit1) %>% fix_format() %>% Output("lmtest::lrtest(fit1): LR-Test")

Pseudo-R-Quadrat 

Goodness <- function(x) {
  cbind(round(broom::glance(x)[, c(6,  3, 4, 5)], 1),
  round(R2(x), 2),
  round(RMSE(x)[2], 2))
}
fit1 %>% Goodness %>% Output()

Weitere Betrachtung ist über Klassifikationstabelle und dazu gibt es die Tests Hosmer Lemeschow-Test und Press Q-Test.

Klassifikation(fit1)

fit1 <- glm(gruppe~tzell+factor(lai), hkarz, family = binomial)
x<- APA2(fit1, include.odds=TRUE )
Nagelkerke<- R2(fit1)[3]

# Interpretation

txt_log_reg <-  paste("Eine logistische Regressionsanalyse zeigt, dass sowohl das Modell als Ganzes (",
 APA(fit1), 
 ") als auch die einzelnen Koeffizienten der Variablen signifikant sind.",
 "Steigen die T-Zelltypisierung  um jeweils eine Einheit, 
 so nimmt die relative Wahrscheinlichkeit eines Krank/Gesund um OR =",  x$odds[2], 
 "zu. Ist die  T-Zelltypisierung positiv so nimmt die  relative Wahrscheinlichkeit um OR= ", x$odds[2],
 "Das R-Quadrat nach Nagelkerke beträgt",round(Nagelkerke,2), 
" was nach Cohen (1992) einem starken Effekt entspricht."  )
# 

#  Quelle Text: http://www.methodenberatung.uzh.ch/de/datenanalyse/zusammenhaenge/lreg.html

Output: r txt_log_reg

# Wahrscheinlichkeiten T-Zell
fit1 <- glm(gruppe ~ tzell , hkarz, family = binomial)
t.zell<-    c(50,55,60,65,70,75,80) #seq(50,80, by=1)  

i<- coef(fit1)["(Intercept)"]
b<-coef(fit1)["tzell"]

z <- i + b*t.zell
p<- 1/(1+exp(z)) 

cbind(t.zell, p=round(p,2))

Poisson Regression

Poisson Regression für Zaehldaten und Kreuztabellen.

#-- SPSS kodiert die Gruppe 3 als Referenz
poisson_sim$prog <-
  factor(poisson_sim$prog, c("vocation", "general",  "academic"))
  fit5 <- glm(num_awards ~ prog + math,
  poisson_sim, family =  poisson())

  poisson_sim %>% Tabelle2(num_awards, prog, math)

  APA_Table(fit5)

The output above indicates that the incident rate for [prog=academic] is 2.042 times the incident rate for the reference group, [prog=vocation]. Likewise, the incident rate for [prog=general] is 0.691 times the incident rate for the reference group holding the other variables at constant. The percent change in the incident rate of num_awards is an increase of 7% for every unit increase in math. http://www.ats.ucla.edu/stat/spss/dae/poissonreg.htm

Kreuztabellen

Intresantes Detail Kreuztabellen, logistische-Regression und poisonn-Regression liefern identische Ergebnisse siehe Beispiel unten.

APA2(xtabs(~ gruppe + lai, hkarz), test = TRUE, type = "fischer")
fit1 <- glm(gruppe ~ lai, hkarz, family = binomial)
thkarz <- as.data.frame(xtabs(~ gruppe + lai, hkarz))
fit2 <- glm(Freq ~ gruppe * lai, thkarz, family = poisson())

APA_Table(fit1, include.odds = TRUE)
APA_Table(
fit1,
fit2,
include.odds = TRUE,
include.b = FALSE,
include.se = FALSE,
type = "long"
)

Hierarchische Modelle

Mixed Designs, Blockdesign, Hierarchische Modelle vs. feste/zufällige Effekte, Mehrebenenanalyse, Messwiederholung die Begriffe bedeuten alle das gleich es gibt keinen Unterschied in der Berechnung.

lmer_fit1<-lmerTest::lmer(score ~ agegrp+trial + (1|id), MMvideo)

lmer_fit2<-lmerTest::lmer(score ~ agegrp*trial + (1|id), MMvideo)
APA_Table(lmer_fit1, lmer_fit2, type="long")


# MySet()
# library(gridExtra)
#   windows(8,4)
#   p1 <- plot(effect("trial",lmer_fit1), multiline=TRUE, main="")
#   p2 <- plot(effect("agegrp*trial",lmer_fit2), multiline=TRUE, main="")
# 
#   grid.arrange(p1,p2,ncol=2)
# 
#  library(coefplot)
#  windows(4,3)
#   coefplot(lmer_fit2, intercept=F, xlab="b (SE)")
# 
#    windows(4,3)
#   multiplot(lmer_fit1, lmer_fit2, intercept=F, xlab="b (SE)")

MANOVA

MANOVA ist eine Methode der Varianzanalyse und prüft multivariate Mittelwertunterschiede. Diskriminanzanalyse ist ein multivariates Verfahren zur Analyse von Gruppenunterschieden.

Vorgehensweise

  1. Deskriptive Analyse mit Box-Plots

  2. Prüfung der Mittelwertunterschiede (MANOVA)

  3. Diskriminanzanalyse

    • Formulierung und Schätzung der Diskriminanzfunktion
    • Prüfung der Diskriminanzfunktion
    • Prüfung der Merkmalsvarianblen
    • Klassifikation neuer Elemente

Beispiel Manova von ucla.edu https://stats.idre.ucla.edu/spss/output/one-waymanova/

m_data<-GetData("C:/Users/wpete/Dropbox/3_Forschung/R-Project/stp25data/extdata/manova.sav")
m_data$GROUP<- factor(m_data$GROUP, 1:3, c("website", "nurse ", "video tape" ))

z<- as.matrix(m_data[,-1])

fit1<- manova(z ~ m_data$GROUP)
APA2(fit1)
#summary(fit1)$Eigenvalues
#library(MASS)
#fit2 <- MASS::lda(GROUP ~ ., data=m_data)
#APA2(fit2)
#plot(fit2)

Voraussetzungen

Im Folgenden die wichtige Voraussetzungen der Regressionsanalyse

Assumptions of the regression model (nach Gellman Hill Seite 45)

We list the assumptions of the regression model in decreasing order of importance. 1. Validity.
2. Additivity and linearity. The most important mathematical assumption 3. Independence of errors 4. Equal variance of errors 5. Normality of errors.

General principles

Our general principles for building regression models for prediction are as follows:

  1. Include all input variables that, for substantive reasons, might be expected to be important in predicting the outcome.
  2. It is not always necessary to include these inputs as separate predictors—for example, sometimes several inputs can be averaged or summed to create a “total score” that can be used as a single predictor in the model.
  3. For inputs that have large effects, consider including their interactions as well.
  4. We suggest the following strategy for decisions regarding whether to exclude a variable from a prediction model based on expected sign and statistical significance (typically measured at the 5% level; that is, a coefficient is “statistically significant” if its estimate is more than 2 standard errors from zero):
  5. If a predictor is not statistically significant and has the expected sign, it is generally fine to keep it in. It may not help predictions dramatically but is also probably not hurting them.
  6. If a predictor is not statistically significant and does not have the expected sign (for example, incumbency having a negative effect on vote share), consider removing it from the model (that is, setting its coefficient to zero).
  7. If a predictor is statistically significant and does not have the expected sign, then think hard if it makes sense. (For example, perhaps this is a country such as India in which incumbents are generally unpopular; see Linden, 2006.) Try to gather data on potential lurking variables and include them in the analysis.
  8. If a predictor is statistically significant and has the expected sign, then by all means keep it in the model.

Multikollinearität, Autokorrelation, Heteroskedastizität Diese speziellen Eigenschaften sind laut Gelman nicht besonders wichtig. Die Prüfung erfolgt über die grafische Darstellung der Residuen. Die Tests für Multikollinearität (VIF) und für Autokorrelation (Durbin-Watson-Test) sind weniger aussagekräftig. Ich verwende statt der sperrigen Begriffe Heteroskedastizität usw. den Begriff "Gleichheit der Varianzen".

Von Multikollinearität, bzw. Kollinearität, spricht man, wenn zwei oder mehrere erklärende x Variablen hoch untereinander korreliert sind. Anzeichen ist Hohes R2 und wenige signifikante Koeffizienten. Test über Variance Inflation Factor VIF VIF>10 Multikollinearitätsproblem

variance inflation factor, aka VIF values over 5 are troubling. should probably investigate anything over 2.5. Quelle: https://hlplab.wordpress.com/2011/02/24/diagnosing-collinearity-in-lme4/

 # Kopie der \link{car::vif} funktion
  library(car)
  fit<-lm(prestige ~ income + education, data=Duncan)
  car::vif(fit)
  VIF2(fit)

Variablenauswahl

Ein Kriterium für die Aufnahme einer Variablen ist das Schrittweise Verfahren (Step). Ein anderes ist der sachlogische Bezug also alle Variablen die sachlogisch Begründet werden gehören ins Model. Ich gehe immer folgenderweise vor: Ich verwende immer beide Varianten (forward und backward). Anschließend überprüfe ich die Ergebnisse Visuell auf einflussreiche Fälle. Das Heisst ich überprüfe ob die auf einzelne Ausreiser beruht. Im nächsten Schritt passe ich das Modell an, es werden also eventuell Variablen entfernt oder hinzugefügt. Schrittweises Verfahren (Choose a model by AIC in a Stepwise Algorithm)

forward selection: die Einflussgröße mit dem kleinsten p-Wert wird als Kandidat ins Modell aufgenommen und der Likelihood- Wert bestimmt. backward elimination: prüft ob auf eine Variable verzichtet werden kann Ich verwende immer beide Varianten (forward und backward)

Regressionsdiagnostik

Überprüft die Linearität, Independence of error, Equal variance of errors usw. Die Regressionsdiagnostik untersucht die Eigenschaften der Residuen. Die Residuen sind dabei die Abweichungen vom Model. Die Untersuchung erfolgt dabei über Grafik zu den Residuen.

loglik Likelihood Ratio Test The log-likelihood from the intercept-only restricted model -2LL: The LL (log-likelihood from the fitted model) llhNull (The log-likelihood from the intercept-only restricted model), G2 (Minus two times the difference in the log-likelihoods)

Homogeneity of Variances Levene Computes Levene's test for homogeneity of variance across groups. Bartlett Test of Homogeneity of Variances

Residual RMSE values should be low (<0.5 and <0.3, respectively). SigmaResidual standard error When the residual standard error is exactly 0 then the model fits the data perfectly (likely due to overfitting) R-Quadrats Cox und Snell R2: [ 0.2 = akzeptabel, 0.4 = gut ] Nagelkerke R2: [ 0.2 = akzeptabel, 0.4 = gut, 0.5 = sehr gut] McFaddens R2: [ 0.2 = akzeptabel, 0.4 = gut ] (see pR2)

Heteroskedastizität

Heteroskedastizität beschreibt einen Sachverhalt, in dem die Varianz der Residuen variiert. Konsequenzen: Trotz Heteroskedastizität sind die Parameter - Schätzer weiterhin unverzerrt, jedoch ineffizient (keine Minimum Varianz mehr), also die Hypothesentests werden unzuverlässig. Abhielfe: robuste Standardfehler.

Heteroskedasticity Breusch-Pagan test against heteroskedasticity.

Autokorrelation

Durbin-Watson Test for Autocorrelated Errors. car::durbinWatsonTest (Autokorrelation erster Ordnung) Autocorrelation Durbin-Watson test for autocorrelation of disturbances. Ist nur bei Zeitreihendaten sinnvoll.

Erkennen von Autokorrelation, der DW d Test weist folgendes Intervall auf:

0<d<4

0: Extrem positive Autokorrelation. 4: Extrem negative Autokorrelation 2: Keine Autokorrelation

Breusch-Godfrey Test der BG Test ist flexibler, da er auf Autokorrelation höherer Ordnung testet und auf autoregressive Modelle angewendet werden kann. (In der Literatur wird dieser Test auch LM-Test auf Autokorrelation genannt.)

Abhielfe: Da die Schätzer nicht verzerrt sind, besteht kein grundsätzlicher Handlungsbedarf.

Multikollinearität

Multikollinearität beschreibt einen Sachverhalt, in dem eine erklärende Variable eine Linearkombination einer anderen oder mehrerer anderer erklärender Variablen ist.

Erkennen von Multikollinearität
Hohes R2 und geringe t-Werte.
Der einfache Korrelationskoeffizient zwischen zwei erklärenden Variablen ist sehr hoch.

Abhilfemaßnahmen Entweder nichts unternehmen oder Überflüssig(st)e Variablen entfernen. VIF variance inflation factor. VIF values over 5 are troubling, should probably investigate anything over 2.5.

Diagnose_Plot <- function(data,
                          x,
                          y,
                          fit,
                          main = "",
                          ylab = "",
                          xlab = "",
                          level = 0.95) {
  q <- data[, x]
  predicted.intervals <-
    predict(fit,
            data[x],
            interval = 'confidence',
            level = level)

  plot(data[, x],
       data[, y],
       # col = 'deepskyblue4',
       xlab = xlab,
       ylab = ylab,
       main = '(a)')

  title("Observed data", line = .5)

  lines(q, predicted.intervals[, 1], col = 'deepskyblue4', lwd = 2)
  lines(q, predicted.intervals[, 2], col = 'black', lwd = 1)
  lines(q, predicted.intervals[, 3], col = 'black', lwd = 1)


  plot(fit, 1, main = "(b)",  cex.caption = .8)
  plot(fit, 2, main = "(c)",  cex.caption = .8)
  plot(fit, 3, main = "(d)",  cex.caption = .8)
  plot(fit, 5, main = "(e)",  cex.caption = .8)
  influence_year <- car::influencePlot(fit, main = '(f)')
  title("influence Plot", line = .5)

  invisible(influence_year)

}

#fit <- lm(Age ~ Concealing, data=df)
#par(mfrow = c(2, 3))
#infFall<- Diagnose_Plot(
#  df, "Concealing","Age",
#  fit,
#  xlab="",  ylab = 'Age' 
#)

Der Diagnose-Plot zeigt Auffälligkeiten
hinsichtlich der Residuen daraus kann abgeleitet werden das die Verwendung des linearen Models nicht vertretbar ist.

Strukturprüfende Verfahren

Multivariate Methoden und Skalen-Analyse

Validität, Objektivität kann nicht mit Statistik bewertet werden. Reliabilität, Interne Konsistenz (Cronbach’s Alpha), Trennschärfe (bei SPSS Korrigierte-Item-Skala-Korrelation). Validität wird in interne (ist durch die Randomisierung gegeben) und externe (Übertragbarkeit) unterschieden. Die Reliabilität kann berechnet werden wenn genug Daten n>60 und mehr als 3 Items vorliegen. Ein alpha>0.08 ist gut.

Cronbach’s Alpha

Verarbeitung <- Reliability(~ F5+F16+F22+F9+F26+F6+F35+F33+F12+F34+F4, fkv, check.keys =TRUE)
Coping <- Reliability(~ F7+F8+F17+F14+F15+F18+F19+F1+F13+F20, fkv, check.keys =TRUE)
Vertrauen <- Reliability(~ F28+F27+F31+F29, fkv, check.keys =TRUE)
Religion <- Reliability(~F21+F25+F30+F23+F24, fkv, check.keys =TRUE)
Distanz <- Reliability(~F3+F2+F10+F11, fkv, check.keys =TRUE)
#Verarbeitung %>% APA2()

Alpha(Verarbeitung, Coping, Vertrauen, Religion, Distanz) %>% Output()

ICC

ICC Intra-Klassen-Korrelation. Kopie von ICC aus dem Packet psych.

 sf <- GetData("
               J1 J2 J3 J4 J5 J6
               1  1  6  2  3  6
               2  2  7  4  1  2
               3  3  8  6  5 10
               4  4  9  8  2  4
               5  5 10 10  6 12
               6  6 11 12  4  8")

# #Quelle  http://www.personality-project.org/r/book/Chapter7.pdf
ICC(sf)

Mediation, Moderation, SEM, Sobel-Test

Der Sobel Test ist der indirekt Effekt also das Produkt aus Mediator und Einfluss-variable. Das Verfahren ist alt (1986) moderneres Verfahren ist bootstrapping oder auch Struktur-Gleichung-Modelle (SEM).

Anmerkungen zu Methoden

Hypothesen vs. Hypothesentest

Hypothesen und Hypothesentest sind nicht das Gleiche. Wenn ein Signifikanztest durchführt wird, ist das nicht automatisch ein Hypothesentest sondern das ist meistens ein explorativer Signifikanztest. Hypothesentests liegen nur vor, wenn man eine A-priori-Hypothese aufstellt (Bortz [3] Seite 378).

Signifikanz, p-Wert

Der p-Wert ist nicht die Signifikanz(-Zahl) und der p-Wert sagt nichts über die Stärke der Unterschiede aus. Der p-Wert ist auch kein Maß oder Beweis dafür, ob eine Unterschied besteht. Wenn p-Werte angegeben werden, liegt nicht automatisch ein Hypothesentest vor. Trotzdem verwende immer aus pragmatischen Gründen die Begriffe Signifikanz-Test und Hypothesentest sowie p-Werte und Signifikanz (fälschlicherweise) als gleichbedeutende Begriffe.

More Examples

Daten aufbereiten

Steko<- GetData("Steko.sav") %>% 
           Drop_NA(key) %>% 
           mutate(jahr = factor(jahr)) %>%
           Label( BMI = "Body-Mass-Index",
                  WHtR =  "Waist-Height-Ratio",
                  WHtR_1 ="Waist-Height-Ratio",
                  bildprof = "Bildungsprofil",
                  jahr = "Jahr")

data<- data.frame(
        g= gl(2, 8, labels = c("Control", "Treat")),
        a=rpois(16, 5),
        b=rep(1:2,8),
        c=rnorm(16))

data$d<-rep(c("a", "b", "c"), 6)[1:16]

upData2(data, labels = c(g="Gruppe", a="A", b="B"))

datadir <- GetData("
          names  labels    levels
            g     Gruppe    NULL
            a     A         numeric
            b     Sex       male;female
            d     DNA       a;b
            c     C         NULL")

upData2(data, datadir)



# oder auch

hkarz$Tzell <- cut(hkarz$tzell, 4)
 upData2(hkarz, data.frame(
   names=c("gruppe", "Tzell", "lai"),
   levels=c("gesund;krank", "Low;;Med;", "neg;pos")
))

Mit separate_multiple_choice() können Strings mit einem Trennzeichen zu einem Mehrfachantworten-Set aufgespalten werden. Erlaubt ist auch ein String mit Zahlen dabei ist aber der Wertebereich von 1 bis 11 begrenzt.

x <- c(123, 23, 456,3)
separate_multiple_choice(x,
                          label = c("Allgemein", "Paro", 
                          "Endo", "Prothetik",
                           "Oral chirurgie",
                          "Kieferorthopedie"))

Drop_NA(data)
Drop_case(data, 1:2)
Select_case(data, g=="Control")
transform2(data, a=a+1)

#library(tidyverse)
#library(stp25aggregate)
#library(stp25data)

hyper1<-hyper[, c("g","chol0","chol1","chol6","chol12")]
hyper_long<- Melt2(hyper1, id.vars=1)
aggregate( value~variable, hyper_long, mean)

hyper_long<-Melt2(chol0+chol1+chol6+chol12~g, 
                  hyper, "Zeit", "Cholesterin")
#-- Spread + aggragate
aggregate(Cholesterin~Zeit+g, hyper_long, mean) %>% 
  spread(g, Cholesterin) 

#-- das gleiche wie aggragate
hyper_long %>% group_by(Zeit, g) %>%
  summarise(Cholesterin=mean(Cholesterin)) %>%
  spread(g, Cholesterin)

#-- Gather das gleiche wie oben aber ohne die Labels
 hyper  %>% 
   tidyr::gather("time", "chol", chol0:chol12) %>%
   dplyr::select(g, time, chol) 

 #-- Recast2 das gleiche wie oben
 Recast2(chol0+chol1+chol6+chol12~g, hyper, mean, 
         "Zeit", "Cholesterin") %>%
   spread(g, Cholesterin)

 #-- Recast2 das gleiche wie oben nur ohne  spread
Recast2(chol0+chol1+chol6+chol12~g, hyper, mean, 
        formula=variable~g) 

 #-- und hier mit margins
Recast2(chol0+chol1+chol6+chol12~g, hyper, 
        mean,formula=variable~g, margins=TRUE)

Messmethoden Weibull-Modul

Beim Weibull-Modul handelt es sich um einen statistischen Wert, der unter anderem auch in Bezug auf die Materialermüdung von spröden Materialien verwendet wird. Dabei gilt: Ist der Weibull-Modul hoch, brechen Proben immer ab der gleichen Belastung. Ist der WeibullModul niedrig, ergeben sich Schwankungen in den Belastungswerten. Anmerkung bei Kunde 674 Nadine Wanda verwendet. Weitere Info unter: https://stats.stackexchange.com/questions/60511/weibull-distribution-parameters-k-and-c-for-wind-speed-data

#Beispiel: Sachs Seite 337
#Scheuerfestigkeit eines Garns
garn  <- c( 550,  760,  830,  890, 1100, 
           1150, 1200, 1350, 1400, 1600, 
           1700, 1750, 1800, 1850, 1850, 
           2200, 2400, 2850, 3200)

Weibull<- function(x){# empirische Verteilungsfunktion
x <- sort(x); n <- length(x); i <- rank(x)
 if(n < 50) F_t <- (i - 0.3) / (n+0.4)  
 else       F_t <- i/(n+1)

 data.frame(data=x, x=log(x),
      y =log(log(1/(1-F_t))))
}

data<-Weibull(garn) 
z <- lm(y ~ x, data)
  res<-round(cbind(shape=coef(z)[2],             
                   scale=exp(-(coef(z)[1]/coef(z)[2]))), 2) 
res


m<-(-(coef(z)[1]/coef(z)[2]))

#library(MASS)
fit <- MASS::fitdistr(garn, densfun="weibull", lower = 0)

library(car)

par(mfrow=c(1,2), lwd=2, font.axis=2, bty="n", ps=11,
    bg="white", las=1, cex=1.1)

plot(data$x, data$y, main="Weibull-Diagram",
     xlab="x=log(Garn)", ylab="y=log(log(1/(1-F)))", 
     xlim=c(5.9,8.5), ylim=c(-4, 2), axes = FALSE, pch=16, cex=1.0)
abline(z)
my.at<- signif(c(500, 1000, exp(m), 5000),4)
axis(1, at = log(my.at), labels = my.at)
axis(2)
i <- rank(data$x)
n <- length(data$x)
v.unten <- 1 / (((n-i+1)/i)*(qf(0.025, 2*(n-i+1), 2*i))+1)
v.oben  <- 1 - 1 / (1 + (i / (n-i+1)*qf(0.025, 2*i, 2*(n-i+1))))
lines(data$x, log(log(1/(1-v.oben))), lty=2)
lines(data$x, log(log(1/(1-v.unten))), lty=2)
abline(h=0, lty=3)
abline(v=m, lty=3)
points(m, 0, cex=5, col="red")



qqPlot(garn, distribution="weibull", main="qqPlot",
       scale=coef(fit)[1], shape=coef(fit)[2], las=1, pch=19)

Literatur



stp4/stp25stat documentation built on Sept. 17, 2021, 2:03 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close