Kapitel 7: Kapitel 7: Statistische Analysen produktiver Kompetenzen

Description

Das ist die Nutzerseite zum Kapitel 7, Statistische Analysen produktiver Kompetenzen, im Herausgeberband Large-Scale Assessment mit R: Methodische Grundlagen der <U+00F6>sterreichischen Bildungsstandard<U+00FC>berpr<U+00FC>fung. Im Abschnitt Details werden die im Kapitel verwendeten R-Syntaxen zur Unterst<U+00FC>tzung f<U+00FC>r Leser/innen kommentiert und dokumentiert. Im Abschnitt Examples werden die R-Syntaxen des Kapitels vollst<U+00E4>ndig wiedergegeben und gegebenenfalls erweitert.

Arguments

Details

Abschnitt 1: Beispieldatens<U+00E4>tze

Der zur Illustration verwendete Datensatz prodRat beinhaltet die beurteilten Schreibkompetenzen im Fach Englisch auf der 8. Schulstufe von 9836 Sch<U+00FC>ler/innen (idstud) die von insgesamt 41 Ratern (rater) beurteilt wurden. Die sechs Schreibaufgaben (aufgabe) wurden auf sechs Testhefte (th) aufgeteilt, wobei jede Aufgabe in genau zwei Testheften vorkommt.

Zur weiteren Analyse verwenden wir auch den Datensatz prodPRat mit sogenannten Pseudoratern.

F<U+00FC>r die Analyse von Varianzkomponenten mittels Linear Mixed Effects (LME) Modellen verwenden wir den urspr<U+00FC>nglichen Datensatz im Long Format (prodRatL).

Abschnitt 2: Beurteiler<U+00FC>bereinstimmung

Listing 1: Berechnen von H<U+00E4>ufigkeitstabellen

Hier werden die Datens<U+00E4>tze prodRat und prodPRat verwendet. Die R-Funktion apply() erm<U+00F6>glicht eine Anwendung einer beliebigen Funktion z.B. prop.table() <U+00FC>ber alle Zeilen (1) oder Spalten (2) eines data.frame.

library(irr) data(datenKapitel07) prodRat <- datenKapitel07$prodRat

# Items ausw<U+00E4>hlen items <- c("TA", "CC", "GR", "VO") # Tabelle erzeugen tab <- apply(prodRat[, items], 2, FUN=function(x){ prop.table(table(x))*100}) print(tab, digits = 2)

# Mittelwert der Ratings berechnen round(apply(prodRat[, items], 2, mean), 2)

Listing 2: Beurteiler<U+00FC>bereinstimmung berechnen

Wir verwenden den Datensatz mit Pseudoratern prodPRat. Die Analysen werden mit dem Paket irr durchgef<U+00FC>hrt.

prodRat <- datenKapitel07$prodRat

items <- c("TA", "CC", "GR", "VO") dfr <- data.frame(items, agree = NA, kappa = NA, wkappa = NA, korr = NA) for(i in 1:length(items)){ dat.i <- prodPRat[, grep(items[i], colnames(prodPRat))] dfr[i, "agree"] <- agree(dat.i, tolerance = 1)["value"] dfr[i, "kappa"] <- kappa2(dat.i)["value"] dfr[i, "wkappa"] <- kappa2(dat.i, weight = "squared")["value"] dfr[i, "korr"] <- cor(dat.i[,1], dat.i[,2]) dfr[i, "icc"] <- icc(dat.i, model = "twoway")["value"] } print(dfr, digits = 3)

Abschnitt 3: Skalierungsmodelle

Listing 1: Skalierungsmodell mit TAM

Der Funktion tam.mm.mfr() muss ein data.frame f<U+00FC>r die Facetten <U+00FC>bergeben werden. Zus<U+00E4>tzlich k<U+00F6>nnen Einstellungen in einer Liste f<U+00FC>r das Argument control = list() <U+00FC>bergeben werden. Hier verwenden wir die Einstellung xsi.start0 = 1, was dazu f<U+00FC>hrt, dass alle Startwerte auf 0 gesetzt werden. Mit fac.oldxsi = 0.1 setzen wir das Gewicht der Parameterwerte aus der vorigen Iteration etwas <U+00FC>ber 0. Damit kann der Algorithmus stabilisiert und Konvergenzprobleme (deviance increase) verhindert werden. Wir definieren noch increment.factor = 1.05 etwas <U+00FC>ber dem default-Wert von 1 um m<U+00F6>gliche Konvergenzprobleme abzufangen. Dieser Wert definiert das Ausma<U+00DF> der Abnahme des maximalen Zuwachs der Parameter pro Iteration (s. TAM-Hilfe).

