Подгонка моделей с избыточной дисперсией (экстра-вариация в распределении Пуассона)

******************************************************
***    ПОДГОНКА МОДЕЛЕЙ С ИЗБЫТОЧНОЙ ДИСПЕРСИЕЙ    ***
***  (распределение Пуассона с экстра-вариацией)   ***
******************************************************
* Подобные модели довольно полезны в биологии при изучении популяций паразитов  *
******************************************************

**************************************************************************************************
* По MJ Moroney 'Facts From Figures'                                                             *
* (London: Penguin Books, 1951; 2nd ed 1953)                                                     *
*                                                                                                *
* В пуассоновских моделях E(x)=V(x)=lambda (только один параметр)                                *
* Функция вероятности:                                                                           *
* P(x)=e^(-lambda)*(lambda^x)/x!                                                                 *
*                                                                                                *
* В распределениях с избыточной дисперсией V(X)>E(x)                                             *
* Отрицательный биномиальный закон распределения - прототип распределений с избыточной дисперсией*
*                                                                                                *
* В моделях с отрицательным биномиальным законом требуется дополнительный параметр (k):          *
* Параметр                     Чем оценивается:                                                  *
* E(x)=lambda ----------------> выборочное среднее                                               *
* V(x)=lambda*(1+lambda/k) ---> выборочная дисперсия                                             *
* k = EІ(x)/(V(x)-E(x))      -> выб. средн.^2/(выб. дисп. - выб. средн.) ('метод моментов')      *
* Комментарий: метод максимального правдоподобия лучше, но гораздо более трудоёмкий!             *
* (см. веб и программу ParaDis)                                                                  *
* Функция вероятности:                                                                           *
* P(x)=[(k/(k+lambda))^k]*[(lambda/(k+lambda))^x]*GAMMA(x+k)/(x!*GAMMA(k))                       *
**************************************************************************************************
 
