Занятие 17

Статистика в R

Считаем данные по успеваемости, с которыми мы уже имели дело.

dat <- read.csv("http://www.ats.ucla.edu/stat/data/hsb2.csv")

Нарисуем гистограмму для оценки по математике

hist(dat$math)

Знак $ служит разделителем имени датафрейма и названия переменной — как точка в Python. Чтобы каждый раз не указывать название датафрейма, можно сделать attach.

attach(dat)
print(math)

##   [1] 41 53 54 47 57 51 42 45 54 52 51 51 71 57 50 43 51 60 62 57 35 75 45
##  [24] 57 45 46 66 57 49 49 57 64 63 57 50 58 75 68 44 40 41 62 57 43 48 63
##  [47] 39 70 63 59 61 38 61 49 73 44 42 39 55 52 45 61 39 41 50 40 60 47 59
##  [70] 49 46 58 71 58 46 43 54 56 46 54 57 54 71 48 40 64 51 39 40 61 66 49
##  [93] 65 52 46 61 72 71 40 69 64 56 49 54 53 66 67 40 46 69 40 41 57 58 57
## [116] 37 55 62 64 40 50 46 53 52 45 56 45 54 56 41 54 72 56 47 49 60 54 55
## [139] 33 49 43 50 52 48 58 43 41 43 46 44 43 61 40 49 56 61 50 51 42 67 53
## [162] 50 51 72 48 40 53 39 63 51 45 39 42 62 44 65 63 54 45 60 49 48 57 55
## [185] 66 64 55 42 56 53 41 42 53 42 60 52 38 57 58 65

hist(math)

Если вы закончили работу с датафреймом, можно сделать detach.

detach(dat)

Как связана успеваемость по математике с успеваемостью по чтению?

plot(math ~ read, data=dat)

Здесь мы сталкиваемся с формулой с тильдой, которая в Python использовалась для регрессий. Здесь она имеет более широкое применение и означает, в основном, сравнить то, что написано слева от тильды, с тем, что написано справа. В частности, plot делает scatter plot двух переменных.

Попробуем теперь понять, насколько успеваемость по разным предметам зависит от пола.

plot(math ~ female, data=dat)

Не очень понятная картинка, конечно. Это потому, что пол на самом деле — категориальная переменная (фактор в терминах R), и рисовать на неё scatter plot довольно странно. А вот если R объяснить, что это фактор, то всё будет хорошо.

dat$female=as.factor(dat$female)
dat$race=as.factor(dat$race)
plot(math ~ female, data=dat)

plot(math ~ race, data=dat)

Это уже знакомые нам ящики с усами.

А вот как можно группировать и аггрегировать данные.

aggregate(read ~ race, data=dat, FUN=mean)

##   race     read
## 1    1 46.66667
## 2    2 51.90909
## 3    3 46.80000
## 4    4 53.92414

aggregate(cbind(read,write,math,science) ~ race, data=dat, FUN=mean)

##   race     read    write     math  science
## 1    1 46.66667 46.45833 47.41667 45.37500
## 2    2 51.90909 58.00000 57.27273 51.45455
## 3    3 46.80000 48.20000 46.75000 42.80000
## 4    4 53.92414 54.05517 53.97241 54.20000

aggregate(read ~ female + race, data=dat, FUN=mean)

##   female race     read
## 1      0    1 47.30769
## 2      1    1 45.90909
## 3      0    2 52.33333
## 4      1    2 51.75000
## 5      0    3 46.85714
## 6      1    3 46.76923
## 7      0    4 54.51471
## 8      1    4 53.40260

aggregate(log(math) ~ female + race, data=dat, FUN=mean)

##   female race log(math)
## 1      0    1  3.887783
## 2      1    1  3.801973
## 3      0    2  4.064116
## 4      1    2  4.020988
## 5      0    3  3.808715
## 6      1    3  3.851613
## 7      0    4  3.975188
## 8      1    4  3.971309

Проверим с помощью t-теста гипотезу о том, как пол связан с успехами по математике.

t.test(math ~ female, data=dat)

## 
##  Welch Two Sample t-test
## 
## data:  math by female
## t = 0.411, df = 187.575, p-value = 0.6816
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.092206  3.193325
## sample estimates:
## mean in group 0 mean in group 1 
##        52.94505        52.39450

p-value в районе 0.6, то есть никак не связан. А вот с письмом другая история.

plot(write ~female, data=dat)

t.test(write ~ female, data=dat)

## 
##  Welch Two Sample t-test
## 
## data:  write by female
## t = -3.6564, df = 169.707, p-value = 0.0003409
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.499159 -2.240734
## sample estimates:
## mean in group 0 mean in group 1 
##        50.12088        54.99083

Ещё немножко регрессий.

fit <- lm(write ~ female + race, data=dat)
summary(fit)

## 
## Call:
## lm(formula = write ~ female + race, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.2665  -5.1350   0.7335   6.5054  20.7029 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  44.2971     1.8729  23.651  < 2e-16 ***
## female1       4.7153     1.2505   3.771 0.000216 ***
## race2        10.2735     3.1983   3.212 0.001541 ** 
## race3         0.8379     2.6556   0.316 0.752703    
## race4         7.2540     1.9272   3.764 0.000221 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.735 on 195 degrees of freedom
## Multiple R-squared:  0.1678, Adjusted R-squared:  0.1507 
## F-statistic: 9.826 on 4 and 195 DF,  p-value: 2.912e-07

fit <- lm (write ~ read, data=dat)
plot(write ~ read, data=dat)
abline(fit)

fit <- lm (read+write+math ~ as.factor(prog), data=dat)
summary(fit)

## 
## Call:
## lm(formula = read + write + math ~ as.factor(prog), data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -59.152 -17.380  -0.111  16.120  61.620 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       151.111      3.256  46.410  < 2e-16 ***
## as.factor(prog)2   18.041      3.892   4.636 6.46e-06 ***
## as.factor(prog)3  -11.731      4.488  -2.614  0.00964 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.84 on 197 degrees of freedom
## Multiple R-squared:  0.257,  Adjusted R-squared:  0.2494 
## F-statistic: 34.07 on 2 and 197 DF,  p-value: 1.964e-13

Занятие 17

Илья Щуров, НИУ ВШЭ

23 Mar 2015

Статистика в R

Программирование в R