Estimation

The estimation of this model using a Poisson likelihood is as follows:

gravity_results <- femlm(Euros ~ log(dist_km), trade, family = "poisson", cluster = c("Origin", "Destination", "Product", "Year"))

Note that alternatively you could have done the same estimation without the argument cluster, by inserting the fixed-effects directly in the formula with a pipe. Further, you can avoid the argument family since the Poisson model is the default:

gravity_results <- femlm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade)

The results can be shown directly with the print method:

print(gravity_results)
#> ML estimation, family = Poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Cluster sizes: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors type: Standard 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5276   0.001925 -793.65 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC:  -768,493.77      Pseudo-R2: 0.74743 
#> Log-likelihood: 1,538,211.80   Squared Cor.: 0.61201

The print reports the coefficient estimates and standard-errors as well as some other information. Among the quality of fit information, the squared-correlation corresponds to the correlation between the dependent variable and the expected predictor; it reflects somehow to the idea of R-square in OLS estimations.

The last line reports the convergence status of the optimization algorithm (the function nlminb from the package stats). In the, rare, event of no convergence, see the troubleshooting section of this Vignette to explore the possible causes and solutions.

Clustering the standard-errors

To cluster the standard-errors, we can simply use the argument se of the summary method. Let’s say we want to cluster the standard-errors according to the first two clusters (i.e. the Origin and Destination variables). Then we just have to do:

summary(gravity_results, se = "twoway")
#> ML estimation, family = Poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Cluster sizes: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors type: Two-way (Origin & Destination) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5276    0.12627 -12.098 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC:  -768,493.77      Pseudo-R2: 0.74743 
#> Log-likelihood: 1,538,211.80   Squared Cor.: 0.61201

The clustering can be done on one (se="cluster"), two (se="twoway"), three (se="threeway") or up to four (se="fourway") variables. If the estimation includes fixed-effects, then by default the clustering will be done using these fixed-effects, in the original order. This is why the Origin and Destination variables were used for the two-way clustering in the previous example. If, instead, you wanted to perform one-way clustering on the Product variable, you need to use the argument cluster:

# Equivalent ways of clustering the SEs:
# - using the vector:
summary(gravity_results, se = "cluster", cluster = trade$Product)
#> ML estimation, family = Poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Cluster sizes: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors type: Clustered 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5276   0.095772 -15.951 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC:  -768,493.77      Pseudo-R2: 0.74743 
#> Log-likelihood: 1,538,211.80   Squared Cor.: 0.61201
# - by reference (only possible because Product has been used as a cluster in the estimation):
summary(gravity_results, se = "cluster", cluster = "Product")
#> ML estimation, family = Poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Cluster sizes: Origin: 15,  Destination: 15,  Product: 20,  Year: 10
#> Standard-errors type: Clustered (Product) 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> log(dist_km)  -1.5276   0.095772 -15.951 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC:  -768,493.77      Pseudo-R2: 0.74743 
#> Log-likelihood: 1,538,211.80   Squared Cor.: 0.61201

Note that you can always cluster the standard-errors, even when the estimation contained no fixed-effect. Buth then you need to use the argument cluster and give explicitely the data to cluster with:

gravity_simple = femlm(Euros ~ log(dist_km), trade)
summary(gravity_simple, se = "twoway", cluster = trade[, c("Origin", "Destination")])
#> ML estimation, family = Poisson, Dep. Var.: Euros
#> Observations: 38,325 
#> Standard-errors type: Two-way 
#>              Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)   10.8940   1.086600 10.0250 < 2.2e-16 ***
#> log(dist_km)  -1.0289   0.152658 -6.7401  1.58e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC: -2,492,169.10      Pseudo-R2: 0.18101 
#> Log-likelihood:  4,984,380.42   Squared Cor.: 0.05511

Using other Likelihoods

To use other likelihood models, you simply have to use the argument family. The fixed-effects Negative Binomial estimation is:

gravity_results_negbin <- femlm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, family = "negbin")

Now to estimate the same relationship with the Gaussian likelihood, we need to put the left hand side in logarithm (since the right-hand-side is not exponentialized):

gravity_results_gaussian <- femlm(log(Euros) ~ log(dist_km)|Origin+Destination+Product+Year, trade, family = "gaussian")

Note that the estimates with the Gaussian likelihood are equivalent to the estimates with Ordinary Least Squares.

