Анализ статистических зависимостей
Alexander Novopoltsev
Sunday, September 14, 2014
Корреляционный анализ
# загрузим данные об авто
cars=read.csv("cars.csv")
#присоединим к списку текущих переменных
attach(cars)Коэффициенты корреляции
Анализ количественных переменных
#коэффициент парной корреляции: линейная зависимость
cor(5:15,7:17)## [1] 1#нарушим зависимость
cor(5:15,c(7:16,50))## [1] 0.6934#данные cars: объем двигателя и мощность
#опция "use" задает метод исключения пропущенных значений
cor(engine.size,horsepower,use="complete.obs")## [1] 0.8108#выберем несколько количественных переменнных
cars.subset=cars[c("horsepower","peak.rpm","price")]
#выведем некоторые строки (в т.ч.c пропущенными значениями)
cars.subset[c(9:11,43:47,127:133),]## horsepower peak.rpm price
## 9 140 5500 23875
## 10 160 5500 NA
## 11 101 5800 16430
## 43 100 5500 10345
## 44 78 4800 6785
## 45 70 5400 NA
## 46 70 5400 NA
## 47 90 5000 11048
## 127 207 5900 32528
## 128 207 5900 34028
## 129 207 5900 37028
## 130 288 5750 NA
## 131 NA NA 9295
## 132 NA NA 9895
## 133 110 5250 11850#корреляционная матрица для нескольких переменных
#по умолчанию use="everething" : для расчета используются все наблюдения
cor(cars.subset)## horsepower peak.rpm price
## horsepower 1 NA NA
## peak.rpm NA 1 NA
## price NA NA 1#use="complete.obs" : только наблюдения, для которых нет проп.зн. по всем переменным
cor(cars.subset,use="complete.obs")## horsepower peak.rpm price
## horsepower 1.0000 0.1079 0.8105
## peak.rpm 0.1079 1.0000 -0.1016
## price 0.8105 -0.1016 1.0000#use="pairwise.complete.obs": для вычисления каждого коэффициента исп. набл., у которых нет проп. зн. только для пары соотв. перем.
cor(cars.subset,use="pairwise.complete.obs")## horsepower peak.rpm price
## horsepower 1.0000 0.1310 0.8105
## peak.rpm 0.1310 1.0000 -0.1016
## price 0.8105 -0.1016 1.0000Анализ порядковых переменных
Используются коэффициенты корреляции Кендалла и Спирмена
#рейтинги кредитоспособности заемциков, вычисленные с помощью разных методик
#меньше число - надежнее заемщик
class1=c(2,3,1,4,2,2,1,2,3,4)
class2=c(2,2,1,4,2,2,2,2,2,4)
#коэффициент корреляции Кендалла
cor(class1,class2,method = "kendall")## [1] 0.7646#для случайных переменных на основе биномиального распределения
bin1=rbinom(n = 100,size = 3,prob = 0.5)+1
bin2=rbinom(n = 100,size = 3,prob = 0.5)+1
bin=cbind(bin1,bin2)
head(bin)## bin1 bin2
## [1,] 2 3
## [2,] 2 3
## [3,] 3 3
## [4,] 3 3
## [5,] 4 2
## [6,] 4 3cor(bin,use="pairwise",method="kendall") # очень малая корреляция## bin1 bin2
## bin1 1.00000 0.02088
## bin2 0.02088 1.00000#коэффициент корреляции Спирмена
cor(bin1,bin2,method="spearman")## [1] 0.02267Статистические тесты
#проверим статистическую значимость на уровне 0.05
#нулевая гипотеза: "коэфф. корр. равняется нулю"
cor.test(engine.size,horsepower) #корреляция не равна нулю##
## Pearson's product-moment correlation
##
## data: engine.size and horsepower
## t = 19.64, df = 201, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7577 0.8532
## sample estimates:
## cor
## 0.8108Другие способы представления корреляции
#линия регрессии
plot(horsepower~engine.size) #построим диаграмму рассеяния
#нанесем линию регресси: коэффициент корреляции - наклон линии регрессии
abline(lm(horsepower~engine.size))
#матрица корреляции
(C=cor(cars.subset,use="complete.obs"))## horsepower peak.rpm price
## horsepower 1.0000 0.1079 0.8105
## peak.rpm 0.1079 1.0000 -0.1016
## price 0.8105 -0.1016 1.0000#представление в виде треугольной таблицы символов
#отмечаются только значимые корреляции
symnum(C)## h p. pr
## horsepower 1
## peak.rpm 1
## price + 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1#представление в виде эллипсов: связь с двумерным норм. распр.
#чем выше корреляция - тем более вытянут эллипс
library(ellipse)
plotcorr(C)
Таблицы сопряженности
#частоты моделей для каждой марки
table(make)## make
## alfa-romero audi bmw chevrolet dodge
## 3 7 8 3 9
## honda isuzu jaguar mazda mercedes-benz
## 13 4 3 17 8
## mercury mitsubishi nissan peugot plymouth
## 1 13 18 11 7
## porsche renault saab subaru toyota
## 5 2 6 12 32
## volkswagen volvo
## 12 11#частоты по цене и типу кузова
price.cut = cut(price,c(0,12000,20000,max(price,na.rm=T))) #разобьем цены по категориям
levels(price.cut)=c("cheap","okay","expensive") #назовем каждую категорию
(t=table(price.cut,body.style))## body.style
## price.cut convertible hardtop hatchback sedan wagon
## cheap 1 4 51 49 13
## okay 3 0 16 28 11
## expensive 2 4 1 17 1#в относительных частотах
options(digits=2)
prop.table(t,1) # 1: по строкам## body.style
## price.cut convertible hardtop hatchback sedan wagon
## cheap 0.0085 0.0339 0.4322 0.4153 0.1102
## okay 0.0517 0.0000 0.2759 0.4828 0.1897
## expensive 0.0800 0.1600 0.0400 0.6800 0.0400#по цене, типу кузова, а также по расходу топлива на трассе (MPG.highway)
mpg = cut(highway.mpg,c(0,20,30,max(highway.mpg))) # преобразуем в категоральную переменную
table(price.cut,body.style,mpg)## , , mpg = (0,20]
##
## body.style
## price.cut convertible hardtop hatchback sedan wagon
## cheap 0 0 0 0 0
## okay 0 0 0 0 0
## expensive 1 1 0 7 0
##
## , , mpg = (20,30]
##
## body.style
## price.cut convertible hardtop hatchback sedan wagon
## cheap 1 3 12 8 4
## okay 3 0 16 23 10
## expensive 1 3 1 10 1
##
## , , mpg = (30,54]
##
## body.style
## price.cut convertible hardtop hatchback sedan wagon
## cheap 0 1 39 41 9
## okay 0 0 0 5 1
## expensive 0 0 0 0 0#можно создать плоскую (двумерную таблицу)
ftable(price.cut,body.style,mpg)## mpg (0,20] (20,30] (30,54]
## price.cut body.style
## cheap convertible 0 1 0
## hardtop 0 3 1
## hatchback 0 12 39
## sedan 0 8 41
## wagon 0 4 9
## okay convertible 0 3 0
## hardtop 0 0 0
## hatchback 0 16 0
## sedan 0 23 5
## wagon 0 10 1
## expensive convertible 1 1 0
## hardtop 1 3 0
## hatchback 0 1 0
## sedan 7 10 0
## wagon 0 1 0#Titanic: данные о пассажирах Титаника
ftable(Titanic, row.vars = 1:3)## Survived No Yes
## Class Sex Age
## 1st Male Child 0 5
## Adult 118 57
## Female Child 0 1
## Adult 4 140
## 2nd Male Child 0 11
## Adult 154 14
## Female Child 0 13
## Adult 13 80
## 3rd Male Child 35 13
## Adult 387 75
## Female Child 17 14
## Adult 89 76
## Crew Male Child 0 0
## Adult 670 192
## Female Child 0 0
## Adult 3 20Регрессионный анализ
Линейная регрессия
Данные
florida:
количество голосов у 13 кандидатов в 67 огругах штата Флорида на выборах в СШАflorida=read.csv("florida.csv")С одним фактором
#x и y - количество голосов для двух кандидатов BUSH и BUCHANAN соответственно во всех огругах
x=florida$BUSH
y=florida$BUCHANAN
#оценить модель регрессии
mod=lm(y~x) # "у зависит от x"
#построить графики: диаграмма рассеяния
plot(x,y)
#нанести линию регрессии
abline(mod)
Свойства оцененной модели
#статистика по оцененной модели
#коэффициент при х значим (Pr близко к 0)
#качество модели плохое: статистика R-квадрат мала (R-squared=0.389)
summary(mod)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -907.5 -46.1 -29.2 12.3 2610.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.53e+01 5.45e+01 0.83 0.41
## x 4.92e-03 7.64e-04 6.43 1.7e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 354 on 65 degrees of freedom
## Multiple R-squared: 0.389, Adjusted R-squared: 0.38
## F-statistic: 41.4 on 1 and 65 DF, p-value: 1.73e-08#оцененные коэффициенты
coefficients(mod)## (Intercept) x
## 45.2899 0.0049#остатки
res=residuals(mod)
plot(res)
#построим все графики остатков
par(mfrow=c(2,2))
plot(mod)
par(mfrow=c(1,1))удаление аномальных наблюдений
plot(x,y)
abline(mod)
На графике видны наблюдения, лежащие далеко в стороне от кривой.
Функция identify(x,y,n=2) позволит определить индексы 2 точек, которые выбираются вручную с помощью клика мышкой по графику после вызова данной функции.
это будут наблюдения с номерами 13,50#удалим указанные точки из выборки и переоценим модель
x=x[-c(13,50)];y=y[-c(13,50)]
mod2=lm(y~x)
#построим график
plot(x,y)
abline(mod2)
#выведем статистику по новой модели
#качество модели удовлетворительное (R-square=0.865),
#все коэффициенты значимы на уровне 0.05
summary(mod2)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -200.9 -28.5 -11.1 27.5 281.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.85e+01 1.31e+01 2.93 0.0047 **
## x 4.40e-03 2.19e-04 20.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 82 on 63 degrees of freedom
## Multiple R-squared: 0.865, Adjusted R-squared: 0.863
## F-statistic: 403 on 1 and 63 DF, p-value: <2e-16Множественная линейная регрессия
#несколько факторов: объема двигателя (engine.size) и макс. кол. оборотов (peak.rpm)
mod.mul=lm(horsepower~engine.size+peak.rpm,data=cars)
summary(mod.mul)##
## Call:
## lm(formula = horsepower ~ engine.size + peak.rpm, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.22 -10.87 -2.34 7.31 100.83
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.52e+02 1.64e+01 -9.32 <2e-16 ***
## engine.size 8.51e-01 3.30e-02 25.78 <2e-16 ***
## peak.rpm 2.90e-02 2.88e-03 10.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19 on 200 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.773, Adjusted R-squared: 0.77
## F-statistic: 340 on 2 and 200 DF, p-value: <2e-16Нелинейная регрессия
n=24
#добавим возмущения в виде БСВ
x = runif(n)
y = x^3 + runif(n, min = -0.1, max = 0.1)
plot(x, y)
s <- seq(from = 0, to = 1, length = 50)
lines(s, s^3, lty = 2)
df <- data.frame(x, y)
#оценим модель
#trace показывает значения power на каждой итерации
m <- nls(y ~ I(x^power), data = df, start = list(power = 1), trace = T)## 1.7 : 1
## 0.32 : 1.9
## 0.091 : 2.6
## 0.082 : 2.9
## 0.082 : 2.9
## 0.082 : 2.9lines(s, predict(m, list(x = s)), col = "green")
summary(m)##
## Formula: y ~ I(x^power)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## power 2.887 0.159 18.1 4.2e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06 on 23 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 8.48e-06