Критерій Кларка-Еванса в R

Критерій Кларка-Еванса дозволяє визначити як розподілені точки у просторі – випадково, рівномірно або груповим зразком. Раніше на 50º North вже був пост про цей критерій. У ньому докладно були описані теоретичні аспекти розрахунку критерію Клака-Еванса та його інтерпретування. Тому зараз ми не будемо вдруге розбирати їх. У цьому пості ще раз звернемо увагу до практичних аспектів. Раніше писали про роботу з критерієм Кларка-Еванса в ArcGIS. Цього разу познайомимося з розрахунком цього критерію в R.

R – це мова програмування та середовище для статистичних розрахунків. У ньому реалізовані всі стандартні процедури статистичного аналізу. А також є безліч додаткових пакетів з функціями, що дозволяють вирішувати різни спеціалізовані завдання, у тому числі в сфері просторової статистики.

Історія R починається в 1996 році в Новій Зеландії. Два вчених з університету Окленда – Росс Іхака та Роберт Джентльмен створили його як вільний аналог S-Plus. Головні переваги R це – безкоштовність та відкритий код. Зо роки свого існування R набув такої популярності, що зараз його можна назвати улюбленою мовою програмування професійних вчених-статистиків. Фактично на наших очах він перетворюється на професійний стандарт у цій галузі. Внаслідок цього R дуже швидко оновлюється. Всі найсучасніші досягнення статистики одразу з’являються в ньому у вигляді додаткових пакетів. Рахунок цих пакетів вже йде на тисячі. Тому не буде перебільшенням сказати, що R перевершує всі комерційні статистичні програми, що зараз існують.

Однією вадою R для багатьох користувачів є відносна складність навчання роботі з ним. Необхідно працювати в режимі командного рядка. Зараз існують віконно-кнопочні графічні інтерфейси користувача, наприклад RCommander, RKWard та інші. Але в них реалізовані стандартні процедури статистичного аналізу. Для роботи з критеріями просторової статистики необхідно вміти писати код на R. І ми хочемо показати, що це не так важко, як може спочатку здаватися.

 

Цей пост починає серію постів про просторову статистку в середовищі статистичних вирахувань R. Як і наступні пости з цієї серії, він буде орієнтований на тих, хто має базові навички роботи в R. Якщо Ви їх ще не маєте, то можете звернутися до підручників та онлайнових курсів з R. Добре, що час, коли вони всі були лише англомовні, вже давно пройшов. Основну увагу ми будемо звертати на функції, що необхідні для розрахунку критеріїв просторової статистики. Але по можливості ми також будемо коротко описувати також і базові функції. Щоб полегшити читання, у тексті посту назви пакетів, функцій, об’єктів та всі елементи коду будуть виокремлюватися шрифтом. Окремо в тексті посту будуть наведені фрагменти коду. У цих фрагментах елементи коду будуть позначені різними кольорами. Функції будуть позначені блакитним кольором, назви аргументів функцій – темно-зеленим, назви об’єктів – жовтогарячим, решта елементів залишаться  чорними.

Весь код на R, якій ми розберемо в цьому пості, наведений на рисунку 1. На цьому рисунку показаний інтерфейс R. Ліворуч в ньому відкритий редактор коду, а праворуч – консоль. Щоб запустити виконання коду, треба в редакторі коду виокремити фрагмент коду та натиснути кнопку, що на рисунку 1 обведена червоною рамкою. Після цього результат з’явиться в консолі. Також можна вводити код безпосередньо в консоль та запускати його натисненням кнопки Enter на клавіатурі.

Рис. 1. Інтерфейс R

 

Розрахунок критерію Кларка-Еванса в R розберемо на прикладі аналізу розміщення населенні пунктів у просторі. Дані для аналізу зображені на рисунку 2. Це міста, села, селища та хутора на території трьох районів, що розташовані у Бєлгородській області – Корочанський, Новооскольский та Чернянський райони.

Рис. 2. Населенні пункті на досліджуваній території

 

1. Імпорт шейп-файлів в R