Viewing the results in R

Now let’s say that we want a compact overview of the results of several estimations. The best way is to use the function res2table. This function summarizes the results of several femlm estimations into a data.frame. To see the fixed-effects results with the three different likelihoods, we just have to type:

res2table(gravity_results, gravity_results_negbin, gravity_results_gaussian, se = "twoway", subtitles = c("Poisson", "Negative Binomial", "Gaussian"))

We added the argument se="twoway" to cluster the standard-errors for all estimations. As can be seen this function gives an overview of the estimates and standard-errors, as well as some quality of fit measures. The argument subtitles is used to add information on each estimation column.

In the previous example, we directly added the estimation results as arguments of the function res2table. But the function also accepts lists of estimations. Let’s give an example. Say you want to see the influence of the introduction of fixed-effects on the estimate of the elasticity of distance. You can do it with the following code:

gravity_subcluster = list()
all_clusters = c("Year", "Destination", "Origin", "Product")
for(i in 1:4){
    gravity_subcluster[[i]] = femlm(Euros ~ log(dist_km), trade, cluster = all_clusters[1:i])
}

The previous code performs 4 estimations with an increasing number of fixed-effects and store their results into the list named gravity_subcluster. To show the results of all 4 estimations, it’s easy:

res2table(gravity_subcluster, se = "twoway", cluster = trade[, c("Origin", "Destination")])

We have a view of the 5 estimations, all reporting two-way clustered standard-errors thanks to the use of the arguments se and cluster.

Exporting the results to Latex

So far we have seen how to report the results of multiple estimations on the R console. Now, with the function res2tex, we can export the results to high quality Latex tables. The function res2tex works exactly as the function res2table, it takes any number of femlm estimations. By default, it reports Latex code on the R console:

res2tex(gravity_subcluster, se = "twoway", cluster = trade[, c("Origin", "Destination")])
#> \global\long\def\sym#1{\ifmmode^{#1}\else\ensuremath{^{#1}}\fi}
#> \begin{table}[htbp]\centering
#> \caption{no title}
#> \begin{tabular}{lcccc}
#>  & & & & \tabularnewline
#> \hline
#> \hline
#> Dependent Variable:&\multicolumn{4}{c}{Euros}\\
#> Model:&(1)&(2)&(3)&(4)\\
#> \hline
#> \emph{Variables}\tabularnewline
#> log(dist_km)&-1\sym{***}&-1\sym{***}&-2\sym{***}&-2\sym{***}\\
#>   &(0.1527)&(0.1975)&(0.1238)&(0.1263)\\
#> \hline
#> \emph{Fixed-Effects}&  & & & \\
#> Year&Yes&Yes&Yes&Yes\\
#> Destination&No&Yes&Yes&Yes\\
#> Origin&No&No&Yes&Yes\\
#> Product&No&No&No&Yes\\
#> \hline
#> \emph{Fit statistics}&  & & & \\
#> Observations& 38,325&38,325&38,325&38,325\\
#> Adj-pseudo $R^2$ &0.18423&0.35047&0.58023&0.74743\\
#> Log-Likelihood & -2,482,334.64&-1,976,479.31&-1,277,291.27&-768,493.77\\
#> BIC & 4,964,785.39&3,953,222.48&2,554,994.15&1,537,599.67\\
#> \hline
#> \hline
#> \multicolumn{5}{l}{\emph{Two-way standard-errors in parenthesis. Signif Codes: ***: 0.01, **: 0.05, *: 0.1}}\\
#> \end{tabular}
#> \end{table}

This function has many optional arguments. The user can export the Latex table directly into a file (argument file), add a title (arg. title) and a label to the table (arg. label).

The coefficients can be renamed easily (arg. dict), some can be dropped (arg. drop) and they can be easily reordered with regular expressions (arg. order).

The significance codes can easily be changed (arg. signifCode) and all quality of fit information can be customized. Among others, the number of fixed-effect per cluster can also be displayed using the argument showClusterSize.