Die Personenparameter werden mit der Funktion tam.wle() gesch<U+00E4>tzt.

Gibt man in der Funktion summary() das Argument file an, so wird der Output direkt in ein Textfile geschrieben.

set.seed(1234) library(TAM)

prodRat <- datenKapitel07$prodRat

# Rater-Facette definieren facets <- prodRat[, "rater", drop = FALSE]

# Response Daten definieren vars <- c("TA", "CC", "GR", "VO") resp <- prodRat[, vars]

# Personen-ID definieren pid <- prodRat$idstud

# Formel f<U+00FC>r Modell formulaA <- ~item*step+item*rater

# Modell berechnen mod <- tam.mml.mfr(resp = resp, facets = facets, formulaA = formulaA, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) summary(mod, file="TAM_MFRM")

# Personenparameter und Rohscores persons.mod <- tam.wle(mod) persons.mod$raw.score <- persons.mod$PersonScores / (persons.mod$N.items)

Listing 1b (Erg<U+00E4>nzung zum Buch): Skalierungsmodell mit TAM

Hier werden alle im Buch besprochenen Modelle berechnet und anschlie<U+00DF>end ein Modellvergleich durchgef<U+00FC>hrt.

f1 <- ~item * rater * step mod1 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f1, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) f2 <- ~item*step+item*rater mod2 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f2, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) f3 <- ~item * step + rater mod3 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f3, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) f4 <- ~item + step + rater mod4 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f4, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) mod4$xsi.facets IRT.compareModels(mod1, mod2, mod3, mod4)

Listing 1c (Erg<U+00E4>nzung zum Buch): Wright-Map

Mit dem Paket WrightMap k<U+00F6>nnen die Ergebnisse f<U+00FC>r die einzelnen Facetten dargestellt werden. Wir machen dies f<U+00FC>r Items und Rater.

library(WrightMap)

item.labs <- vars rater.labs <- unique(prodRat$rater) item.labs <- c(item.labs, rep(NA, length(rater.labs) - length(item.labs)))

pars <- mod$xsi.facets$xsi facet <- mod$xsi.facets$facet item.par <- pars[facet == "item"] rater.par <- pars[facet == "rater"] item_rat <- pars[facet == "item:rater"] len <- length(item_rat) item.long <- c(item.par, rep(NA, len - length(item.par))) rater.long <- c(rater.par, rep(NA, len - length(rater.par)))

wrightMap(persons.mod$theta, rbind(item.long, rater.long), label.items = c("Items", "Rater"), thr.lab.text = rbind(item.labs, rater.labs), axis.items = "", min.l=-3, max.l=3, axis.persons = "Personen")

Listing 2: Fit-Indices berechnen

Die Fit-Indices werden mit der Funktion msq.itemfitWLE f<U+00FC>r die Raterparameter und Itemparameter gesondert berechnet. Der Funktion muss ein Vektor mit Parameterbezeichnungen <U+00FC>bergeben werden so wie sie im Modell-Objekt vorkommen. Im Paket TAM gibt es noch die Funktion tam.fit(), diese basiert auf einer Simulation der individuellen Posterior-Verteilung. Die Funktion msq.itemfitWLE wertet dagegen die individuelle Posterior-Verteilung direkt aus (s. TAM-Hilfe f<U+00FC>r weitere Beispiele) und f<U+00FC>hrt keine Simulation durch.

# Infit/Outfit berechnen pseudo_items <- colnames(mod$resp) pss <- strsplit(pseudo_items , split="-") item_parm <- unlist(lapply(pss, FUN = function(ll){ll[1]})) rater_parm <- unlist(lapply(pss, FUN = function(ll){ll[2]}))

