• Program R

• RStudio

• Materiały

Metoda naukowa

Labirynt

• niezależność

• rozkład normalny

• podobna wariancja (zmienność)

• jednostajny

• Bernoulliego

• Poissona

• studenta

• \(\chi^2\)

• …

jednostajny = seq(1:100)
normalny = rnorm(100, mean = 50, sd = 1)
par(mfrow=c(1,2))
qqnorm(jednostajny)
qqnorm(normalny)

jednostajny = seq(1:100)
shapiro.test(jednostajny)

## 
##  Shapiro-Wilk normality test
## 
## data:  jednostajny
## W = 0.95472, p-value = 0.001722

normalny = rnorm(100, mean = 50, sd = 1)
shapiro.test(normalny)

## 
##  Shapiro-Wilk normality test
## 
## data:  normalny
## W = 0.98377, p-value = 0.258

data("PlantGrowth")

boxplot(weight ~ group, data = PlantGrowth,
        ylab = "Sucha masa roslin")

data("PlantGrowth")

bartlett.test(weight ~ group, data = PlantGrowth)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  weight by group
## Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371

Labirynt

• Porównujemy 2 hodowle A. thaliana rosnące przy niskim lub wysokim natężeniu światła przez 21 dni.

• Po upływie 3 tygodni zbieramy liście, ekstrahujemy barwniki i porówujemy zawartość chlorofilu całkowitego.

• Czy warunki oświetlenia mają wpływ na zawartość chlorofili w liściach ?

• \(H_{0}: brak \ istotnych \ różnic \ pomiędzy \ średnimi\)

• Test t (Studenta)

\(\ t=\frac{\overline{x}_{1}-\overline{x}_{2}}{\sqrt{\frac{S_{1}^{2}}{n_{1}}+\frac{S_{2}^{2}}{n_{1}}}}\)

• Jeżeli wartość t jest większa od wartości krytycznej odczytanej z tablic odrzucamy \(H_{0}\)!

\(\ t=\frac{\overline{x}_{LL}-\overline{x}_{HL}}{\sqrt{\frac{S_{LL}^{2}}{n_{LL}}+\frac{S_{HL}^{2}}{n_{HL}}}}\)

\(\ \overline{x}_{LL}\) - średnia zawartość chlorofilu w warunkach LL

\(\ S_{LL}\) - odchylenie standardowe zawartości chlorofilu w warunkach LL

\(\ n_{LL}\) - liczebność roślin w warunkach LL

t.test(LL, HL)

1. Przypisz 2 grupom roślin HL i LL zawartości chlorofili

HL = c(8.2, 6.5, 9, 8.7, 8.13)
LL = c(10.5, 9.5, 10, 7.91, 9.2)

1. Narysuj wykres porównujący grupy

boxplot(HL, LL)

1. Sprawdź testem t czy masz rację

t.test(HL, LL)

## 
##  Welch Two Sample t-test
## 
## data:  HL and LL
## t = -2.1379, df = 7.9987, p-value = 0.065
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.7355354  0.1035354
## sample estimates:
## mean of x mean of y 
##     8.106     9.422

set.seed(1234)
dat <- data.frame(cond = factor(rep(c("A","B"), each=200)), 
                  rating = c(rnorm(200),rnorm(200, mean=.8)))

hist(subset(dat, dat$cond=="A")$rating)

hist(subset(dat, dat$cond=="B")$rating)

a = sample(subset(dat, dat$cond == "A")$rating, 3)
b = sample(subset(dat, dat$cond == "B")$rating, 3)

t.test(a, b)

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -0.216, df = 2.5562, p-value = 0.845
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.336276  4.719061
## sample estimates:
## mean of x mean of y 
## 0.4207416 0.7293491

a = sample(subset(dat, dat$cond == "A")$rating, 3)
b = sample(subset(dat, dat$cond == "B")$rating, 3)

t.test(a, b)

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -1.5844, df = 2.822, p-value = 0.2169
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.341895  1.173371
## sample estimates:
## mean of x mean of y 
##  1.122729  2.206991

a = sample(subset(dat, dat$cond == "A")$rating, 3)
b = sample(subset(dat, dat$cond == "B")$rating, 3)