Для тесту Кларка-Еванса в R як первинні дані можна використати векторні просторові дані у популярному форматі шейп-файлу. Треба мати два файли. Один з них – це шейп із точковою геометрією, для якого необхідно виконати тест Кларка-Еванса. Другий – це шейп із полігональною геометрією, що окреслює межі досліджуваної території. Він потрібний для оцінки площі цієї території. В нашому прикладі використовуємо файли під назвами settlement.shp (точки розташування населених пунктів) та rayon.shp (територія адміністративного району).

Для імпортування шейп-файлів до середовища R є різни засоби. Ми зробимо це за допомогою пакету maptools. Цей пакет містить функції для читання просторових даних та маніпуляцій з ними.

1) Для імпортування шейп-файлів спочатку встановлюємо робочу папку, з якої ми їх і будемо імпортувати. Щоб дізнатися, яка папка є робочою за замовчуванням, треба запустити функцію getwd. Жодних аргументів для неї вказувати не треба.

getwd()

Після виконання функції getwd у консоль буде виведений шлях до робочої папки. Якщо ця папка для Вас є зручною, то треба скопіювати в неї шейп-файли та переходити до їх імпортування в R. Будьте уважними під час копіювання, шейп-файл – це насправді декілька файлів з одним іменем та різними розширеннями (shp, shx, dbf та інші, детальніше читайте тут). Обов’язково треба копіювати всі ці файли разом, інакше може бути втрачена частина інформації або взагалі шейп-файл буде неможливо прочитати.

Якщо ви маєте бажання працювати з іншою папкою, зробіть її робочою за допомогою функції setwd. Вона має один аргумент – dir. Він вказує шлях до робочої папки (його треба брати у лапки).

setwd (dir = 'D:/SpatialData')

Зверніть увагу, на те яка риска стоїть у шляху до робочої папки. В R треба використати риску з правим нахилом.
2) Пакет maptools відноситься до додаткових пакетів R. Тобто він відсутній в стандартній комплектації R, його треба інсталювати додатково за допомогою функції install.packages, яка завантажує пакет з Інтернету та інсталює його в папку з програмою R. Аргумент pkgs – це назва пакету. Якщо завдати аргументу dependencies значення TRUE, програма інсталює не тільки пакет maptools, але й інші пакети, що потрібні для його функціонування.

install.packages (pkgs = 'maptools', dependencies = TRUE)

Якщо пакет вже є на комп’ютері, треба його запустити. Деякі пакети при запуску R запускаються автоматично. Але всі додаткові пакети треба запускати вручну. Для цього використовують функцію library. Її аргумент package – це ім’я пакету (брати його в лапки не обов’язково).

library (package = maptools)

3) Після запуску пакету maptools вже можна імпортувати в R шейп-файли. Для цього в maptools призначені функції readShapePoints та readShapePoly. Перша функція імпортує шейпи з точковою геометрією, а друга – з полігональною геометрією. Ім’я файлу, що ми імпортуємо, для обох функцій задається аргументом fn.

Кожний шейп-файл ми записуємо в окремий об’єкт (шар з населеними пунктами в об’єкт Setl, а шар з полігоном досліджуваної території – в об’єкт Rayon). Для цього спочатку ми пишемо ім’я об’єкту, далі знак присвоєння ( <- ), а потім вже саму функцію.

Setl <- readShapePoints (fn = 'settlement.shp')

Rayon <- readShapePoly (fn = 'rayon.shp')

Пакет maptools працює зі специфічними класами об’єктів. Для точкових просторових даних це клас SpatialPointDataFrame. Для полігональних даних це клас SpatialPoligonDataFrame. Але пакет spatstat не працює з цими класами об’єктів. Тому треба дані конвертувати в об’єкти інших класів. Для точкових об’єктів в spatstat є клас ppp. А межі досліджуваної території вказуються об’єктом класу owin.

 

2. Створення об’єкту класу owin

Для аналізу точкових просторових даних в R існує пакет spatstat. Це теж додатковий пакет, як і maptools. І його також необхідно запускати вручну за допомогою функції library.

library (package = spatstat)

Спочатку треба створити об’єкт класу owin, а після нього об’єкт класу ppp. Бо owin використовується в процесі створення ppp. Для нашого прикладу зручно перетворити на owin об’єкт Rayon, який ми зробили з шейп-файлу rayon.shp. Для цього є функція as.owin. Її аргумент W – це об’єкт, з якого ми робимо об’єкт класу owin.