# Fit Items res.items <- msq.itemfitWLE(mod, item_parm) summary(res.items)

# Fit Rater res.rater <- msq.itemfitWLE(mod, rater_parm) summary(res.rater)

Listing 2a (Erg<U+00E4>nzung zum Buch): Abbildung Fit-Indices

# Abbildung: Histogramm, Rohscores

par(mfcol=c(1,2))

hist(persons.mod$theta, col="grey", breaks=40, main = "", xlab = "Theta (logits)", ylab = "H<U+00E4>ufigkeit") with(persons.mod, scatter.smooth(raw.score, theta, pch = 1, cex = .6, xlab = "Roscores", ylab = "Theta (logits)", lpars = list(col = "darkgrey", lwd = 2, lty = 1)))

# Abbildung: Fit-Statistik par(mfcol=c(1,2)) fitdat <- res.rater$fit_data fitdat$var <- factor(substr(fitdat$item, 1, 2)) boxplot(Outfit~var, data=fitdat, ylim=c(0,2), main="Outfit") boxplot(Infit~var, data=fitdat, ylim=c(0,2), main="Infit")

Listing 2b (Erg<U+00E4>nzung zum Buch): Korrelationen

Pearson und Spearman Korrelationskoeffizient wird zwischen Rohscores und Theta berechnet.

korr <- c(with(persons.mod, cor(raw.score, theta, method = "pearson")), with(persons.mod, cor(raw.score, theta, method = "spearman"))) print(korr)

Listing 3: Q3-Statistik berechnen

Die Q3-Statistik f<U+00FC>r lokale stochastische Unabh<U+00E4>ngigkeit wird mit der Funktion tam.modelfit() berechnet. Der Output enth<U+00E4>lt eine Vielzahl an Fit-Statistiken, f<U+00FC>r weitere Details sei hier auf die TAM-Hilfeseite verwiesen. Die adjustierte aQ3-Statistik berechnet sich aus den Q3-Werten abz<U+00FC>glich des Gesamtmittelwerts von allen Q3-Werten.

Mit tam.modelfit() werden Fit-Statistiken f<U+00FC>r alle Rater x Item Kombinationen berechnet. Diese werden im Code unten anschlie<U+00DF>end aggregiert um eine <U+00DC>bersicht zu erhalten. Dazu werden zuerst nur Paare gleicher Rater ausgew<U+00E4>hlt, somit wird die aggregierte Q3-Statistik nur Rater-spezifisch berechnet. Das Objekt rater.q3 beinhaltet eine Zeile pro Rater x Item Kombination. Kombinationen ergeben sich nur f<U+00FC>r einen Rater, nicht zwischen unterschiedlichen Ratern.

Anschlie<U+00DF>end kann man mit aggregate() separat <U+00FC>ber Rater und Kombinationen mitteln und diese als Dotplot darstellen (Paket lattice).

# Q3 Statistik mfit.q3 <- tam.modelfit(mod) rater.pairs <- mfit.q3$stat.itempair

# Nur Paare gleicher Rater w<U+00E4>hlen unique.rater <- which(substr(rater.pairs$item1, 4,12) == substr(rater.pairs$item2, 4,12)) rater.q3 <- rater.pairs[unique.rater, ]

# Spalten einf<U+00FC>gen: Rater, Kombinationen rater.q3$rater <- substr(rater.q3$item1, 4, 12) rater.q3 <- rater.q3[order(rater.q3$rater),] rater.q3$kombi <- as.factor(paste(substr(rater.q3$item1, 1, 2), substr(rater.q3$item2, 1, 2), sep="_"))

# Statistiken aggregieren: Rater, Kombinationen dfr.raterQ3 <- aggregate(rater.q3$aQ3, by = list(rater.q3$rater), mean) colnames(dfr.raterQ3) <- c("Rater", "Q3") dfr.itemsQ3 <- aggregate(rater.q3$aQ3, by = list(rater.q3$kombi), mean) colnames(dfr.itemsQ3) <- c("Items", "Q3") dfr.itemsQ3

Listing 3 (Erg<U+00E4>nzung zum Buch): Lattice Dotplot

library(lattice) library(grid)

