stp4/stp25APA2
Resultate Formatiern

Package index
Search the stp4/stp25APA2 package
Vignettes
Functions
674
Source code
75
Man pages
39
Home
/
GitHub
/
stp4/stp25APA2
/

In stp4/stp25APA2: Resultate Formatiern

Link:

#vignette: |
 # %\VignetteIndexEntry{Statistische Methoden} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8}#
#setwd("C:/Users/wpete/Dropbox/3_Forschung/R-Project/stp25APA2/vignettes")
knitr::opts_chunk$set(echo = TRUE)
#owd = setwd('vignettes')

Einleitung

Die von mir verwendeten statistischen Methoden basieren vor allem auf den Empfehlungen von Bortz [4] sowie Sachs [7]. Die Darstellung der Ergebnisse entspricht wissenschaftlicher Vorgaben, insbesondere halte ich mich bei den Tabellen und Grafiken sowie der Darstellung statistischer Kennzahlen an die Vorgaben von APA-Style[2]. (APA-Style ist im Kontext sozialwissenschaftlicher Forschung quasi der Gold-Standard hinsichtlich des Berichtens von Ergebnissen.)

 #getwd()
library(stpvers)
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 (Lienert [5] 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 (Bortz [4] 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 Prop_Test2().

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],]

Prop_Test2(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)
#  APA2(~g+g2, data, type="freq.ci")
#  APA2(g~g2, data, type="freq.ci")

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 ass 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))
# ANOVA 
 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))

require(broom)
require(nlme)
res <- summary(gls(m1 ~ geschl, varana,  weights=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.

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


xtab <- xtabs(~ gruppe+LAI, hkarz)
APA2(xtab, caption="Harnblasenkarzinom", test=FALSE)
APA2(xtab, type="sens", test=TRUE, caption = "type=sens")
APA2(xtab, type="sens", caption = "geht nur mit teat=TRUE +  type=sens")
APA2(xtabs(~ gruppe+Tzell, hkarz), caption="APA_Xtabs: 2x3 Tabelle", test=FALSE)
APA2(xtabs(~ gruppe+LAI+Tzell, hkarz), caption="APA_Xtabs: 2x2x3 Tabelle", test=FALSE)


APA2(xtab, include.total.columns=TRUE, caption = "include.total.columns")
APA2(xtab, include.total.sub=TRUE, caption = "include.total.sub")




xtab <- xtabs(~ gruppe+Tzell, hkarz)
APA2(xtab, test=FALSE, caption="APA2: 2x3 Tabelle")
APA_Xtabs(xtab, caption="APA_Xtabs: 2x3 Tabelle")

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.

xtab <- xtabs(~ gruppe+LAI, hkarz)
fit1<- glm(gruppe~lai, hkarz, family = binomial)

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

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}$$

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.

J = max(Sensitivity + Specificity - 1)

Sachs 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))

# require(pROC)
# data %>% Tabelle(Blutzucker, by=~Gruppe) 

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)
# windows(8,8)
# plot(roc_curve, print.thres = "best",
#      print.auc=TRUE)
# abline(v = 1, lty = 2)
# abline(h = 1, lty = 2)

Weitere Info zu ROC findet sich unter https://www.hranalytics101.com/how-to-assess-model-accuracy-the-basics/

fit1<- glm(gruppe~lai, hkarz, family = binomial)
x1 <- Klassifikation(fit1)
x1$statistic[c(1,7,8),]

roc_curve <- roc(x1$response, x1$predictor)
auc(roc_curve)
#plot(roc_curve, print.thres = "best", print.auc=TRUE)
#abline(v = 1, lty = 2)
#abline(h = 1, lty = 2)
#text(.90, .97, labels = "Ideal Model")
#points(1,1, pch = "O", cex = 0.5)

# Compute confidence intervals at every point of the curve (looks nicer when plotting)
ciobj <- ci.se(rocobj, # CI of sensitivity...
               specificities = seq(0,100,5)) # ...at every 5% of Sensitivity from 0-100

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)

### Not run: 
# The latter used Delong's test. To use bootstrap test:
#roc.test(roc1, roc2, method="bootstrap")
# Increase boot.n for a more precise p-value:
#roc.test(roc1, roc2, method="bootstrap", boot.n=10000)