* Пример набора данных (из руководства к программе ParaDis):
* (http://www.bondy.ird.fr/~pichon/paradis/parad2.html).

DATA LIST FREE/ number freq (2 F8.0).
BEGIN DATA
0 70
1 38
2 17
3 10
4  9
5  3
6  2
7  1
END DATA.
VARIABLE LABEL number 'Клещи/лист растения'.
 
* Добавим среднюю, дисперсию и число наблюдений к данным и проведём проверку Фишера на избыточность вариации *.
CACHE.
EXECUTE.
MATRIX.
PRINT /TITLE='ДАННЫЕ И СТАТИСТИКИ'.
GET x/VAR=number.
GET freq/VAR=freq.
COMPUTE ndata=NROW(x).
COMPUTE nt=MSUM(freq).
COMPUTE mean=MSUM(freq&*x)/nt.
COMPUTE variance=(MSUM(freq&*(x&**2))-nt*mean&**2)/(nt-1).
COMPUTE k=(mean&**2)/(variance-mean).
PRINT {x,freq}
 /FORMAT='F8.0'
 /TITLE='Данные выборки'
 /CLABELS='X','Част'.
PRINT nt
 /FORMAT='F8.0'
 /TITLE='Общий объём выборки'
 /RLABEL='Nt'.
PRINT {mean,variance,(variance/mean),k}
 /FORMAT='F8.2'
 /TITLE='Статистики'
 /CLABELS='Среднее','Дисперсия','DP (**)','K(*)'.
PRINT /TITLE='(*) Подсчитан по методу моментов.'.
PRINT /TITLE='(**) Отношение дисперсии к средней.'.
PRINT /TITLE='ПРОВЕРКА ФИШЕРА НА ИЗБЫТОЧНОСТЬ ДИСПЕРСИИ'.
COMPUTE fd=(MSUM(freq&*(x-mean)&**2))/mean.
COMPUTE fdsig=1-CHICDF(fd,MSUM(freq)-1).
PRINT {fd,fdsig}
 /FORMAT='F8.4'
 /TITLE='Статистика проверки (df=Nt-1)'
 /CLABELS='Хи^2','Знач.'.
DO IF fdsig LE 0.05.
- PRINT /TITLE='Присутствует избыточная дисперсия! Обнаружена экстра-вариация в распр. Пуассона'.
END IF.
COMPUTE means=MAKE(ndata,1,mean).
COMPUTE n=MAKE(ndata,1,nt).
COMPUTE kp=MAKE(ndata,1,k).
COMPUTE namevec={'number','n','mean','k'}.
SAVE {x,n,means,kp} /OUTFILE='c:\\temp\\temp.sav' /NAMES=namevec.
END MATRIX.
MATCH FILES /FILE=*
 /FILE='C:\\Temp\\temp.sav'
 /BY number.
EXECUTE.
 
* Подсчёт ожидаемых частот в предположении отрицательного биномиального распределения *
  (Для критеря хи-квадрат происходит потеря 3-х степеней свободы: используется N, среднее и дисперсия) *.
COMPUTE expect=n*((mean/(k+mean))**number)*((k/(k+mean))**k)*EXP(LNGAMMA(number+k))/(EXP(LNGAMMA(number+1))*EXP(LNGAMMA(k))).
COMPUTE sumexp=expect.
* Небольшая хитрость: делаем сумму ожидаемых частот равной N, корректируя ожидаемую частоту последней категории *.
DO IF $casenum GT 1.
- COMPUTE sumexp=sumexp+LAG(sumexp).
END IF.
SORT CASES BY number(D).
DO IF $casenum EQ 1.
- COMPUTE expect=expect+(n-sumexp).
END IF.
SORT CASES BY number(A).
 
* Табличный и графический вывод наблюдаемых и ожидаемых частот *.
VAR LABEL freq 'Наблюдаемые' /expect 'Ожидаемые'.
REPORT
 /FORMAT=LIST
 /TITLE='Наблюдаемые и ожидаемые частоты (с учётом избыточной дисперсии)'
 /VAR=number freq expect.
GRAPH /BAR(GROUPED)=VALUE(freq expect) BY number
 /TITLE='Наблюдаемые и ожидаемые(*) частоты'
 /Footnote='(*) С учётом избыточной дисперсии'.
 
* Проверяем, нет ли малых ожидаемых частот и объединяем ячейки, чтобы избавиться от таких частот *.
COMPUTE id=$casenum.

DO IF (expect LT 2).
- COMPUTE id=LAG(id).
END IF.
AGGREGATE
 /OUTFILE=*
 /BREAK=id
 /observed = SUM(freq)
 /expected = SUM(expect).
STRING groups (A2).
COMPUTE groups = STRING(id,F2.0) .
 
* Проверка согласия по критерию хи-квадрат *.
MATRIX.
PRINT /TITLE='ПРОВЕРКА СОГЛАСИЯ ХИ-КВАДРАТ'.
GET groups/VAR=groups.
GET obs/VAR=observed.
GET expect/VAR=expected.
PRINT {obs,expect,(obs-expect);MSUM(obs),MSUM(expect),0}
 /FORMAT='F8.2'
 /TITLE='Частоты (после объединения категорий)'
 /CLABELS='Наблюдаемые','Ожидаемые','Разность'
 /RNAMES=groups.
COMPUTE k=NROW(obs).
PRINT {k-3}
 /FORMAT='F8.0'
 /TITLE='Число степеней свободы (k-3)'.
COMPUTE chi2=MSUM((obs-expect)&**2/expect).
COMPUTE chisig=1-CHICDF(chi2,k-3).
PRINT {chi2,chisig}
 /FORMAT='F8.4'
 /TITLE='Статистика проверки'
 /CLABELS='Хи^2','Знач.'.
COMPUTE minexp=CMIN(expect).
COMPUTE flag=0.
LOOP i=1 TO k.
- DO IF expect(i) LT 5.
-  COMPUTE flag=flag+1.
- END IF.
END LOOP.
COMPUTE pflag=100*flag/k.
DO IF flag GT 0.
- PRINT pflag
 /FORMAT='F8.1'
 /TITLE='ВНИМАНИЕ: присутствуют ячейки с ожидаемой частотой, меньшей, чем 5.'
 /RLABEL='% таких ячеек='.
- PRINT minexp
 /FORMAT='F8.1'
 /TITLE='Минимальная ожидаемая частота ячейки:'
 /RLABEL='Ожид.='.
END IF.
END MATRIX.

**********************
* КОНЕЦ СИНТАКСИСА.