# Lattice Dotplot mean.values <- aggregate(rater.q3$aQ3, list(rater.q3$kombi), mean)[["x"]] dotplot(aQ3~kombi, data=rater.q3, main="Q3-Statistik", ylab="Q3 (adjustiert)", col="darkgrey", panel = function(x,...){ panel.dotplot(x,...) panel.abline(h = 0, col.line = "grey", lty=3) grid.points(1:6, mean.values, pch=17) })

Abschnitt 4: Generalisierbarkeitstheorie

Listing 1: Varianzkomponenten mit lme4 berechnen

Mit der Funktion lmer() aus dem Paket lme4 sch<U+00E4>tzen wir die Varianzkomponenten. In der Formel definieren wir dabei die Facetten als random effects.

library(lme4)

prodRatL <- datenKapitel07$prodRatL

# Zusammenfassung der Ergebnisse summary(mod.vk)

Listing 1a (Erg<U+00E4>nzung zum Buch): Summary-Funktion f<U+00FC>r Varianzkomponenten

Wir generieren eine Funktion summary.VarComp(), die den Output des Modells mod.vk in einen ansprechenden data.frame schreibt. Hier werden auch die prozentualen Anteile der Varianzkomponenten berechnet.

# Helper-Function um die Varianzkomponenten zu extrahieren summary.VarComp <- function(mod){ var.c <- VarCorr(mod) var.c <- c(unlist(var.c) , attr(var.c , "sc")^2) names(var.c)[length(var.c)] <- "Residual" dfr1 <- data.frame(var.c) colnames(dfr1) <- "Varianz" dfr1 <- rbind(dfr1, colSums(dfr1)) rownames(dfr1)[nrow(dfr1)] <- "Total" dfr1$prop.Varianz <- 100 * (dfr1$Varianz / dfr1$Varianz[nrow(dfr1)]) dfr1 <- round(dfr1,2) return(dfr1) } summary.VarComp(mod.vk)

Listing 2: Berechnung des G-Koeffizienten

Den G-Koeffizienten berechnen wir nach der Formel im Buch.

vk <- summary.VarComp(mod.vk) n.p <- length(unique(prodRatL$idstud)) # Anzahl Sch<U+00FC>ler n.i <- 4 # Anzahl Items n.r <- c(1,2,5) # Anzahl Rater

# Varianzkomponenten extrahieren sig2.p <- vk["idstud", "Varianz"] sig2.i <- vk["item", "Varianz"] sig2.r <- vk["rater", "Varianz"] sig2.ri <- vk["rater:item", "Varianz"] sig2.pr <- vk["idstud:rater", "Varianz"] sig2.pi <- vk["idstud:item", "Varianz"] sig2.pir <- vk["Residual", "Varianz"]

# Fehlervarianz berechnen sig2.delta <- sig2.pi/n.i + sig2.pr/n.r + sig2.pir/(n.i*n.r)

# G-Koeffizient berechnen g.koeff <- sig2.p / (sig2.p + sig2.delta) print(data.frame(n.r, g.koeff), digits = 3)

Listing 2a (Erg<U+00E4>nzung zum Buch): Phi-Koeffizient berechnen

sig2.D <- sig2.r/n.r + sig2.i/n.i + sig2.pi/n.i + sig2.pr/n.r + sig2.ri/(n.i*n.r) + sig2.pir/(n.i*n.r) phi.koeff <- sig2.p / (sig2.p + sig2.D) print(data.frame(n.r, phi.koeff), digits = 3)

# Konfidenzintervalle 1.96*sqrt(sig2.D)

Listing 2c (Erg<U+00E4>nzung zum Buch): Variable Rateranzahl

F<U+00FC>r eine variable Rateranzahl (hier 1 bis 10 Rater) werden die G-Koeffizienten berechnet.

n.i <- 4 # Anzahl Items dn.r <- seq(1,10) # 1 bis 10 m<U+00F6>gliche Rater delta.i <- sig2.pi/n.i + sig2.pr/dn.r + sig2.pir/(n.i*dn.r) g.koeff <- sig2.p / (sig2.p + delta.i) names(g.koeff) <- paste("nR", dn.r, sep="_") print(g.koeff[1:4])