# we set the dictionary once and for all
myDict = c("log(dist_km)" = "$\\ln (Distance)$", "(Intercept)" = "Constant")
# 1st export: we change the signif code and drop the intercept
res2tex(gravity_subcluster, signifCode = c("a" = 0.01, "b" = 0.05), drop = "Int", dict = myDict, file = "Estimation Table.tex", append = FALSE, title = "First export -- normal Standard-errors")
# 2nd export: clustered S-E + distance as the first coefficient
res2tex(gravity_subcluster, se = "cluster", cluster = trade$Product, order = "dist", dict = myDict, file = "Estimation Table.tex", title = "Second export -- clustered standard-errors (on Product variable)")

Extracting the fixed-effects

To obtain the fixed-effects of the estimation, the function getFE must be performed on the results. This function returns a list containing the fixed-effects for each cluster. The summary method helps to have a quick overview:

fixedEffects <- getFE(gravity_results)
summary(fixedEffects)
#> Fixed-effects coefficients.
#>                         Origin Destination Product   Year
#> Number of fixed-effects     15          15      20     10
#> Number of references         0           1       1      1
#> Mean                      9.53        3.09  0.0129  0.157
#> Variance                  1.63        1.23    1.86 0.0129
#> 
#> COEFFICIENTS:
#>   Origin:    BE    LU    NL    DE    GB                 
#>           9.744 6.414 10.62 10.89 9.934 ... 10 remaining
#> -----
#>   Destination: LU    BE   NL    DE    IE                 
#>                 0 2.696 3.23 4.322 2.588 ... 10 remaining
#> -----
#>   Product: 1     2      3     4      5                 
#>            0 1.414 0.6562 1.449 -1.521 ... 15 remaining
#> -----
#>   Year: 2007    2008     2009    2010  2011                
#>            0 0.06912 0.005226 0.07331 0.163 ... 5 remaining

We can see that the fixed-effects are balanced across clusters. Indeed, apart from the first cluster, only one fixed-effect per cluster needs to be set as reference (i.e. fixed to 0) to avoid collinearity across the fixed-effects of the different clusters. This ensures that the fixed-effects coefficients can be compared within cluster. Had there be strictly more than one reference per cluster, their interpretation would have not been possible at all. If this was the case though, a warning message would have been prompted. Note that the mean values are meaningless per se, but give a reference points to which compare the fixed-effects within a cluster. Let’s look specifically at the Year fixed-effects:

fixedEffects$Year
#>        2007        2008        2009        2010        2011        2012 
#> 0.000000000 0.069122369 0.005225534 0.073308290 0.163013422 0.192605196 
#>        2013        2014        2015        2016 
#> 0.230629312 0.242605390 0.282800746 0.310325679

Finally, the plot method helps to distinguish the most notable fixed-effects:

plot(fixedEffects)

For each cluster, the fixed-effects are first centered, then sorted, and finally the most notable (i.e. highest and lowest) are reported. The exponential of the coefficient is reported in the right hand side to simplify the interpretation for models with log-link (as the Poisson model). As we can see from the country of destination cluster, trade involving France (FR), Italy (IT) and Germany (DE) as destination countries is more than 2.7 times higher than the EU15 average. Further, the highest heterogeneity come from the product category, where trade in product 4 (dairy products) is roughly 2.7 times the average while product 14 (vegetable plaiting materials) represents a negligible fraction of the average.

Note however that the interpretation of the fixed-effects must be taken with extra care. In particular, here the fixed-effects can be interpreted only because they are perfectly balanced.

Non-linear in parameters example

The function femlm also allows to have non-linear in parameters right-hand-sides (RHS). First an example without fixed-effects, the one with fixed-effects is given later. Let’s say we want to estimate the following relation with a Poisson model:

E(z_i) = a × x_i + b × y_i.

In fact, this type of model is non-linear in the context of a Poisson model because the sum is embedded within the log:

E(z_i) = exp (log (a × x_i + b × y_i)).

So let’s estimate such a relation. First we generate the data:

# Generating data:
n = 1000
# x and y: two positive random variables
x = rnorm(n, 1, 5)**2
y = rnorm(n, -1, 5)**2
# E(z) = 2*x + 3*y and some noise
z = rpois(n, 2*x + 3*y) + rpois(n, 1)
base = data.frame(x, y, z)

To estimate the non-linear relationship, we need to use the argument NL.fml where we put the non-linear part. We also have to provide starting values with the argument NL.start. Finally, to ensure the RHS can be evaluated in any situation, we add lower bounds for the parameters with the argument lower.

