<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Area level model
Step1:  EBLUP Predictor
As shown in the milk expenditure datasets, some of the coefficients of variation are not small which indicates unstability of the direct survey estimates. Hence, we can try to reduce the variability of the estimates by smoothing them through modeling. For illustration purpose, we will model the average expenditure on milk using the major areas as auxiliary variables.
First, we use the method fit() to estimate the model parameters. The pandas's method get_dummies() create a matrix with dummy values (0 and 1) from the categorical variable major_area.
Step2:  Now the the model has been fitted, we can obtain the EBLUP average expenditure on milk by running predict() which is a method of EblupAreaModel class. This run will produce two main attributes that is area_est and area_mse which are python dictionaries pairing the small areas to the eblup estimates and the MSE estimates, respectively.
Step3:  We can use the utility method to_dataframe() to output the estimates as a pandas dataframe. The function provides the area, the estimate and its MSE estimates. We can use col_names to customize the name of the columns. For example, using col_names = ["small_area", "eblup_estimate", "eblup_mse"]. Otherwise, if col_names is not provided, "_area", "_estimates" and "_mse" are used as defaults.
Step4:  We could also fit the model parameters using the maximum likelihood (ML) method which will impact the MSE estimation as well. To estimate the area means using the ML methdo, we only need to set method="ML" then run the prediction as follows.
Step5:  Similar, we can use the Fay-Herriot method as follows
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
import samplics 
from samplics.datasets import ExpenditureMilk
from samplics.sae import EblupAreaModel
# Load Expenditure on Milk sample data
milk_exp_cls = ExpenditureMilk()
milk_exp_cls.load_data()
milk_exp = milk_exp_cls.data
nb_obs = 15
print(f"\nFirst {nb_obs} observations of the Milk Expendure dataset\n")
milk_exp.tail(nb_obs)
Explanation: Area level model: empirical best linear unbiased predictor (EBLUP)
Small area estimation (SAE) are useful techniques when the sample sizes are not sufficient to provide reliable direct domain estimates given the sampling design. In this tutorial, the direct estimates refer to estimates obtained from the design-based approach. It usually consists of applying adjusted design weights to the variable of interest to compute sample parameters as estimates of equivalent population parameters. When auxiliary information is available, we can use model assisted survey methods can be used to estimate population parameters. 
In this tutorial, we will go futher and use modeling techniques to produce domains estimates. For the area level model, the modeling is done at the area level using generalized linear mixed models. The sections below shows how to use the EblupAreaModel class from the samplics package to produce area level estimates. 
Milk Expenditure data
To illustrate the EblupAreaModel class, we will use the Milk Expenditure dataset used in Rao and Molina (2015). As mentioned in the book, this dataset was originally used by Arora and Lahiri (1997) and later by You and Chapman (2006). For the R users, this dataset is also used by the R package sae (https://cran.r-project.org/web/packages/sae/index.html). 
The Milk Expenditure data contains 43 observations on the average expenditure on fresh milk for the year 1989. The datasets has the following values: major area representing (major_area), small area (small_area), sample size (samp_size), direct survey estimates of average expenditure (direct_est), standard error of the direct estimate (std_error), and coefficient of variation of the direct estimates (coef_variance).
End of explanation
area = milk_exp["small_area"]
yhat = milk_exp["direct_est"]
X = pd.get_dummies(milk_exp["major_area"],drop_first=True)
sigma_e = milk_exp["std_error"]
## REML method
fh_model_reml = EblupAreaModel(method="REML")
fh_model_reml.fit(
    yhat=yhat, X=X, area=area, error_std=sigma_e, intercept=True, tol=1e-8,
)
print(f"\nThe estimated fixed effects are: {fh_model_reml.fixed_effects}")
print(f"\nThe estimated standard error of the area random effects is: {fh_model_reml.re_std}")
print(f"\nThe convergence statistics are: {fh_model_reml.convergence}")
print(f"\nThe goodness of fit statistics are: {fh_model_reml.goodness}\n")
Explanation: EBLUP Predictor
As shown in the milk expenditure datasets, some of the coefficients of variation are not small which indicates unstability of the direct survey estimates. Hence, we can try to reduce the variability of the estimates by smoothing them through modeling. For illustration purpose, we will model the average expenditure on milk using the major areas as auxiliary variables.
First, we use the method fit() to estimate the model parameters. The pandas's method get_dummies() create a matrix with dummy values (0 and 1) from the categorical variable major_area.
End of explanation
fh_model_reml.predict(
    X=X, area=area, intercept=True
)
import pprint
pprint.pprint(fh_model_reml.area_est)
Explanation: Now the the model has been fitted, we can obtain the EBLUP average expenditure on milk by running predict() which is a method of EblupAreaModel class. This run will produce two main attributes that is area_est and area_mse which are python dictionaries pairing the small areas to the eblup estimates and the MSE estimates, respectively.
End of explanation
milk_est_reml = fh_model_reml.to_dataframe(col_names = ["small_area", "eblup_estimate", "eblup_mse"])
print(f"\nThe dataframe version of the area level estimates:\n\n {milk_est_reml}")
Explanation: We can use the utility method to_dataframe() to output the estimates as a pandas dataframe. The function provides the area, the estimate and its MSE estimates. We can use col_names to customize the name of the columns. For example, using col_names = ["small_area", "eblup_estimate", "eblup_mse"]. Otherwise, if col_names is not provided, "_area", "_estimates" and "_mse" are used as defaults.
End of explanation
## ML method
fh_model_ml = EblupAreaModel(method="ML")
fh_model_ml.fit(
    yhat=yhat, X=X, area=area, error_std=sigma_e, intercept=True, tol=1e-8,
)
milk_est_ml = fh_model_ml.predict(
    X=X, area=area, intercept=True
)
milk_est_ml = fh_model_ml.to_dataframe(col_names = ["small_area", "eblup_estimate", "eblup_mse"])
print(f"\nThe dataframe version of the ML area level estimates:\n\n {milk_est_ml}")
Explanation: We could also fit the model parameters using the maximum likelihood (ML) method which will impact the MSE estimation as well. To estimate the area means using the ML methdo, we only need to set method="ML" then run the prediction as follows.
End of explanation
## FH method
fh_model_fh = EblupAreaModel(method="FH")
fh_model_fh.fit(
    yhat=yhat, X=X, area=area, error_std=sigma_e, intercept=True, tol=1e-8,
)
milk_est_fh = fh_model_fh.predict(
    X=X, area=area, intercept=True
)
milk_est_fh = fh_model_fh.to_dataframe(col_names = ["small_area", "eblup_estimate", "eblup_mse"])
print(f"\nThe dataframe version of the ML area level estimates:\n\n {milk_est_fh}")
Explanation: Similar, we can use the Fay-Herriot method as follows
End of explanation
<END_TASK>