plot(g.koeff, type = "b", pch = 19, lwd = 2, bty = "n", main = "G-Koeffizient: Raters", ylab = "G-Koeffizient", xlab = "Anzahl Raters", xlim = c(0,10)) abline(v=2, col="darkgrey")

Abschnitt 5: Strukturgleichungsmodelle

In R setzen wir das Struktugleichungsmodell mit dem Paket lavaan um. Das Modell wird als Textvariable definiert, welche anschlie<U+00DF>end der Funktion sem() <U+00FC>bergeben wird. Latente Variablen im Messmodell werden in lavaan mit der Form latente Variable =~ manifeste Variable(n) definiert, die Ladungen werden dabei auf den Wert 1 fixiert, was mittels der Multiplikation der Variable mit dem Wert 1 umgesetzt werden kann (z.B. 1*Variable). Varianzen und Kovarianzen werden mit Variable ~~ Variable gebildet, wobei hier die Multiplikation mit einem Label einerseits den berechneten Parameter benennt, andererseits, bei mehrmaligem Auftreten des Labels, Parametersch<U+00E4>tzungen von verschiedenen Variablen restringiert bzw. gleichstellt (z.B. wird f<U+00FC>r die Within-Varianz von TA <U+00FC>ber beide Rater nur ein Parameter gesch<U+00E4>tzt, n<U+00E4>mlich Vta_R12). Die ICC wird f<U+00FC>r jede Dimension separat direkt im Modell spezifiziert, dies geschieht durch abgeleitete Variablen mit der Schreibweise Variable := Berechnung. Die Modellspezifikation und der Aufruf der Funktion sem() ist wie folgt definiert:

Listing 1 (mit Erg<U+00E4>nzung zum Buch): SEM

library(lavaan)

prodPRat <- datenKapitel07$prodPRat

# SEM Modell definieren lv.mod <- " # Messmodell TA =~ 1*TA_R1 + 1*TA_R2 CC =~ 1*CC_R1 + 1*CC_R2 GR =~ 1*GR_R1 + 1*GR_R2 VO =~ 1*VO_R1 + 1*VO_R2

# Varianz Between (Personen) TA ~~ Vta * TA CC ~~ Vcc * CC GR ~~ Vgr * GR VO ~~ Vvo * VO

# Varianz Within (Rater X Personen) TA_R1 ~~ Vta_R12 * TA_R1 TA_R2 ~~ Vta_R12 * TA_R2 CC_R1 ~~ Vcc_R12 * CC_R1 CC_R2 ~~ Vcc_R12 * CC_R2 GR_R1 ~~ Vgr_R12 * GR_R1 GR_R2 ~~ Vgr_R12 * GR_R2 VO_R1 ~~ Vvo_R12 * VO_R1 VO_R2 ~~ Vvo_R12 * VO_R2

# Kovarianz Within TA_R1 ~~ Kta_cc * CC_R1 TA_R2 ~~ Kta_cc * CC_R2 TA_R1 ~~ Kta_gr * GR_R1 TA_R2 ~~ Kta_gr * GR_R2 TA_R1 ~~ Kta_vo * VO_R1 TA_R2 ~~ Kta_vo * VO_R2 CC_R1 ~~ Kcc_gr * GR_R1 CC_R2 ~~ Kcc_gr * GR_R2 CC_R1 ~~ Kcc_vo * VO_R1 CC_R2 ~~ Kcc_vo * VO_R2 GR_R1 ~~ Kgr_vo * VO_R1 GR_R2 ~~ Kgr_vo * VO_R2

# ICC berechnen icc_ta := Vta / (Vta + Vta_R12) icc_cc := Vcc / (Vcc + Vcc_R12) icc_gr := Vgr / (Vgr + Vgr_R12) icc_vo := Vvo / (Vvo + Vvo_R12) " # Sch<U+00E4>tzung des Modells mod1 <- sem(lv.mod, data = prodPRat) summary(mod1, standardized = TRUE)

# Inspektion der Ergebnisse show(mod1) fitted(mod1) inspect(mod1,"cov.lv") inspect(mod1, "free")

Listing 2: Kompakte Darstellung der Ergebnisse

