Аналіз відкритих даних в R, частина 1

Введення

На момент написання статті більшість додатків на основі відкритих даних (на офіційних сайтах data.mos.ru/apps і data.gov.ru) представляють собою інтерактивні довідники по інфраструктурі міста чи поселення з наочною візуалізацією і часто з опцією вибору оптимального маршруту. Мета цієї і наступних публікацій полягає в тому, щоб привернути увагу спільноти до обговорення стратегій аналізу відкритих даних, у т. ч. спрямованих на прогнозування, побудова статистичних моделей та вилучення інформації, яка не представлена у явному вигляді. В якості інструментарію використовується мова R і середовище розробки RStudio.

В якості ілюстрації простоти використання R розглянемо невеликий приклад і пропонуємо на карті реєстру залучення КК, ТСЖ за 2013 рік в Ульяновську. У цьому наборі даних міститься інформація про штрафи різних керуючих компаній. Для візуалізації вибрано кількість усних зауважень, кількість штрафів і штрафи в рублях за все представлене в датасете час. Цих даних(як і всім іншим з відкритих джерел, що мені попадалися) властиво мати багато пропусків, помилок, і т. д. тому в R завантажуються предобработанные дані:


ulyanovsk.data <- read.csv('./data/formated_tsj.csv')


Для візуалізації використовується інтерфейс google maps і бібліотека ggplot2, для організації графіків на пристрої виводу використовується gridExtra:

Кодlibrary(ggmap)
library(ggplot2)
library(gridExtra)

У датасете є стовпці lat і lon, в яких міститься інформація про широті і довготі. Нам потрібна карта не всього Ульяновська а лише його частини, координати якої знаходяться у вибірці:

Кодbox < — make_bbox(lon, lat, data = ulyanovsk.data)
ulyanovsk.map < — get_map(location = 'ulyanovsk', zoom=calc_zoom(box))
p0 < — ggmap(ulyanovsk.map)
p0 < — p0 + ggtitle(«Частина карти Ульяновська») + theme(legend.position=«none»)

Наносимо на карту штрафи в рублях(солбец penalties_rub_total):

Кодp1 < — ggmap(ulyanovsk.map)
p1 < — p1 + geom_point(data = ulyanovsk.data, aes(x = lon, y = lat, color = penalties_rub_total, size = penalties_rub_total))
p1 < — p1 + scale_color_gradient(low = 'yellow', high = 'red')
p1 < — p1 + ggtitle(«Штрафи, в рублях») + theme(legend.position=«none»)

Кількість штрафів у штуках:

Кодp2 < — ggmap(ulyanovsk.map)
p2 < — p2 + geom_point(data = ulyanovsk.data, aes(x = lon, y = lat, color = penalties_total, size = penalties_total))
p2 < — p2 + scale_color_gradient(low = 'yellow', high = 'red')
p2 < — p2 + ggtitle(«Кількість штрафів») + theme(legend.position=«none»)

Кількість усних зауважень у штуках:

Кодp3 < — ggmap(ulyanovsk.map)
p3 < — p3 + geom_point(data = ulyanovsk.data, aes(x = lon, y = lat, color = oral_comments_total, size = oral_comments_total))
p3 < — p3 + scale_color_gradient(low = 'yellow', high = 'red')
p3 < — p3 + ggtitle(«Кількість усних зауважень») + theme(legend.position=«none»)

Будуємо разом:


grid.arrange(p0, p1, p2, p3, ncol=2)




Розмір точки і величина червоного монотонно зростають із параметром(наприклад, чим більше штрафів, тим точка більше і червонішими). З графіків видно, що в досліджуваному датасете керуючі компанії з правого берега Ульяновська мають більшу кількість штрафів та зауважень, ніж з лівого.

Дослідження статистики ДТП

В цьому розділі використані дані єдиної міжвідомчої інформаційно-статистичної системи. У датасете міститься показник смертності від ДТП на 100 тис. осіб по регіонах РФ за 2011-2014 рік:

Кодaccidents < — read.csv('./data/accidents.csv')
head(accidents)




Середнє для кожного року по областях:


apply(accidents[,2:5],2,FUN = mean)


## y2011 y2012 y2013 y2014
## 21.36942 21.55694 20.18582 20.34727

Гістограми:

Кодlibrary(reshape2)
accidents.melt < — melt(accidents[,2:5])
ggplot(accidents.melt, aes(x = value, fill = variable)) + geom_histogram(alpha = 0.3)




Функція melt реорганізує датасет таким чином, що стовпці value зберігається значення показника, у стовпці variable рік, за який показник був отриманий:


head(accidents.melt)




Цікаво порівняти, чи є статистично значущим відмінність показників для різних років. Приймаються наступні модельні припущення:

1) Значення показника є випадковою величиною, що підкоряється розподілу Стьюдента
2) Для кожного року є свої параметри розподілу, nu — параметр нормальності, sigma — дисперсія, mu — середнє
3) Параметри nu, sigma, mu так само є випадковими величинами, зі своїми розподілами, не залежними від року

Формальна постановка завдання: використовуючи дані встановити, чи є параметр розподілу mu для різних років різним. Якщо виписати для нашого випадку формулу Байєса(яку я не буду тут наводити, щоб не захаращувати статтю), то виходить ієрархічна модель:

1) показник смертності (поле value в accidents.melt) ~ t(mu,sigma,nu)
2) mu | group_id ~ N(mu_hyp, sigma_hyp )
3) sigma | group_id ~ Unif(a_hyp, b_hyp | group_id)
4) nu | group_id ~ Exp(lambda_hyp)

Де mu | group_id ~ N(mu_hyp, sigma_hyp ) означає, що параметр mu для групи group_id має нормальний розподіл з гиперпараметрами mu_hyp, sigma_hyp.
Змінна group_id приймає значення в залежності від року, 1 для 2011, 2 для 2012 і т. д.:


accidents.melt$group_id <- as.numeric(accidents.melt$variable)


Для моделювання R використовується інтерфейс до мови Jags і допоміжна функція plotPost з файлу DBDA2E-utilities.R, який можна знайти на тут(автор — John Kruschke):

Кодlibrary(rjags)
source('DBDA2E-utilities.R')


Опис моделі відбувається на мові jags, тому спочатку створюється файл або символьна змінна всередині R за моделлю:

КодmodelString < — '

model {

for(i in 1:N)
{

y[i] ~ dt(mu[group_id[i]], 1/sigma[group_id[i]]^2, nu[group_id[i]])

}

for(j in 1:NumberOfGroups)
{

mu[j] ~ dnorm(PriorMean, 0.01)
sigma[j] ~ dunif(1E-1, 1E+1)
nuTemp[j] ~ dexp(1/19)
nu[j] < — nuTemp[j] + 1

}

}

'


Після того, як модель специфікована, вона передається jags разом з даними:

КодmodelJags < — jags.model(file = textConnection(modelString),
data = list(y = accidents.melt$value, group_id = accidents.melt$group_id, N = length(accidents.melt), NumberOfGroups = 4, PriorMean = (20)),
n.chains = 8, n.adapt = 1000)


Jags використовує MCMC і реалізує модель випадкових блукань у просторі наших параметрів mu, sigma, nu. Причому ймовірність того, що параметри візьмуть конкретне значення пропорційна їх спільної ймовірності, хоча спільна ймовірність не враховується у явному вигляді. Таким чином jags повертає список значень mu, sigma, nu, тому можна оцінити їх ймовірності використовуючи гістограми. Наступні рядки запускають нашу модель і повертають результат:

Кодupdate(modelJags,2000)
summary(samplesCoda < — coda.samples(model = modelJags, variable.names = c('nu','mu','sigma'), n.iter = 50000, thin = 500))


Для відповіді на запитання про значення параметра mu для різних років розглянемо гістограми різниці значень параметрів:

КодsamplesCodaMatrix < — as.matrix(samplesCoda)
samplesCodaDataFrame < — as.data.frame(samplesCodaMatrix)
samplesCodaDataFrame < — samplesCodaDataFrame[,1:4]
plotPost(samplesCodaMatrix[,1]- samplesCodaMatrix[,2], cenTend = 'mean', main = '2011 vs. 2012')
plotPost(samplesCodaMatrix[,1]- samplesCodaMatrix[,3], cenTend = 'mean', main = '2011 vs. 2013')
plotPost(samplesCodaMatrix[,1]- samplesCodaMatrix[,4], cenTend = 'mean', main = '2011 vs. 2014')
plotPost(samplesCodaMatrix[,2]- samplesCodaMatrix[,3], cenTend = 'mean', main = '2012 vs. 2013')
plotPost(samplesCodaMatrix[,2]- samplesCodaMatrix[,4], cenTend = 'mean', main = '2012 vs. 2014')
plotPost(samplesCodaMatrix[,3]- samplesCodaMatrix[,4], cenTend = 'mean', main = '2013 vs. 2014')









На графіках позначений HDI інтервал, значення в якому мають велику ймовірність, ніж значення поза інтервалу, а сума ймовірностей цих значень дорівнює 0,95. Видно, що значення 0 потрапляє в HDI, що дає підставу припустити, що параметр mu один і той же для всіх представлених в датасете років. Тут є спірний момент — значення HDI підбирається згідно знання предметної області, тобто можливе для цієї задачі значення 95% не саме підходяще.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.