t.test(a, b)

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -2.2632, df = 3.7663, p-value = 0.09051
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.9881229  0.3406797
## sample estimates:
##  mean of x  mean of y 
## -0.1738133  1.1499083

• tworzymy populację o średniej mi1 i odchyleniu standardowym si1

• tworzymy populację o średniej mi2 i odchyleniu standardowym si2

• z każdej populacji losujemy po kilka (NP) próbek

• sprawdzamy czy populacje są istotnie różne testem t

• procedurę losowania i testowania powtarzamy 1000 razy

pval <- function(NP=3, mi1=0, mi2=1, si1=1, si2=1)
{
  a = rnorm(NP, mi1, si1)
  b = rnorm(NP, mi2, si2)
  t = t.test(a, b)
  t$p.value
}

Wyniki 1

Wyniki 2

Wyniki 3

Wyniki 4

• Moc testu (moc statystyczna) to prawdopodobieństwo niepopełnienia błędu drugiego rodzaju – nieodrzucenia hipotezy zerowej, gdy w rzeczywistości jest ona fałszywa. Im większe jest to prawdopodobieństwo, tym lepszy jest dany test jako narzędzie do różnicowania między hipotezą prawdziwą i fałszywą.

• Moc zależy bezpośrednio i przede wszystkim od poniższych czynników:

  • wielkości próby użytej w badaniu,

  • rzeczywistej wielkości efektu na tle losowej zmienności w populacji,

  • przyjętego poziomu istotności α (najczęściej 0,05).

• pwr.t.test(d=1, n=3, sig.level = 0.05)

  • d - \(\ \frac{\mu_1 - \mu_2}{\sigma}\)

  • n - liczba powtórzeń

  • sig.level - prawdopodobieństwo, że średnie są równe (\(H_0\)!)

• pwr.anona.test

• pwr.norm.test

library(pwr)

pwr.t.test(d=1, n=3, sig.level = 0.05)

## 
##      Two-sample t test power calculation 
## 
##               n = 3
##               d = 1
##       sig.level = 0.05
##           power = 0.1587909
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

pwr.t.test(d=1, n=6, sig.level = 0.05)

## 
##      Two-sample t test power calculation 
## 
##               n = 6
##               d = 1
##       sig.level = 0.05
##           power = 0.3473537
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Moc testu t

Labirynt

Labirynt

\(Hipoteza: \mu_{0} = \mu_{1} = \mu_{2}\)

\(P_{całkowite(3)} = prawdziwe \ jest \ pierwsze \ i \ drugie \ porównanie.\)

\(P_{całkowite(3)} = P_{1} P_{2}\)

\(P_{całkowite(n)} = \prod_{k=0}^{n-1} P\)

Labirynt

Labirynt

Labirynt

Labirynt

Labirynt

Wykresy

Wykresy

a_thaliana <- read.csv("~/Dokumenty/repetytorium/a_thaliana.csv",
                       sep=";")
library(readxl)
a_thaliana2 <- read_excel("a_thaliana.xls")

library(ggplot2)

p = ggplot(data = a_thaliana)

p2 = p + geom_point(aes(x = as.factor(nutrient), y = total.fruits))

p3 = p2 + geom_boxplot(aes(x = as.factor(nutrient), y = total.fruits))

p4 = p + geom_jitter(aes(x = as.factor(nutrient), y = total.fruits))

p5 = p4 + xlab("skład pozywki") + ylab("liczba owocow")

p6 = p5 + theme_classic()

p4 = p + geom_jitter(aes(x = as.factor(nutrient),
                         y = total.fruits,
                         color = reg))
p5 = p4 + xlab("skład pozywki") + ylab("liczba owocow")
p6 = p5 + theme_classic()

p4 = p + geom_jitter(aes(x = as.factor(nutrient),
                         y = total.fruits,
                         color = reg,
                         shape = status))
p5 = p4 + xlab("skład pozywki") + ylab("liczba owocow")
p6 = p5 + theme_classic()

• Przykłady zastosowania

• Oficjalna dokumentacja

• R for data science

• R graphics cookbook