#ciobj <- ci.se(roc2)
#plot(ciobj, type = "shape", col="#D3D3D3", alpha = .5)

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

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")))

 APA_Correlation(~a+b+c, DF, caption="~a+b+c")

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. Literatur: Bortz, J. & Doering, N. (2006). Forschungsmethoden und Evaluation fuer Human-und Sozialwissenschaftler (4. Auflage). Berlin: Springer. Seite 155

DF <-structure(list(
        Geschlecht = structure(c(1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L),
                               .Label = c("Maennlich", "Weiblich"), class = "factor"),
        Alter = structure(c(2L, 4L, 2L, 4L, 2L, 2L, 2L, 3L, 3L, 2L, 1L, 1L, 3L, 4L, 4L, 4L, 2L, 1L, 2L, 1L, 4L, 4L, 3L, 4L, 2L, 2L, 1L, 4L, 4L, 3L, 3L, 3L, 3L, 2L, 3L, 4L, 3L, 3L, 1L, 3L, 1L, 1L, 2L, 1L, 1L, 4L, 3L, 1L, 4L, 2L, 2L, 1L, 3L, 3L, 2L, 3L, 4L, 4L, 1L, 2L, 3L, 2L, 1L, 2L, 1L, 2L, 3L),
                          .Label = c("20 - 29", "30 - 39", "40 - 49", "50 - 59"), class = "factor"),
        Konsum = structure(c(1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 3L, 1L, 1L, 1L, 2L, 3L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 3L, 1L, 2L, 2L, 3L, 2L),
                                .Label = c("weniger als 3 T.", "3 bis 6 T.", "mehr als 6 T."), class = "factor"),
        Kaffeeform = structure(c(3L, 1L, 3L, 2L, 3L, 3L, 3L, 1L, 3L, 1L, 3L, 3L, 1L, 2L, 3L, 3L, 3L, 3L, 1L, 3L, 3L, 2L, 2L, 3L, 3L, 3L, 3L, 2L, 3L, 1L, 2L, 2L, 2L, 1L, 3L, 3L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 3L, 3L, 3L, 2L, 3L, 3L, 2L, 3L, 3L, 3L, 3L, 3L, 2L, 3L, 3L, 2L, 1L, 3L, 3L, 3L, 2L, 2L),
                               .Label = c("Espresso", "Filterkaffee", "Milchkaffee"), class = "factor"),
        FavA = structure(c(3L, 1L, 2L, 1L, 3L, 3L, 4L, 1L, 2L, 2L, 1L, 1L, 1L, 4L, 3L, 4L, 3L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 3L, 4L, 1L, 1L, 4L, 4L, 1L, 1L, 1L, 2L, 1L, 2L, 4L, 3L, 2L, 4L, 1L, 1L, 2L, 2L, 2L, 4L, 2L, 2L, 2L, 1L, 3L, 2L, 4L, 2L, 4L, 1L, 4L, 4L, 2L, 1L, 1L, 4L, 2L, 1L, 3L, 2L, 3L),
                         .Label = c("Cubanischer Arabica Filter", "Cubanischer Arabica Kaltextrakt", "Dallmayr Prodomo Kaltextrakt", "Dallmayr Prodomo Filter"), class = "factor"),
        FavB = structure(c(4L, 2L, 1L, 3L, 2L, 1L, 3L, 2L, 1L, 1L, 4L, 4L, 2L, 2L, 2L, 2L, 4L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 4L, 2L, 2L, 3L, 4L, 4L, 1L, 3L, 1L, 2L, 4L, 4L, 1L, 3L, 3L, 1L, 3L, 1L, 1L, 1L, 3L, 1L, 2L, 2L, 1L, 3L, 3L, 3L, 2L, 2L, 3L, 3L, 2L, 4L, 1L, 1L, 2L, 2L, 1L, 2L),
                         .Label = c("Cubanischer Arabica Filter", "Cubanischer Arabica Kaltextrakt", "Dallmayr Prodomo Kaltextrakt", "Dallmayr Prodomo Filter"), class = "factor"),
        FavC = structure(c(2L, 3L, 3L, 4L, 1L, 2L, 1L, 4L, 4L, 3L, 2L, 3L, 3L, 3L, 4L, 3L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 3L, 3L, 3L, 3L, 1L, 4L, 3L, 2L, 3L, 2L, 3L, 1L, 2L, 3L, 2L, 4L, 4L, 4L, 4L, 3L, 2L, 3L, 1L, 3L, 4L, 4L, 3L, 2L, 1L, 1L, 3L, 3L, 2L, 1L, 4L, 2L, 3L, 3L, 4L, 1L, 3L, 1L),
                         .Label = c("Cubanischer Arabica Filter", "Cubanischer Arabica Kaltextrakt", "Dallmayr Prodomo Kaltextrakt", "Dallmayr Prodomo Filter"), class = "factor"),
        FavD = structure(c(1L, 4L, 4L, 2L, 4L, 4L, 2L, 3L, 3L, 4L, 3L, 2L, 4L, 1L, 1L, 1L, 1L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L, 1L, 4L, 2L, 1L, 3L, 2L, 2L, 3L, 4L, 4L, 4L, 3L, 1L, 1L, 3L, 2L, 2L, 3L, 1L, 4L, 3L, 4L, 4L, 4L, 3L, 1L, 4L, 1L, 4L, 2L, 4L, 1L, 1L, 4L, 3L, 3L, 2L, 4L, 3L, 4L, 4L, 4L),
                         .Label = c("Cubanischer Arabica Filter", "Cubanischer Arabica Kaltextrakt", "Dallmayr Prodomo Kaltextrakt", "Dallmayr Prodomo Filter"), class = "factor")),
        .Names = c("Geschlecht", "Alter", "Konsum", "Kaffeeform", "FavA", "FavB", "FavC", "FavD"), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L), class = "data.frame")

