Читать онлайн Проведение статистического анализа на языке программирования R в медико-биологических исследованиях. Учебное пособие бесплатно

При проведении статистического анализа аналитик сталкивается с тремя основными типами данных:

1) количественными;

2) качественными;

3) смешанными (сочетающими количественные и качественные характеристики).

Количественные данные38 – численные данные, имеющие шкалу измерения. Они могут быть разделены на дискретные и непрерывные.

Дискретные количественные данные – это числовые данные, полученные путем подсчета какой-либо величины.

Примером таких данных являются: вес пациента, измеряющийся в килограммах; рост пациента, измеряемый в сантиметрах; количество пациентов; концентрация различных веществ в крови или моче пациента и т. д. (результаты антропометрических, лабораторных и функциональных исследований пациентов).

Непрерывные количественные данные – изменение какой-либо количественной величины, измеренной за определенный период времени.

Примерами непрерывных количественных данных являются электроэнцефалограмма головного мозга, электрокардиограмма сердечной мышцы и др. В таблице 1 представлен пример дискретных количественных данных.

Качественные данные – это данные, описывающие признак предмета исследования естественным языком. Они могут быть номинальными и порядковыми.

Номинальные данные – это подгруппа качественных данных, используемая для именования переменных, не имеющих числового значения.

Примером качественного номинального признака предмета исследования являются пол пациента, наличие или отсутствие вредных привычек, занятие пациентом физической культурой и т. д.

Порядковые данные – это подгруппа качественных данных, имеющих порядок или масштаб. Часто качественные порядковые данные могут быть представлены числами, которым соответствует качественное описание. Примером порядковых качественных данных могут являться порядковый номер пациента, порядковый номер ответа на вопрос и т. д.

Наиболее распространенный пример качественных порядковых данных – ответ на вопрос по пяти-, десяти- и т. д. балльной шкале. Пример пятибалльной шкалы уверенности в каком-либо вопросе, применяющейся при проведении опросов:

1. Да.

2. Скорее да.

3. Затрудняюсь ответить.

4. Скорее нет.

5. Нет.

Каждому из пяти представленных ответов может соответствовать число от 1 до 5 в случае применения порядковой шкалы или от 0 до 1 в случае вероятностной шкалы. Выбор шкалы соответствия ограничивается только фантазией исследователя, и единственное требование, которое возникает при планировании сбора качественных данных – единообразие используемых шкал.

В языке программирования R качественное описание объекта, содержащееся в наборе данных, представляется как тип данных «фактор»39, а значения, принимаемые переменной, имеющей тип «фактор», называются уровнями фактора. В таблице 2 представлен пример факторных данных.

В данном случае идентификатор представлен численно буквенным кодом, идентифицирующим пациента в некоторой базе данных, а пол пациента обозначен буквами «М» (мужской) и «Ж» (женский).

В практике анализа и проведения эксперимента редко встречаются наборы данных, содержащие только количественные или только качественные данные. Чаще всего наборы данных имеют смешанный характер. В таблице 3 представлен пример смешанного набора данных.

Фактически набор данных представляет собой совокупность переменных (название каждого столбца), каждая из которых описывает качественную или количественную характеристику исследуемого объекта или явления. Соответственно, переменные, содержащие только качественные значения, называются факторами. Создание наборов медицинских данных детально рассматривается в курсе «Создание наборов данных»40.

Для каждого типа переменных или их совокупности существует свой набор статистических тестов, позволяющих проводить детальный анализ наличия или отсутствия различий между всевозможными уровнями факторов, наличие связи или отсутствие связей между данными, принадлежащими к разным уровням фактора, которые будут рассмотрены далее. Алгоритм анализа, представленный на рисунке 2, является верным для количественных и смешанных наборов данных. Методы анализа качественных данных не рассматриваются в настоящем учебном пособии.

Ниже представлены примеры наборов данных, содержащих смешанные данные, но с превалирующим количеством качественных переменных и с превалирующим числом количественных данных.

1.1. Пример количественных и качественных данных

В практике анализа наборы данных, содержащие только качественные переменные, возникают при сборе методом опроса или анкетирования (даже в этом случае результаты опроса или анкетирования, скорее всего, будут содержать возраст пациента и время (календарную дату) проведения опроса). То же самое относится и к наборам данных, содержащим только количественные переменные: как правило, присутствует порядковый номер образца, и/или идентификатор пациента и/или гендерный признак пациента. Рассмотрим на примерах некоторые наборы данных, содержащиеся в пакетах языка R, в частности в пакете MASS. Знак «#» применяется для экранирования (компилятор языка не будет воспринимать текст, находящийся после данного знака, как текст программы) комментариев в тексте программы.

Листинг 141

library (MASS) #Подключаем библиотеку, содержащую набор данных bacteria

head (bacteria) # выводим первую часть набора данных

#=======================================================

# Результат вывода первой части набора данных

#=====================================================

y ap hilo week ID trt

1 y p hi 0 X01 placebo

2 y p hi 2 X01 placebo

3 y p hi 4 X01 placebo

4 y p hi 11 X01 placebo

5 y a hi 0 X02 drug+

6 y a hi 2 X02 drug+

#=======================================================

#Проводим определение структуры набора данных

#=======================================================

str (bacteria) # Выводим структуру данных, содержащуюся в наборе данных

# bacteria

#=======================================================

# Результаты применения функции

#=======================================================

'data.frame’: 220 obs. of 6 variables:

$ y: Factor w/ 2 levels «n», «y»: 2 2 2 2 2 2 1 2 2 2…

$ ap: Factor w/ 2 levels «a», «p»: 2 2 2 2 1 1 1 1 1 1…

$ hilo: Factor w/ 2 levels «hi», «lo»: 1 1 1 1 1 1 1 1 2 2…

$ week: int 0 2 4 11 0 2 6 11 0 2…

Продолжение листинга 1

$ ID: Factor w/ 50 levels «X01», «X02», «X03»,..: 1 1 1 1 2 2 2 2 3 3…

$ trt: Factor w/ 3 levels «placebo», «drug»,..: 1 1 1 1 3 3 3 3 2 2…

