Мне нужно применить модель GLM к моим данным, которая очень похожа на пример в пакете lme4 R для метода glmer. Используя их данные cbpp, пример в R выглядит примерно так:
library(lme4)
head(cbpp)
m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
family=binomial, data=cbpp)
summary(m1)
herd incidence size period
1 2 14 1
1 3 12 2
1 4 9 3
1 0 5 4
2 3 22 1
2 1 18 2
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: cbind(incidence, size - incidence) ~ period + (1 | herd)
Data: cbpp
AIC BIC logLik deviance df.resid
194.1 204.2 -92.0 184.1 51
Scaled residuals:
Min 1Q Median 3Q Max
-2.3816 -0.7889 -0.2026 0.5142 2.8791
Random effects:
Groups Name Variance Std.Dev.
herd (Intercept) 0.4123 0.6421
Number of obs: 56, groups: herd, 15
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3983 0.2312 -6.048 1.47e-09 ***
period2 -0.9919 0.3032 -3.272 0.001068 **
period3 -1.1282 0.3228 -3.495 0.000474 ***
period4 -1.5797 0.4220 -3.743 0.000182 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) perid2 perid3
period2 -0.363
period3 -0.340 0.280
period4 -0.260 0.213 0.198
Все мои данные обрабатываются в Python, поэтому вместо того, чтобы записывать данные на диск, создайте Rscript для запустить на данных, а затем снова прочитать результат обратно к Python, думал, возможно ли будет сделать модель GLM в Python, вероятно, с пакетом statsmodels . Однако, пытаясь реализовать это с помощью statsmodels, я получаю результаты, отличные от результатов, полученных с помощью R, поэтому я считаю, что неправильно определяю термины формулы. То, что я пробовал:
import pandas as pd
import statsmodels.api as sm
df = pd.read_csv('cbpp.Rdataset.lme4.tsv')
df['prop'] = df['incidence'] / df['size'] # the incidence proportion
df.head()
formula = 'prop ~ period + (1 | herd)'
mod = sm.GLM.from_formula(formula=formula, data=df,
family=sm.families.Binomial())
res = mod.fit()
print(res.summary())
, что приводит к:
herd incidence size period prop
0 1 2 14 1 0.142857
1 1 3 12 2 0.250000
2 1 4 9 3 0.444444
3 1 0 5 4 0.000000
4 2 3 22 1 0.136364
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: prop No. Observations: 56
Model: GLM Df Residuals: 53
Model Family: Binomial Df Model: 2
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -13.743
Date: Fri, 17 Apr 2020 Deviance: 8.0863
Time: 12:52:50 Pearson chi2: 8.61
No. Iterations: 6
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -0.2391 1.228 -0.195 0.846 -2.645 2.167
period -0.6000 0.442 -1.358 0.174 -1.466 0.266
1 | herd -0.0700 0.101 -0.695 0.487 -0.267 0.127
==============================================================================
Перехваты сильно различаются между двумя методами ...
Как можно указать модель R выше в statsmodels?