ans <- Rangreihe(~FavA+FavB+FavC+FavD, DF )
APA2(ans, caption="Alle")

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

Understanding Bland Altman analysis, Davide Giavarina, Biochemia medica 2015;25(2) 141-51 doi: 10.11613/BM.2015.015

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")))
)

MetComp2(~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)

MetComp2(~A+B, Botulinum)
# x <- BlandAltman(~A+E, DF)
# x %>% Output( )

 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 <- BlandAltman(~A+B+E, Giavarina)
# x %>% Output("BA-Analyse der Messwertreihe")
plot(x)

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

x <- BlandAltman(A+E+B~group, DF)
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 = round( A + rnorm(n,5,10) + (100-A/10) )
                E = A + ifelse(A<cutA, A/5, -A/5 )+ rnorm(n, 0, 10)
)

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

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

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

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

Regressionsanalyse

Regressionsanalysen sind eine Gruppe von statistische Verfahren, die Beziehungen zwischen Variablen modellieren. 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 [12]). Unstandardisierte Effektgrößen sind Differenzen von Gruppenmittelwerten (raw mean difference) und unstandardisierte Regressionskoeffizienten. Der p-Wert ist keine Effektstärke. Kreuztabelle: Cohen’s w, φ, Cramer’s V, C ANOVA/Regression: f, η², η²part und d, Δ, Hedge’s g (bei Vergleich der Stichproben) sowie Odds Ratio, Risk Ratio, r, R²

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)

α-Adjustierungen

R bietet α-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)

#  ANOVA
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||||

# a. 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) %>%
  APA_Table(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"
)

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

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)
APA_Table(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.

 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("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
 lmtest::lrtest(fit1) %>% fix_format() %>% Output("LR-Test")
# Fuer F-Test :  lmtest::waldtest(fit1)

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.

Klassifikation2(fit1)

fit1 <- glm(gruppe~tzell+factor(lai), hkarz, family = binomial)

x<- APA2(fit1, include.odds=TRUE )

Nagelkerke<- R2(fit1)[3]




# Interpreztation

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$OR[2], 
 "zu. Ist die  T-Zelltypisierung positiv so nimmt die  relative Wahrscheinlichkeit um OR= ", x$OR[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

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))


#xyplot(p~t.zell)

Poisson Regression

Poisson Regression für Zahldaten 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)")

Kaplan-Meier

Zum Vergleich von zwei Therapien werden häufig Überlebenszeiten herangezogen. Diese sollten https://www.uni-kiel.de/medinfo/lehre/seminare/methodik/Dtsch%20Arztebl%2015%20%C3%9Cberlebenszeitanalyse.pdf

Survival analysis

Beispiel MagenkarzinomBuehl

Die Datei mkarz.sav Seite 553 ist ein Datensatz mit 106 Patientan mit Magenkarzinom über einen Zeitraum von 5 Jahren

library(survival)
mkarz <- GetData("C:/Users/wpete/Dropbox/3_Forschung/1 Statistik/BspDaten/SPSS/_Buehl/MKARZ.SAV")
mkarz$status<- ifelse(mkarz$status=="tot", 1, 0)
mkarz %>% Tabelle2(survive="median", status, lkb)