Приведенный набор данных содержит в основном «факторы», его уровни представлены естественным языком. Таким представлением в языке R описываются качественные данные, а сам набор данных содержит только один количественный показатель (week), описывающий порядковый номер недели проведения исследования.

Листинг 2

library (MASS) #Подключаем пакет, содержащий набор данных muscle

head (muscle) #Выводим начало набора данных

#=======================================================

#Результат применения команды

#=======================================================

Strip Conc Length

3 S01 1 15.8

4 S01 2 20.8

5 S01 3 22.6

6 S01 4 23.8

9 S02 1 20.6

10 S02 2 26.8

#=======================================================

#Выводим структуру набора данных

#=======================================================

str (muscle)

#=======================================================

'data.frame’: 60 obs. of 3 variables:

$ Strip: Factor w/ 21 levels «S01», «S02», «S03»,..: 1 1 1 1 2 2 2 2 3 3…

Продолжение листинга 2

$ Conc: num 1 2 3 4 1 2 3 4 0.25 0.5…

$ Length: num 15.8 20.8 22.6 23.8 20.6 26.8 28.4 27 7.2 15.4…

Набор данных muscle содержит значения концентрации хлорида кальция, кратные величине 2,2 ммоль (Conc), и длину полосы сокращения сердечной мышцы, измеряемую в миллиметрах (Length) – представленные данные являются количественными. Одновременно набор данных содержит и качественный параметр, а именно обозначение линии сердечной мышцы при проведении измерений (Strip).

Каждая переменная, входящая в набор данных, имеет совокупность значений, полученных в результате наблюдения или путем прямого измерения. Взаимная зависимость и влияние изменения одной величины на другую являются предметом исследования, в том числе методами статистического анализа.

Если на исследуемую величину действует большое количество независимых переменных и отсутствует возможность выделить доминирующую величину, то распределение исследуемой величины стремится к так называемому нормальному или Гауссову закону распределения (при проведении исследования данных выбор применяемого для анализа теста или критерия зависит от того, принадлежит ли исследуемая величина нормальному закону распределения или нет).

На рисунке 3 представлен пример графика (гистограммы42) случайной величины, значительно приближенной к нормальному закону распределения.

Рисунок 3 – Пример частоты встречаемости наблюдаемого и/или измеренного значения для нормального закона распределения

Опишем нахождение основных метрик и охарактеризуем их для данных, распределенных в соответствии с нормальным законом.

2.1. Базовый статистический анализ количественных данных

Базовый статистический анализ составляет основу любой аналитической работы с данными (независимо от того, какого типа данные анализируются).

В случае, если набор данных содержит в основном качественные данные, то при проведении базового статистического анализа определяют количество данных, содержащихся в каждом факторе, количество уровней факторов и количество данных, содержащихся в каждом уровне. На основании этих данных вычисляются:

1. Доля каждого уровня в факторе.

2. Среднее квадратичное отклонение доли в факторе.

В случае, если набор данных содержит в основном количественные данные, то при проведении базового статистического анализа определяют:

1. Выборочное среднее и/или среднее взвешенное значение.

2. Медиану.

3. Среднее квадратическое отклонение.

4. Доверительный интервал выборочного среднего на основании первой и третьей величины.

5. Максимальное значение.

6. Минимальное значение.

7. «Размах» значений на основании максимального и минимального значений.

8. Первый квартиль в распределении данных.

9. Последний квартиль в распределении данных.

10. Межквартильный интервал на основании первого и последнего квартиля в распределении данных.

Перечисленные величины описывают основные статистические свойства исследуемых данных, но не дают ответа на вопрос о наличии различий в группах данных и их взаимном влиянии. Рассмотрим более детально способ вычисления величин, входящих в перечень базовых статистических величин.

2.1.1. Выборочное среднее

При проведении анализа данных первая величина, которую вычисляет исследователь – это выборочное среднее значение43. В данном учебном пособии не рассматривается вопрос математически строгого представления средних величин, а представлены наиболее часто использующиеся на практике. Выборочная средняя величина вычисляется по уравнению (1):

где N – количество исследований, вычисленное по одному параметру в наборе данных; X_i – фактическое значение усредняемой величины.

Она является обобщенной характеристикой в случае однородности данных, описывающей явления, имеющие одну и ту же размерность. Например, если врач проводит исследование веса или роста пациентов, возраст которых составляет 20 полных лет, то ему необходимо описать вес или рост всех пациентов в возрасте 20 полных лет. Для этих целей он будет использовать выборочное среднее значение. В случае данных, распределенных в соответствии с нормальным законом распределения, выборочное среднее значение будет соответствовать максимуму частоты встречаемости исследуемого значения, как представлено на рисунке 4.

Рисунок 4 – Выборочное среднее (зеленая вертикальная линия)

для нормального закона распределения44

На языке программирования R выборочное среднее значение вычисляется с помощью функции mean () из пакета base.

В случае, если врачу необходимо более точно описать средний вес пациентов (например) в возрасте 20 лет с учетом их (например) роста, то для этих целей нужно использовать среднее взвешенное значение (2):

где N – количество исследований, вычисленное по одному параметру в наборе данных; X_i – фактическое значение усредняемой величины; ω_i – вес i-го значения усредняемой величины в общей выборке.

В качестве весового коэффициента ω_i в этом случае будет выступать величина, описывающая рост i-го пациента в возрасте 20 лет. Средняя взвешенная величина выступает в качестве характеристики, описывающей совокупность явлений, имеющих одну и ту же размерность с учетом влияния сторонних признаков (например, вес пациентов с учетом их роста).

На языке программирования R выборочное среднее значение вычисляется с помощью функции weighted.mean () из пакета base.

При исследованиях долей возникновения признака может возникнуть необходимость вычисления средней доли признака по отношению ко всем признакам. В этом случае среднее значение доли вычисляется как (3):

где p – доля единиц, обладающих исследуемым признаком; q – доля единиц, не обладающих исследуемым признаком, равна 1-p.

Пример вычисления выборочного среднего и средневзвешенного на языке R

Для примера вычисления средних и средневзвешенных значений используем набор данных anorexia из пакета MASS. Вычислим среднее и средневзвешенное значение веса пациента. В качестве весового коэффициента (уравнение 2 и текст под ним) будем использовать долю пациентов в каждой группе.