RayonOwin <- as.owin (W = Rayon)

 

3. Створення об’єкту класу ppp

Пакет spatstat працює зі спеціальним класом об’єктів – ppp. Об’єкти цього класу створюються за допомогою функції, яка теж має назву ppp. Для нашого прикладу ми задаємо значення для трьох аргументів цієї функції. Аргумент x та y – це два вектори з координатами точок (відповідно X та Y координати у прямокутній системі координат). Аргумент window – це межа території, яку ми досліджуємо. Щоб завдати координати точок ми будемо використати об’єкт Setl, який ми створили з шейп-файлу settlement.shp. Цей об’єкт класу ppp вміщує інформацію про координати точок. Але треба знати як її з нього прочитати.

Об’єкт класу ppp – це складний об’єкт. Він складається з окремих частин (слотів). Щоб прочитати інформацію зі слоту, треба без пробілів написати та запустити код – ім’я об’єкту, знак @, ім’я слоту. Слот, що містить інформацію про координати в будь якому об’єкті класу ppp має назву coords.

Перший рядок коду, що наведений нижче, виводить в консоль разом координати X та Y у вигляді двох стовпчиків. Результат його виконання наведений на рисунку 3 ліворуч. Щоб звернутися до окремого стовпчика, треба дописати до цього коду (без пробілів) у квадратних дужках номер стовпчика. Приклад цього – другий рядок коду, що наведений нижче.

Стовпчик із номером 1 – це координата X, стовпчик з номером 2 – це координата Y. Зверніть увагу на те, що перед номером стовпчика стоїть кома. Вона відділяє номер рядку від номера стовпчика. У нашому випадку ми звертаємося не до окремого рядка, а до всіх рядків з обраного стовпчика. Тому номер рядка писати не треба, але кому все одно треба обов’язково ставити.

Setl@coords

Setl@coords[,1]

Приклад виведення координат у консоль R показаний на рисунку 3.

Рис. 3. Координати, що виведені у консоль R (ліворуч – обидві координати, X та Y, праворуч зверху– координата X, праворуч внизу – координата Y)

 

Як значення для аргументу window можна вказати об’єкт класу owin.В нашому випадку це об’єкт RayonOwin.

Setl_ppp <- ppp (x = Setl@coords[,1],

                y = Setl@coords[,2],

                window = RayonOwin)

Після того, як ми створили об’єкт класу ppp, його можна зобразити графічно (рис. 4). Але робити картографічні зображення краще в ГІС. А в R використати візуалізацію лише для перевірки правильності створення об’єкту.

plot (x = Setl_ppp)

Рис. 4. Результат візуалізації об’єкту класу ppp

 

Якщо все зроблено правильно, ми побачимо всі точки с таким просторовим розташуванням, яке було в шейп-файлі. А навколо точок буде проведена межа досліджуваної території. Порівняйте рисунок 2 та рисунок 4, і Ви бачите, що в нашому прикладі все на своєму місці.

 

4. Виконання тесту Кларка-Еванса

Коли ми маємо об’єкт класу ppp, ми можемо вже нарешті виконати тест Кларка-Еванса. Робиться він за допомогою функції clarkevans.test. Її аргумент X – це назва об’єкту класу ppp, для якого ми виконуємо тест.

clarkevans.test (X = Setl_ppp)

Після виконання функції clarkevans.test у консоль будуть виведеними значення критерію Кларка-Еванса та імовірність помилки першого роду (p-значення). Для нашого прикладу критерій Кларка-Еванса (R) приблизно рівний 1,11, що статистично значно більше ніж 1,00 (p<0,05). Ці результати свідчать, що населенні пункти на території Корочанського, Новоскольского та Чернянського району розподілені у просторі рівномірно (рис. 5).

 

Рис. 5. Результат тесту Кларка-Еванса в консолі R (зверху) та результат тесту, що отриманий за допомогою ArcGIS (внизу)

 

Якщо порівняти результати, отримані за допомогою R та ArcGIS, то бачимо, що вони однакові. Але в R є опції розрахунку критерію Кларка-Еванса, які відсутні в ArcGIS. Але про це вже буде окремий пост.