# Kaplan-Meier estimator without grouping

m0 <- Surv(survive, status) ~ 1
res0<- survfit(m0, mkarz)
APA2(res0)

APA2(summary(res0, times= c(5, 10, 20, 60)),
     percent=TRUE,
     #Statistik Anfordern und ander Schreibweise
     include=c( time ="time", n.risk ="n.risk", 
                n.event ="n.event", surv = "survival",
                lower = "lower 95% CI",upper ="upper 95% CI"),
     caption="Kaplan-Meier" )

m1 <- Surv(survive, status) ~ lkb
res1<- survfit(m1, mkarz)
fit1<- coxph(m1, mkarz)
logrank1<- survdiff(m1, mkarz)
APA2(res1, caption="Kaplan-Meier")
APA2(logrank1)
APA2(coxph(m1,mkarz))

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.

Heteroskedasticity Breusch-Pagan test against heteroskedasticity.

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

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

VIF variance inflation factor. VIF values over 5 are troubling, should probably investigate anything over 2.5.

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)

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
ICC2(sf)

Explorative Faktorenanalyse (EFA)

Die EFA ist eine Gruppe von Verfahren für die Analyse unbekannten korrelativen Strukturen. Zu diesen Methoden zählen die Faktorenanalyse (FA) und die Hauptkomponenten-Analyse (PCA) . Ziel der Faktorenanalyse ist aus einer Vielzahl von Variablen unabhängige Beschreibung und Erklärungsvariablen zu extrahieren. Die Abgrenzung von FA zu PCA ist je nach Lehrbuch und verwendeter Software fließend. So spricht [15] dass die PCA ein Spezialfall der FA ist dagegen grenzt [16] die FA stärker von der PCA ab. Verwirrung schafft zudem SPSS welches keinen “Button” für PCA besitzt. Aus diesen Gründen verwende ich hier die Begriffe Faktorenanalyse (FA) und PCA als synonym.

Der Unterschied ist, das bei der PAC die Anzahl der Faktoren sich aus der Berechnung ergeben (Anwendungen sind sozialwissenschaftliche Fragestellungen) und bei der (FA) die Faktoren vor der Berechnung bekannt sein müssen (Anwendungen sind biologische Fragestellungen bei denen bestimmte Gradienten wie Temperatur, pH Wert bekannte Komponenten darstellen).

Interpretation FA Eine Auswertung die sich an den Tabellen Output von SPSS orientiert, kann mit der library library(psych) und den Funktionen principal() und fa() bewerkstelligt werden. Der Output entspricht dem was SPSS unter Faktorenanalyse berechnet. Ein anderer Zugang sind die Funktionen princomp() und prcomp() oder die library library(vegan) mit der Funktion rda().

Faktorladung = Korrelation zwischen Variable und Faktor Kommunalität h² = Ausmaß an Varianz der Variable Eigenwert = Varianz des Faktors