Листинг 3

library (MASS) #45Подключаем библиотеку, содержащую набор данных anorexia

#=====================================================

#Исследуем структуру набора данных

#=====================================================

str (anorexia)

#=======================================================

# Выводим структуру данных

#=======================================================

'data.frame’: 72 obs. of 3 variables:

$ Treat: Factor w/ 3 levels «CBT», «Cont», «FT»: 2 2 2 2 2 2 2 2 2 2..

Продолжение листинга 3

$ Prewt: num 80.7 89.4 91.8 74 78.1 88.3 87.3 75.1 80.6 78.4…

$ Postwt: num 80.2 80.1 86.4 86.3 76.1 78.1 75.1 86.7 73.5 84.6…

#========================================================

# Выводим заголовочную часть набора данных

#========================================================

head (anorexia)

#=======================================================

# Выводим заголовочную часть файла

#========================================================

Treat Prewt Postwt

1 Cont 80.7 80.2

2 Cont 89.4 80.1

3 Cont 91.8 86.4

4 Cont 74.0 86.3

5 Cont 78.1 76.1

6 Cont 88.3 78.1

#========================================================

# Вычисляем выборочное среднее в контрольной группе пациентов до

# проведения лечения

#=====================================================

mean (anorexia [anorexia$Treat==«Cont», ] $Prewt)

#========================================================

# результат вычисления среднего веса пациента в контрольной группе

#========================================================

[1] 81.55769

#========================================================

mean (anorexia [anorexia$Treat==«CBT», ] $Prewt)

#========================================================

# результат вычисления среднего веса пациента в группе с

Продолжение листинга 3

# когнитивно-поведенческой терапией

#========================================================

[1] 82.68966

#========================================================

mean (anorexia [anorexia$Treat==«FT», ] $Prewt)

#========================================================

# результат вычисления среднего веса пациента в группе с семейной терапией

#========================================================

[1] 83.22941

#========================================================

#Вычисление взвешенного среднего значения

#========================================================

# Определяем количество пациентов в каждой группе и общее количество пациентов

numInCont <– length46 (anorexia [anorexia$Treat==«Cont», ] $Prewt)

numInCBT <– length (anorexia [anorexia$Treat==«CBT», ] $Prewt)

numInFT <– length (anorexia [anorexia$Treat==«FT», ] $Prewt)

numAll <– length (anorexia$Prewt)

#=======================================================

# Формируем вектор весов и вектор данных

# В данном примере в качестве веса выступает доля пациентов в каждой

# группе

#=======================================================

weightData <– c (numInCont/numAll, numInCBT/numAll, numInFT/numAll)

Продолжение листинга 3