parameterEstimates(mod1, ci = FALSE, standardized = TRUE)

Listing 2a (Erg<U+00E4>nzung zum Buch): Schreibe Ergebnisse in Latex-Tabelle

library(xtable)

xtable(parameterEstimates(mod1, ci = FALSE, standardized = TRUE), digits = 3)

Abschnitt 7: <U+00DC>bungen

<U+00DC>bung 1: MFRM M3 und M4 umsetzen und Vergleichen

Wir setzen die Modelle separat in TAM um und lassen uns mit summary() die Ergebnisse anzeigen. Einen direkten Zugriff auf die gesch<U+00E4>tzen Parameter bekommt man mit mod$xsi.facets. Dabei sieht man, dass im Modell 4 keine generalized items gebildet werden, da hier kein Interaktionsterm vorkommt. Den Modellvergleich machen wir mit IRT.compareModels(mod3, mod4). Modell 3 weist hier kleinere AIC-Werte auf und scheint etwas besser auf die Daten zu passen als Modell 4. Dies zeigt auch der Likelihood-Ratio Test, demnach sich durch das Hinzuf<U+00FC>gen von Parametern die Modellpassung verbessert.

library(TAM) prodRatEx <- datenKapitel07$prodRatEx

# Rater-Facette definieren facets <- prodRatEx[, "rater", drop = FALSE]

# Response Daten definieren vars <- c("TA", "CC", "GR", "VO") resp <- prodRatEx[, vars]

# Personen-ID definieren pid <- prodRatEx$idstud

# Modell 3 f3 <- ~item * step + rater mod3 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f3, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) # Modell 4 f4 <- ~item + step + rater mod4 <- tam.mml.mfr(resp = resp, facets = facets, formulaA = f4, pid = pid, control=list(xsi.start0 = 1, fac.oldxsi = 0.1, increment.factor = 1.05)) summary(mod3, file = "TAM_MFRM") summary(mod4, file = "TAM_MFRM")

mod3$xsi.facets mod4$xsi.facets

IRT.compareModels(mod3, mod4)

$IC Model loglike Deviance Npars Nobs AIC BIC AIC3 AICc CAIC 1 mod3 -60795.35 121590.7 69 9748 121728.7 122224.5 121797.7 121729.7 122293.5 2 mod4 -61041.47 122082.9 51 9748 122184.9 122551.4 122235.9 122185.5 122602.4

$LRtest Model1 Model2 Chi2 df p 1 mod4 mod3 492.2264 18 0

<U+00DC>bung 2: Varianzkomponentenmodell

Das Varianzkomponentenmodell setzen wir f<U+00FC>r die short prompts nach den Vorgaben im Buchkapitel um. Dabei ver<U+00E4>ndern wir die Anzahl der m<U+00F6>glichen Rater durch n.r <- c(2,10,15). Der Phi-Koeffizient kann laut Gleichung 6.9 und 6.10 berechnet werden. Die Ergebnisse zeigen einen prozentuellen Anteil der Interaktion Person und Rater von ca. 15%, dieser scheint auf Halo-Effekte hinzuweisen.

library(lme4) prodRatLEx <- datenKapitel07$prodRatLEx

# Zusammenfassung der Ergebnisse summary(mod.vk) print(vk <- summary.VarComp(mod.vk))

Varianz prop.Varianz idstud:item 0.10 2.45 idstud:rater 0.64 15.21 idstud 2.88 67.94 rater:item 0.01 0.22 rater 0.19 4.39 item 0.00 0.02 Residual 0.41 9.78 Total 4.24 100.00

# Ver<U+00E4>ndern der Rateranzahl n.p <- length(unique(prodRatLEx$idstud)) # Anzahl Sch<U+00FC>ler n.i <- 4 # Anzahl Items n.r <- c(2,10,15) # Anzahl Rater

# Phi-Koeffizient berechnen sig2.D <- sig2.r/n.r + sig2.i/n.i + sig2.pi/n.i + sig2.pr/n.r + sig2.ri/(n.i*n.r) + sig2.pir/(n.i*n.r) phi.koeff <- sig2.p / (sig2.p + sig2.D) print(data.frame(n.r, phi.koeff), digits = 3)