result_NL = femlm(z~0, base, NL.fml = ~ log(a*x + b*y), NL.start = list(a=1, b=1), lower = list(a=0, b=0))

Note that the arguments NL.start and lower are named lists. Setting lower = list(a=0, b=0) means that the optimization algorithm will never explores parameters for a and b that are lower than 0. The results obtained can be interpreted similarly to results with linear RHS. We can see them with a print:

print(result_NL)
#> Non-linear ML estimation, family = Poisson, Dep. Var.: z
#> Observations: 1,000 
#> Standard-errors type: Standard 
#>   Estimate Std. Error z value  Pr(>|z|)    
#> a   2.0384   0.011144  182.92 < 2.2e-16 ***
#> b   2.9921   0.012829  233.23 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>            BIC: -3,595.77      Pseudo-R2: 0.93459 
#> Log-likelihood:  7,219.17   Squared Cor.: 0.99181

We can see that we obtain coefficients close to the generating values.

# the class of each observation
id = sample(20, n, replace = TRUE)
base$id = id
# the vector of fixed-effects
gamma = rnorm(20)**2
# the new vector z_bis
z_bis = rpois(n, gamma[id] * (2*x + 3*y)) + rpois(n, 1)
base$z_bis = z_bis

# we add the fixed-effect in the formula
result_NL_fe = femlm(z_bis~0|id, base, NL.fml = ~ log(2*x + b*y), NL.start = list(b=1), lower = list(b=0))
# The coef should be around 3
coef(result_NL_fe)
#>        b 
#> 2.978882
# the gamma and the exponential of the fixed-effects should be similar
rbind(gamma, exp(getFE(result_NL_fe)$id))
#>               1           2        3           4         5         6          7
#> gamma 0.1415361 0.004566642 2.173635 0.004215653 0.4888401 0.7229866 0.02089668
#>       0.1580355 0.016142141 2.192493 0.010182497 0.5047119 0.7328707 0.03136660
#>               8        9          10       11       12       13       14
#> gamma 0.3215559 1.033220 0.009370039 1.599764 1.873162 2.858219 3.559306
#>       0.3346130 1.030818 0.018401000 1.603864 1.923171 2.871548 3.595552
#>              15        16       17       18        19       20
#> gamma 0.1656070 0.9312762 5.197757 4.236893 0.1572261 1.505646
#>       0.1601142 0.9407736 5.249692 4.255693 0.1689799 1.511539

Using multi-cores

The package FENmlm integrates multi-platform parallelism to hasten the estimation process. To use the multi-core facility, just use the argument cores:

# Sample of results:
# 1 core: 2.7s
system.time(femlm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, family = "negbin", cores = 1))
# 2 cores: 1.74s
system.time(femlm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, family = "negbin", cores = 2))
# 4 cores: 1.31s
system.time(femlm(Euros ~ log(dist_km)|Origin+Destination+Product+Year, trade, family = "negbin", cores = 4))

As you can see, the efficiency of increasing the number of cores is not 1 to 1. Two cores do not divide the computing time by 2, nor four cores by 4. However it still reduces significantly the computing time, which might be valuable for large sample estimations.

Although parallelism has been implemented for all likelihood models, the only ones for which it is really efficient are the Negative Binomial and the Logit models.

Indeed, the underlying femlm code alternate between R and C and the parallelism has been implemented with the C library OpenMP. Setting the parallelism with OpenMP incurs “fixed-costs”, and since there is many back and forth between R and C, these fixed-costs add up significantly. Thus to be worth, the C-code sections with parallelism must be “complex”.

Because the Negative Binomial and Logit models are the ones with the most computationnaly intensive C-code chunks, they are the ones that benefit the most from the parallelism.

Troubleshooting

This section deals with possible error/warning messages stemming from the estimation. There are two identified – possible – problems: collinearity and precision.

base_coll = trade
base_coll$constant_variable = 1
res <- femlm(Euros ~ log(dist_km) + constant_variable|Origin+Destination+Product+Year, base_coll)
#> Warning in femlm(Euros ~ log(dist_km) + constant_variable | Origin +
#> Destination + : [femlm]: The optimization algorithm did not converge, the
#> results are not reliable. The information matrix is singular (likely presence
#> of collinearity). Use function diagnostic() to pinpoint collinearity problems.
diagnostic(res)
#> [1] "Variable 'constant_variable' is collinear with cluster Origin."