using MixedLM as lme REML in R

[fortozs]

I am trying to convert some R code over to Python. I have a random effects model in R:

data.model = lme(log(Prey_Num) ~ distance + factor(species), random = ~1 + distance|month,
           control=lmeControl(opt='optim'),data=mydata, method='REML', na.action=na.exclude)

This gives me an AIC value of 327.25. I would like to implement this in Python. I tried following some examples here. So far I have:

model = sm.MixedLM.from_formula("np.log(Prey_Num) ~ distance + C(species)", model_data.dropna(),
                                re_formula="distance", groups=model_data.dropna()['month'])

result = model.fit()

print(result.summary())

aic = 2 * (len(result.params) + 1) - 2 * result.llf

print(aic)

I end up with an AIC of 271.88. In the examples it looks like I should be getting the same values. Can anyone tell me if and how it would be possible to translate this R code to Python? Thanks.

[fortozs]

I managed to get it. I hadn't noticed that the number of observations were different between the two analyses. I guess lme and statsmodels handle missing data differently. Now the AIC and likelihood are identical to the hundredths place. The code I used to make it work is:

model = sm.MixedLM.from_formula("np.log(Prey_Num) ~ distance + C(species)",

                               model_data[model_data['Prey_Num'].notnull()],

                               re_formula = "distance",

                               groups=model_data[model_data['Prey_Num'].notnull()]['month'])

result = model.fit()

print(result.summary())

aic = 2 * (len(result.params) + 1) - 2 * result.llf

print(aic)

Time to model!