meanData <– c (mean (anorexia [anorexia$Treat==«Cont», ] $Prewt),

mean (anorexia [anorexia$Treat==«CBT», ] $Prewt),

mean (anorexia [anorexia$Treat==«FT», ] $Prewt)

weighted.mean (meanData, weightData)

#=======================================================

# Результат вычисления средневзвешенного значения веса пациента

#=======================================================

[1] 82.40833

Полученные значения среднего и средневзвешенного веса пациентов отражают типичный уровень веса пациента в исследуемой группе, формирующийся под воздействием доминирующих неслучайных параметров. Установление неслучайных параметров, влияющих (в данном случае на вес пациента) на среднее значение, требует дальнейшего исследования данных с использованием соответствующих статистических критериев, в том числе критериев, позволяющих сравнивать средние величины (например, t-критерий Стьюдента). Оценка среднего отклонения отдельных значений (например, веса пациента) от среднего значения описывается с помощью среднего квадратического отклонения.

2.1.2. Вычисление среднего квадратического отклонения

Среднее квадратическое отклонение является статистической характеристикой, показывающей степень близости отдельных значений к средней величине (чем меньше среднее квадратическое отклонение, тем лучше средняя величина описывает исследуемую переменную). Среднее квадратическое отклонение генеральной совокупности значений вычисляется как (4):

где N – полное количество исследований в генеральной совокупности; x̅ – среднее значение исследуемой величины; x_i – фактическое значение, полученное при измерении величины.

В большинстве случаев исследователь имеет дело не с генеральной совокупностью (всей популяцией), а с выборочной ее частью. Среднее квадратическое отклонение для выборки из генеральной совокупности вычисляется как (5):

Оценки верны для случаев, когда исследуемые переменные имеют количественную природу, и данные принадлежат к нормальному закону распределения. Иллюстрация среднего квадратического отклонения представлена на рисунке 5.

Рисунок 5 – Белой вертикальной линией отмечено выборочное среднее значение данных, распределенных по нормальному (Гауссову) закону. Расстояние от белой вертикальной линии до зеленой (влево или вправо) задается средним квадратическим отклонением

В случае, если исследователю необходимо оценить среднее квадратическое отклонение вероятности возникновения признака, то составляется уравнение (6):

где p̂ – вероятность возникновения исследуемого признака:

где N – общее количество значений, обладающих различными признаками; n – количество исследований, обладающих анализируемым признаком.

На языке программирования R для вычисления среднего квадратического значения применяется функция sd () из пакета stats. В случае данных, распределенных по нормальному закону, выборочное среднее квадратическое отклонение будет соответствовать рисунку 5.

Пример вычисления среднего квадратического отклонения

Для примера проведения вычислений среднего квадратического отклонения используем набор данных anorexia из пакета MASS.

Листинг 4

library (MASS) #Подключаем пакет MASS, содержащий набор данных anorexia

#========================================================

#Выделяем из набора данных подгруппы пациентов с разным типом лечения

#========================================================

contData <– anorexia [anorexia$Treat==«Cont», ] $Prewt

cbtData <– anorexia [anorexia$Treat==«CBT», ] $Prewt

ftData <– anorexia [anorexia$Treat==«FT», ] $Prewt

#=====================================================

#Вычисляем среднее квадратическое отклонение веса пациента до

# прохождения лечения

Продолжение листинга 4

#=====================================================

sd (contData)

#====================================================

# Результат вычисления в контрольной группе

#====================================================

[1] 5.70706

#====================================================

sd (cbtData)

#====================================================

# Результат вычисления в группе с когнитивно-поведенческой терапией

#=====================================================

[1] 4.845495

#========================================================

sd (ftData)

#=====================================================

# Результаты вычислений в группе с семейной терапией

#=====================================================

[1] 5.016693

#=====================================================

# Вычисляем количество пациентов до прохождения лечения

#=====================================================

lenPrewt <– length (anorexia$Prewt) #Общее количество данных

lenContPrew <– length (anorexia [anorexia$Treat==«Cont», ] $Prewt) # Количество

# данных в контрольной группе

lenCBTPrewt <– length (anorexia [anorexia$Treat==«CBT», ] $Prewt) # Количество

# данных в группе с когнитивно-поведенческой терапией

lenBTPrew <– length (anorexia [anorexia$Treat==«FT», ] $Prewt) # Количество

# данных в группе с семейной терапией

Продолжение листинга 4

#========================================================

# Вычисляем доли пациентов в каждой группе

#=====================================================

pCP <– lenContPrew/ lenPrewt # Доля в контрольной группе

pCBTP <– lenCBTPrewt/ lenPrewt # Доля в когнитивно-поведенческой группе

pBTP <– lenBTPrew/lenPrewt # Доля в семейной группе

#========================================================

sdCP <– sqrt47 (pCP* (1-pCP)) # Среднее квадратичное отклонение долей в

# контрольной группе.

sdCBT <– sqrt (pCBTP* (1- pCBTP))

sdBT <– sqrt (pBTP* (1- pBTP))

#=====================================================

# Результаты вычисления

#=====================================================

print (sdCP)

#=====================================================

[1] 0.4803227

#=====================================================

print (sdCBT)

#=====================================================

[1] 0.4904568

#=====================================================

print (sdBT)

#=====================================================

[1] 0.4246912

2.1.3. Вычисление доверительного интервала

Вычислив значения выборочной средней величины (или средневзвешенной величины)48 и среднего квадратического отклонения, можно построить доверительный интервал средней величины. Вычисление доверительного интервала для средней величины тесно связано с понятием доверительной вероятности, ошибки средней величины и предельной ошибкой выборки.

Доверительная вероятность определяет степень уверенности того факта, что измеренная величина находится вблизи среднего значения. Средняя ошибка выборки показывает объективно возникающее расхождение между характеристиками выборки и генеральной совокупностью, а предельная ошибка выборки – это ошибка выборки, исчисляемая с заданной степенью вероятности. Ошибка средней величины в выборке для количественного признака вычисляется по уравнению (8):

где S – выборочное среднее квадратическое отклонение; N – количество исследований.

Ошибка средней величины для долей вхождения признака в данные также вычисляется по уравнению (8) с той лишь разницей, что вместо выборочного среднего квадратического отклонения в уравнение подставляется среднее квадратическое отклонение долей исследуемого признака (уравнение 5).

Предельная ошибка выборки вычисляется по уравнению (9):

где значение t представлено в таблице 4.

В медицинских исследованиях доверительный интервал принимается равным (10)49:

где X̅ – среднее значение исследуемой величины; Δ_lim – предельная ошибка выборки исследуемой величины.

В случае оценки доли фактора и отклонения долей величина доверительного интервала вычисляется по уравнению (11):

где X̅_e – средняя доля исследуемого признака, определенная по уравнению (3); Δ^e_lim – предельная ошибка выборки долей вхождения признака в данные.

На рисунке 6 изображено графическое представление доверительного интервала при условии данных, распределенных в соответствии с нормальным (Гауссовым) законом.

Рисунок 6 – Доверительный интервал (уравнение 10) выборочного среднего значения (расстояние от белой вертикальной линии до зеленой вертикальной линии)

Пример вычисления доверительных интервалов

На примере значений веса пациентов с подтвержденным диагнозом «анорексия» (набор данных anorexia пакета MASS) проведем анализ доверительных интервалов весов и долей пациентов в контрольной группе, группе с когнитивно-поведенческой терапией и в группе с назначенной семейной терапией.

Листинг 5

library (MASS) # Подключаем библиотеку, содержащую набор данных anorexia

library (ggplot2)

#========================================================

# Выделяем из набора данных подгруппы пациентов с разным типом лечения

#========================================================

contData <– anorexia [anorexia$Treat==«Cont», ] $Prewt

cbtData <– anorexia [anorexia$Treat==«CBT», ] $Prewt

ftData <– anorexia [anorexia$Treat==«FT», ] $Prewt

#=====================================================

# Вычисляем средние значения

#=====================================================

meanCont <– mean (contData) # Средний вес в контрольной группе

meanCBT <– mean (cbtData) # Средний вес в когнитивно-поведенческой группе

meanFT <– mean (ftData) # Средний вес в группе семейной терапии

#========================================================

# Вычисляем среднее квадратическое отклонение

#========================================================

sdCont <– sd (contData) / (sqrt (length (contData))) # Среднее квадратическое

# отклонение веса в контрольной группе

Продолжение листинга 5

sdCBT <– sd (cbtData) / (sqrt (length (cbtData))) # Среднее квадратическое отклонение

# веса в когнитивно-поведенческой группе

sdFT <– sd (ftData) / (sqrt (length (ftData))) # Среднее квадратическое отклонение

# веса в группе семейной терапии

#=====================================================

# Вычисление доверительных интервалов в группах

#=====================================================

upCont <– meanCont+2*sdCont # Верхняя граница доверительного интервала в

# контрольной группе

downCont <– meanCont-2*sdCont # Нижняя граница доверительного интервала в

# контрольной группе

#=====================================================

upCBT <– meanCBT+2*sdCBT # Верхняя граница доверительного интервала в

# группе с когнитивно-поведенческой терапией

downCBT <– meanCBT-2*sdCBT # Нижняя граница доверительного интервала в

# группе с когнитивно-поведенческой терапией

#=====================================================

upFT <– meanFT+2*sdFT # Верхняя граница доверительного интервала в

# группе с семейной терапией

downFT <– meanFT-2*sdFT # Нижняя граница доверительного интервала в

# группе с семейной терапией

#=====================================================

# Построение графика средних значений и доверительных интервалов

#=====================================================

Продолжение листинга 5

dataAn <– data.frame (Type = c («Cont», «CBT», «FT»), Mean=c (meanCont, meanCBT, meanFT), Sd = c (sdCont, sdCBT, sdFT))

#=======================================================

grPP <– ggplot (data = dataAn, mapping = aes (Type, Mean))

grPP <– grPP + geom_point (mapping = aes (Type, Mean),colour=«blue», size=5)

grPP <– grPP + geom_errorbar (mapping = aes (ymin=Mean-2*Sd, ymax=Mean+2*Sd), colour=«red», linewidth=1.5,width = 0.5)

grPP <– grPP + theme_light ()

grPP <– grPP + theme(panel.border = element_rect (linewidth = 6, colour = «black»),

axis. h2 = element_text (size = 20, face = «bold», colour = «black»),

axis. text. x = element_text (size = 20, face = «bold», colour = «black»),

axis. text. y = element_text (size = 20, face = «bold», colour = «black»),

axis. ticks = element_line (linewidth = 3, colour = «black»),

legend. h2 = element_text (size = 20, face = «bold», color = «black»),

legend. text = element_text (size = 20, face = «bold», color = «black»))

grPP <– grPP + labs (x=«Therapy type», y=«Weight, kg»)

print (grPP)

#=====================================================

# Вывод графика

#=====================================================

На рисунке 7 изображено графическое представление средних значений веса с доверительными интервалами трех исследуемых групп пациентов до проведения терапии.

Рисунок 7 – Средние значения (синие кружочки) и доверительный интервал для доверительной вероятности 95,5% (красные линии, нижние и верхние горизонтальные линии обозначают границы доверительных интервалов) веса пациентов в трех группах (Cont – контрольной группе; CBT – группе с назначенной когнитивно-поведенческой терапией; FT – с назначенной семейной терапией) до проведения терапии

На примере того же набора данных продемонстрируем вычисление вероятности нахождения пациента с определенным весом в определенной группе и вычислим доверительный интервал для этой вероятности.

Листинг 6

# Вычисляем количество пациентов до прохождения лечения

#=====================================================

lenPrewt <– length (anorexia$Prewt) #Общее количество данных

Продолжение листинга 6

lenContPrew <– length (anorexia [anorexia$Treat==«Cont», ] $Prewt) # Количество

# данных в контрольной группе

lenCBTPrewt <– length (anorexia [anorexia$Treat==«CBT», ] $Prewt) # Количество

# данных в группе с когнитивно-поведенческой терапией

lenBTPrew <– length (anorexia [anorexia$Treat==«FT», ] $Prewt) # Количество

# данных в группе с семейной терапией

#=====================================================

# Вычисляем доли пациентов в каждой группе

#=====================================================

pCP <– lenContPrew/ lenPrewt # Доля в контрольной группе

pCBTP <– lenCBTPrewt/ lenPrewt # Доля в когнитивно-поведенческой группе

pBTP <– lenBTPrew/lenPrewt # Доля в семейной группе

#=====================================================

sdCP <– sqrt (pCP* (1-pCP)) # Среднее квадратичное отклонение долей в

# контрольной группе

sdCBT <– sqrt (pCBTP* (1- pCBTP))

sdBT <– sqrt (pBTP* (1- pBTP))

#=====================================================

# Результаты вычисления

#=====================================================

dataAn <– data.frame (Type = c («Cont», «CBT», «FT»), Mean=c (pCP, pCBTP, pBTP),

Sd = c (sdCP, sdCBT, sdBT))

#=====================================================

grPP <– ggplot (data = dataAn, mapping = aes (Type, Mean))

grPP <– grPP + geom_point (mapping = aes (Type, Mean),colour=«blue», size=5)

Продолжение листинга 6

grPP <– grPP + geom_errorbar (mapping = aes (ymin=Mean-2*Sd, ymax=Mean+2*Sd), colour=«red», linewidth=1.5,width = 0.5)

grPP <– grPP + theme_light ()

grPP <– grPP + theme(panel.border = element_rect (linewidth = 6, colour = «black»),

axis. h2 = element_text (size = 20, face = «bold», colour = «black»),

axis. text. x = element_text (size = 20, face = «bold», colour = «black»),

axis. text. y = element_text (size = 20, face = «bold», colour = «black»),

axis. ticks = element_line (linewidth = 3, colour = «black»),

legend. h2 = element_text (size = 20, face = «bold», color = «black»),

legend. text = element_text (size = 20, face = «bold», color = «black»))

grPP <– grPP + labs (x=«Therapy type», y=«Proportion of patients»)

print (grPP)

#=====================================================

# Вывод графика

#=====================================================

На рисунке 8 графически изображены средние значения долей пациентов со средними квадратичными отклонениями в каждой из исследуемых групп пациентов до прохождения лечения.

Рисунок 8 – Доля пациентов в каждой группе (синяя точка) и доверительный интервал для 95,5% доверительной вероятности (красные линии обозначают ширину доверительного интервала, горизонтальными линиями обозначены границы доверительного интервала) доли пациентов в каждой группе пациентов, с назначенными различными типами терапии (Cont – контрольной группе; CBT – группе с назначенной когнитивно-поведенческой терапией; FT – с назначенной семейной терапией)

Сравнение доверительных интервалов и долей показывает, что все три группы достаточно близки друг к другу, что позволяет выдвинуть гипотезу об отсутствии статистически значимых различий. Данная гипотеза требует дальнейшего подтверждения или опровержения посредством статистических критериев.

2.1.4. Поиск максимального, минимального значения и размах

Анализ максимального значения, минимального значения и размаха количественной величины дополняет анализ средних значений и доверительных интервалов в части оценки фактического разброса значений. Рассмотрим оценку максимального, минимального значений и размаха количественной величины на примере анализа веса пациентов с диагнозом «анорексия» в трех группах с различным видом назначенного лечения.

Листинг 7

library (MASS) #Подключаем библиотеку, содержащую набор данных anorexia

library (ggplot2) #Подключаем библиотеку построения графиков

#=====================================================

# Выделяем из набора данных подгруппы пациентов с разным типом лечения

#=====================================================

contData <– anorexia [anorexia$Treat==«Cont», ] $Prewt

cbtData <– anorexia [anorexia$Treat==«CBT», ] $Prewt

ftData <– anorexia [anorexia$Treat==«FT», ] $Prewt

#=====================================================

# Определяем максимальные значения в данных

#=====================================================

maxCont <– max50 (contData) # Максимальный вес до проведения лечения

# в контрольной группе пациентов

maxCBT <– max (cbtData) # Максимальный вес до проведения лечения

# в группе с назначенной когнитивно-поведенческой терапией

maxFT <– max (ftData) # Максимальный вес до проведения лечения

# в группе с назначенной семейной терапией

Продолжение листинга 7

#=====================================================

# Определяем минимальные значения в данных

#=====================================================

minCont <– min51 (contData) # Минимальный вес пациентов до прохождения лечения

# в контрольной группе

minCBT <– min (cbtData) # Минимальный вес пациентов до прохождения лечения

# в группе с назначенной когнитивно-поведенческой терапией

minFT <– min (ftData) # Минимальный вес пациентов до прохождения лечения

# в группе с назначенной семейной терапией

#=====================================================

# Определяем размах значений веса в каждой группе

#=====================================================

deltaCont <– maxCont – minCont # Размах значений веса пациентов до прохождения

# лечения в контрольной группе

deltaCBT <– maxCBT – minCBT # Размах значений веса пациентов до прохождения

# лечения в группе с назначенной когнитивно-поведенческой терапией

deltaFT <– maxFT – minFT # Размах значений веса пациентов до прохождения

# лечения в группе с назначенной семейной терапией

#=====================================================

# Группируем результаты вычислений

#=====================================================

Продолжение листинга 7

dataAn <– data.frame (Type = c («Cont», «CBT», «FT»),

Max=c (maxCont, maxCBT, maxFT),

Min = c (minCont, minCBT, minFT),

Delta = c (deltaCont, deltaCBT, deltaFT))

#=====================================================

grPP <– ggplot (data = dataAn, aes (Type, Max))

grPP <– grPP + geom_point (mapping = aes (Type, Min),colour=«magenta», size=7)

grPP <– grPP + geom_point (mapping = aes (Type, Max),colour=«blue», size=7)

grPP <– grPP + geom_errorbar (mapping = aes (ymin=Min, ymax=Max), colour=«red», linewidth=1.5,width = 0.5)

grPP <– grPP + theme_light ()

grPP <– grPP + theme(panel.border = element_rect (linewidth = 6, colour = «black»),

axis. h2 = element_text (size = 20, face = «bold», colour = «black»),

axis. text. x = element_text (size = 20, face = «bold», colour = «black»),

axis. text. y = element_text (size = 20, face = «bold», colour = «black»),

axis. ticks = element_line (linewidth = 3, colour = «black»),

legend. h2 = element_text (size = 20, face = «bold», color = «black»),

legend. text = element_text (size = 20, face = «bold», color = «black»))

grPP <– grPP + labs (x=«Therapy type», y=«Value range»)

print (grPP)

На рисунке 9 графически изображены размах значений массы тела в каждой из исследуемых групп пациентов до прохождения лечения.

Рисунок 9 – Размах значений веса пациентов (синие закрашенные кружочки – максимальный вес; фиолетовые кружочки – минимальный вес пациентов) до проведения лечения в трех группах (CBT – группа с назначенным когнитивно-поведенческим лечением; Cont – контрольная группа; FT – назначенная семейная терапия)

Сравнение размаха значений и доверительных интервалов позволяет определить наличие выбросов в измерениях. Значения выше или ниже доверительного интервала должны быть проанализированы отдельно от основной выборки и не анализироваться в основной совокупности данных.

2.1.5. Понятие о квантилях, децилях, квартилях распределения

При проведении статистического анализа данных возникает необходимость оценки значения, соответствующего 5%52 от общего количества данных или 95% от общего количества данных. Для этих целей применяются различные уровни квантилей.

Квантили – это значения, которые делят упорядоченную выборку на равные доли.

Допустим, имеется произвольный ряд десятичных чисел:

85.50 99.69 37.19 47.43 40.26 93.58 68.97 19.50 87.39 43.81 42.96 6.985.22 62.67 47.16 76.30 69.56 59.40 92.33 31.22 41.59 33.27 30.25 36.7749.43 33.85 27.44 8.45 99.84 52.81

Нам необходимо определить, какое значение соответствует 5% в представленном ряду значений. Для решения данной задачи необходимо выстроить числа ряда в порядке возрастания (от меньшего к большему):

5.21 6.98 8.45 19.50 27.44 30.25 31.22 33.27 33.85 36.77 37.18 40.26 41.59 42.96 43.81 47.16 47.43 49.43 52.81 59.40 62.67 68.97 69.56 76.30 85.50 87.39 92.33 93.58 99.69 99.84

В данном ряду содержится 30 значений, номер значения, соответствующего 5% данного ряда, равен 1.5, т.е. между первым и вторым значением соответственно квантиль будет равен:

Децили – значения, которые делят упорядоченную выборку на десять примерно равных частей. Допустим, необходимо разделить ряд, представленный выше, на 10 равных частей с шагом в 10%, тогда децили будут равны (процедура поиска значения, соответствующего каждому проценту, такая же, как в описанном выше примере):

Наиболее часто применяемыми на практике являются квартили – значения, которые делят упорядоченную выборку на четыре примерно равные части. Для приведенного ранее ряда значения квартилей равны:

На языке программирования R вычисление квартилей проводится с помощью функции quantile (), входящей в пакет stats. На практике наиболее часто применяются первый и третий квартили, межквартильный интервал представлен на рисунке 10.

Рисунок 10 – Межквартильный интервал (между первым и третьим квартилем) представлен зелеными вертикальными линиями

В межквартильном интервале находятся 50% всех встречаемых значений измеряемой величины. Первый квартиль ограничивает слева 25% (первая зеленая вертикальная линия при просмотре рисунка 10 слева направо) значений, а третий – 75% (вторая зеленая вертикальная линия при просмотре рисунка 10 слева направо).

Пример вычисления квантилей, децилей и квартилей

На примере набора данных anorexia из пакета MASS, содержащих информацию о весе пациентов с подтвержденным диагнозом «анорексия», рассмотрим вычисления квантилей, децилей и квартилей массы тела пациентов в контрольной группе, в группе с назначенной когнитивно-поведенческой терапией и семейной терапией.

Листинг 8

library (MASS) # Подключаем библиотеку, содержащую набор данных anorexia

Продолжение листинга 8

library (ggplot2) #Подключаем библиотеку построения графиков

#=======================================================

# Выделяем из набора данных подгруппы пациентов с разным типом

# назначенного лечения

#=======================================================

contData <– anorexia [anorexia$Treat==«Cont», ] $Prewt #Вес пациентов в

# контрольной группе до проведения лечения

cbtData <– anorexia [anorexia$Treat==«CBT», ] $Prewt # Вес пациентов в группе

# с назначенным когнитивно-поведенческим лечением до проведения лечения

ftData <– anorexia [anorexia$Treat==«FT», ] $Prewt # Вес пациентов в группе

# с назначенной семейной терапией

#=======================================================

#Задаем значения интересующих нас квантилей

propQ <– c (0.05, 0.95)

propDec <– c (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0)

propQvar <– c (0.25, 0.5, 0.75, 1.0)

#=======================================================

contQuan <– quantile (contData, propQ, type = 4)

print (contQuan)

cbtQuan <– quantile (cbtData, propQ, type = 4)

print (cbtQuan)

ftQuan <– quantile (ftData, propQ, type = 4)

print (ftQuan)

2.1.6. Вычисление медианы

Выборочное среднее значение при большом объеме данных, распределенных в соответствии с нормальным законом, стремится к моде53 или наиболее вероятному значению в данной выборке. Такая оценка значений не всегда бывает адекватной, особенно в случае наличия аномалий (значений существенно меньших или существенно больших, чем большинство значений из выборочной совокупности). Для более адекватного описания всей совокупности при наличии выбросов или данных, распределенных отлично от нормального закона распределения, лучше использовать медианную оценку, или медиану – это значение, которое разделяет ранжированную по возрастанию выборку пополам (50% квантиль).

Предположим, имеется упорядоченный по возрастанию ряд чисел 1, 4, 6, 9, 11. Центр этого ряда составляет число 6, это и будет медианой данного ряда. В случае, если ряд имеет четное число значений, то медианой будет являться среднее арифметическое значение между двумя центральными значениями. Например, имеется ряд чисел, упорядоченных по возрастанию: 1, 4, 6, 9, 11, 12 – двумя центральными значениями будут являться 6 и 9, их среднее арифметическое значение равно 7,5 – оно и будет являться медианой.

Пример вычисления медианы

Рассмотрим вычисление медианы на наборе данных anorexia из пакета MASS. Рассчитаем медиану веса пациентов до прохождения лечения в контрольной группе, группе пациентов с предписанной когнитивно-поведенческой терапией и в группе с семейной терапией.

Листинг 9

library (MASS) # Подключаем пакет, содержащий набор данных anorexia

#Вычисляем медиану в контрольной группе пациентов

median (anorexia [anorexia$Treat==«Cont», ] $Prewt)

#=======================================================

# Результат вычисления

#=======================================================

[1] 80.65

#Вычисляем медиану в группе с когнитивно-поведенческой терапией

median (anorexia [anorexia$Treat==«CBT», ] $Prewt)

#=======================================================

#Результат вычисления

#=======================================================

[1] 82.6

#Вычисляем медиану в группе с семейной терапией

median (anorexia [anorexia$Treat==«FT», ] $Prewt)

#=======================================================

#Результат вычисления

#=======================================================

[1] 83.3

Результаты вычислений показывают, что медианная оценка веса пациентов с подтвержденным диагнозом «анорексия» во всех трех группах примерно равна выборочному среднему значению (см. пример в п. 2.1.1), что согласуется с одним из свойств нормального закона распределения количественной величины – медиана, среднее и мода приблизительно равны между собой.

Все значения проведенного базового статистического анализа обобщаются в виде таблицы значений. В таблице 5 представлен пример сводной таблицы базового статистического анализа, проведенного для пациентов с подтвержденным диагнозом anorexia.

Хорошей иллюстрацией для обобщения базового статистического анализа является так называемый ящик с усами – диаграмма размаха, на которой в графическом виде представлен базовый статистический анализ результатов, рассчитанный на основе межквартильного интервала (рисунок 11).

Рисунок 11 – Диаграмма размаха «ящик с усами» веса пациентов: а – до прохождения лечения; б – после прохождения лечения

Центральная (жирная линия) ящика означает медиану, нижняя граница ящика (прямоугольника) означает первый квартиль распределения значений, верхняя граница ящика (прямоугольника) – третий квартиль распределения исследуемой величины (в данном случае веса пациентов), «усы ящика» вычисляются по уравнениям (12, 13):

где U_down – значение нижней границы «усов ящика»; – Q₁ первый квартиль; IQR – межквартильный интервал.

где U_up – значение верхней границы усов ящика; Q₃ – третий квартиль; IQR – межквартильный интервал.

Точки, выходящие за пределы длины «усов», могут считаться аномальными (значения, сильно отличающиеся в большую или меньшую сторону от основной совокупности значений).

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

Контрольные вопросы и задачи

1. Напишите уравнение для вычисления выборочного среднего значения.

2. Запишите уравнение для вычисления выборочного среднего квадратического отклонения.

3. Опишите способ поиска медианы в числовом векторе, расположенном в случайном порядке.

4. Используя набор данных Indometh и Theoph из пакета MASS, вычислите медиану, среднее квадратическое отклонение и выборочное среднее значение концентраций индометацина и теофиллина в плазме крови пациентов.

Одним из подходов предварительной формулировки нулевой гипотезы является построение гистограмм распределения случайной величины. На рисунке 13 представлена гистограмма распределения абстрактной случайной величины X.

Рисунок 13 – Гистограмма распределения абстрактной случайной величины X

По виду диаграммы можно сделать предположение о близости типа распределения случайной величины к нормальному закону. Однако не стоит забывать о том, что данное предположение будет являться нулевой гипотезой, которая требует дальнейшей проверки59.

Требование проверки нулевой гипотезы связано с особенностью построения гистограммы распределения случайной величины, а именно с выбором ширины интервалов, в которых подсчитывается число попаданий случайных значений. На рисунке 14 представлены изменения гистограммы случайной величины X (рисунок 13).

Рисунок 14 – Изменение вида гистограммы распределения случайной величины X в зависимости от выбора ширины интервалов, в которых подсчитывается количество случайных величин, попавших в данный интервал: а – зауженный интервал; б – интервал увеличенной ширины; в – интервал средней ширины; г – широкий интервал

В результате предварительного анализа гистограммы, представленной на рисунке 14а, можно выдвинуть нулевую гипотезу о принадлежности распределения случайной величины X близко к дискретному типу распределения.

Результатом визуального анализа гистограммы, представленной на рисунке 14б, может служить постановка нулевой гипотезы о мультимодальности распределения случайной величины Х (в этом случае необходимо проверить сложную статистическую гипотезу; рассмотрение анализа сложных статистических гипотез выходит за рамки данного пособия).

Анализ гистограммы, представленной на рисунке 14в, может стать основой для формулирования нулевой гипотезы о принадлежности распределения случайной величины к нормальному закону распределения.

Таким образом, гистограмма является точкой начала отсчета для формулировки нескольких нулевых гипотез, которые должны быть проверены различными статистическими критериями.

4.1. Расчет оптимального интервала на гистограмме

Расчет оптимального количества интервалов и диапазона значений, попадающих в интервал гистограммы распределения, является важной составляющей базовой статистики и может служить основой для выделения групп равномерно распределенных данных и выработки классификации на их основе. В настоящем пособии рассматривается несколько простых подходов расчета оптимального количества интервалов и диапазона изменения значений. Первый подход основан на эвристической формуле Стерджесса (14)60:

где N – объем анализируемой выборки (количество значений); n – количество интервалов разбиения данных.

В работе61, описывающей применение статистического анализа при контроле качества продукции, применяют уравнение Брукса—Каррузера (15):

Наиболее простым методом оценки количества интервалов является (16)62:

При больших объемах данных рекомендуется использовать уравнение (17)63:

После расчета количества интервалов разбиения гистограммы распределения проводят вычисления ширины интервала по уравнению (18):

где max (X) – максимальное значение исследуемой переменной X; min (X) – минимальное значение исследуемой переменной X; n – количество интервалов на гистограмме.

Соответственно, с использованием уравнения (18) вычисляются первый и последующие интервалы по уравнению (19):

где X_i+1 – правая граница i-го интервала, не превышающая max (X); X_i – левая граница i-го интервала, начало отсчета которой равно min (X).

Высота столбца диаграммы определяется по уравнению (20):

где N – общее количество исследований; N_i – количество исследований в i-м интервале; H_i – высота интервала.

Более точная (и она же более сложная) оценка размера интервала может быть получена при оценках статистической мощности64, данный подход не рассматривается в настоящем пособии.

4.2. Построение гистограмм распределения на языке R

В языке программирования R существует несколько способов построения гистограмм распределения случайной величины:

– с применением функции hist (), входящей в пакет graphics;

– посредством функции geom_histogram (), входящей в состав пакета ggplot2;

– с помощью столбчатой диаграммы с предварительным расчетом частот и интервалов распределения величин. Построение столбчатых диаграмм можно произвести с помощью:

– функции barplot (), входящей в пакет graphics;

– функции geom_bar (), входящей в пакет ggplot2.

В примерах данной главы приведены все перечисленные способы построения гистограммы распределений.

Примеры построения гистограмм распределений на языке R

Для построения диаграмм распределения использовались данные по выживаемости заболевших вирусом иммунодефицита, собранных в Австралии после 1 июля 1991 года. Все данные содержатся в пакете MASS, наборе данных Aids2.

Листинг 10

#Построение гистограммы распределения с помощью функции hist () пакета graphics

library («MASS») #Подключаем пакет MASS, содержащий набор данных Aids2

N <– length (Aids2$Age) #Вычисляем количество пациентов в наборе данных

dBreak <-5*log10 (N) #Определяем количество интервалов для гистограммы

dBreak <– round (dBreak) #Округляем значение до целых чисел

dX <– (max (Aids2$Age) -min (Aids2$Age)) / dBreak #Рассчитываем шаг

step <– seq (min (Aids2$Age),max (Aids2$Age), dX) # Вектор интервалов

hist (x = Aids2$age, break = step, xlab= «Возраст выявления, полных лет»,

ylab = «Количество, чел», main = «Выживаемость, заболевших СПИДом»,

col = «blue») #Строим гистограмму распределения возрастов пациентов с

# с выявленным вирусом иммунодефицита

# x – вектор значений исследуемых распределений, полученный из Data Frame

Продолжение листинга 10

#;

# break – количество участков на гистограмме;

# xlab – подпись оси ОХ;

# ylab – подпись оси OY;

# main – название гистограммы;

# col – задание цвета гистограммы распределения.

Результат выполнения данной команды представлен на рисунке 15.

Рисунок 15 – Результат применения функции

hist () из пакета grahics

Функция hist () хороша для применения при необходимости быстро построить гистограмму и посмотреть на распределение данных, однако более гибким инструментом для построения и дальнейшей публикации гистограмм является функция geom_histogram (), входящая в пакет ggplot2.

Листинг 11

# Построение гистограммы распределения с помощью функции

# geom_histogram () пакета ggplot2

library («MASS») # Подключаем пакет MASS, содержащий набор данных Aids2

library («ggplot2») # Подключаем пакет ggplot2, содержащий функцию # geom_histogram ()

p <– ggplot () #Создаем объект p, содержащий65 слои графика

p <– p + geom_histogram (mapping = aes (x= Aids2$age), fill= «blue», binwith=30)

p <– p+labs (x = «Возраст выявления, полных лет»,

y= «Количество, чел»,

h2 = «Выживаемость, заболевших СПИДом»)

print (p)

Результат выполнения кода листинга 11 представлен на рисунке 16.