Kaiser-Gutman-Kriterium = nur solche Faktoren sind zu interpretieren deren Eigenwert größer 1 ist (Scree-Test.

  # APA2( ~., fkv, test=TRUE)
  # library(arm)
  # windows(5,5)
  # corrplot(fkv, abs=TRUE, n.col.legend=7)

 fit1<- Principal(fkv, 5, cut=.35)
 names(fit1$Loadings ) <- c("Item", "nr", 
                           "PC_1", "PC_2", "PC_3", "PC_4", "PC_5", "h2"  )
 fit1$Loadings %>% Output()

 fit1$Eigenvalue %>% Output()
 fit1$Test %>% Output()

Test der Eignung für die PCA Kaiser-Meyer-Olkin (KMO) KMO>0.80 ist gut. (KMO ist die partiellen Korrelationen zwischen den Itempaaren) Der Bartlett Test testet die Spherizität, überprüft die Nullhypothese, ob die Korrelationsmatrix eine Identitätsmatrix ist(gut ist wenn p < .050). 

 fit1$KMO %>% Output()

Faustregel nach Bortz:

unter n<60 ist keine Faktorenanalyse durchführbar, ausreichend sind n=100 die kann Faktorenstruktur interpretiert werden... - wenn auf einem Faktor mindestens vier Variablen Ladungen h>0.60 aufweisen
- oder wenn für zehn oder mehr Variablen Ladungen h>0.40 vorliegen - oder wenn weniger als zehn Variablen eine h>0.40 aufweisen und die Stichprobe n>300 Probanden enthält

Beispiel FKV

Freiburger Fragebogen zur Krankheitsverarbeitung (Buehl Seite 473)

Quelle:http://statistikguru.de/spss/hauptkomponentenanalyse/auswerten-und-berichten.html Es wurde eine Hauptkomponentenanalyse durchgeführt, um die wichtigsten, unabhängigen Faktoren zu extrahieren. Das Kaiser-Meyer-Olkin-Kriterium war 0.71 und der Bartlett-Test hoch signifikant (p < .001), was eine ausreichend hohe Korrelation zwischen Items darstellt, um eine Hauptkomponentenanalyse durchzuführen. Nur Faktoren mit Eigenwerten ≥ 1 wurden in Betracht gezogen (Guttman, 1954; Kaiser, 1960). Eine Überprüfung des Kaiser–Kriteriums und Scree-Plots rechtfertigte die Extraktion von zwei Faktoren, jeweils mit Eigenwerten über 1, die eine Gesamtvarianz von 43% aufklären. Unter den Lösungen lieferte die Varimax rotierte zweifaktor-Lösung die Lösung, die am besten zu interpretieren war, bei der die meisten Items nur auf einen der beiden Faktoren hohe Ladungen zeigten. Die Faktoren konnten den theoretischen Begriffen Verarbeitung, Coping, Vertrauen, Religion, Distanz und Verarbeitung zugeordnet werden

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).

SEM

Backhaus Multivariate Analysemethoden 11 AuflageSeite 383

Moosbrugger, Kelava 2012 Testtheorie 2. Auflage Seite 339

Fit-Mass | Guter Fit | Akzeptabler Fit -------- | --------- | --------------- Chi2/df | .00-2.00 | 2.01-3.00 RMSEA | .000-.050 | 051.-.080 CFI | .970-1.00 | .950-.969 NFI | .950-1.00 | .900-.949

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.

Literatur

[1] Achim Bühl, (2014), SPSS 22 Einführung in die moderne Datenanalyse, 14. aktualisierte Auflage, Pearson

[2] APA, 2009, Publication Manual of the American Psychological Association

[3] Daniel Wollschläger (2012), Grundlagen der Datenanalyse mit R: Eine anwendungsorientierte Einführung 2. Aufl., Heidelberg: Springer

[4] Jürgen Bortz, Nicola Döring, (2006), Forschungsmethoden und Evaluation, Heidelberg: Springer

[5] Jürgen Bortz, Christof Schuster, (2010), Statistik für Human- und Sozialwissenschaftler, 7. Aufl., Heidelberg: Springer [6] John Fox, Sanford Weisberg, (2011), An R Companion to Applied Regression, Second Edition, Sage

[7] Lothar Sachs, Jürgen Hedderich, (2006), Angewandte Statistik, 12.Aufl. Heidelberg: Springer

[9] Data Analysis Using Regression and Multilevel/Hierarchical Models; Cambridge;2009; Gelman, Hill

[10] Torsten Hothorn, Brian S., Everitt, (2009), A Handbook of Statistical Analyses Using R, Second Edition, Chapman and Hall/CRC

[11] Frank Harrell, PhD An Introduction to S and The Hmisc and Design Libraries https://cran.r-project.org/doc/contrib/Alzola+Harrell-Hmisc-Design-Intro.pdf (18-02-2016)

[12] Hemmerich MatheGuru http://matheguru.com/stochastik/effektstaerke.html (18-02-2016)

[13] R Core Team (2015). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.

[14] "Frank E Harrell Jr f.harrell@vanderbilt.edu, (2017), Hmisc: Harrell Miscellaneous, URL https://CRAN.R-project.org/package=Hmisc"

Literatur [15] Stoyan, D., Stoyan, H., Jansen, U. Umweltstatistik. Teubner, Stuttgart, Leipzig, 1997.

[16] Stier, W. Empirische Forschungsmethoden. 2. Au age, Springer, Berlin, 1999.

[17] Backhaus, K., Erichson, B., Plinke, W. & Weiber, R. Multivariate Analysemethoden. 12. Au age, Springer, Berlin, 2008.

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]))
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")


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

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


stp4/stp25APA2 documentation built on May 24, 2019, 9:59 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