femlm: Fixed effects maximum likelihood models

In FENmlm: Fixed Effects Nonlinear Maximum Likelihood Models

Fixed effects maximum likelihood models

Description

This function estimates maximum likelihood models (e.g., Poisson or Logit) and is efficient to handle any number of fixed effects (i.e. cluster variables). It further allows for nonlinear in parameters right hand sides.

Usage

femlm(fml, data, family = c("poisson", "negbin", "logit", "gaussian"),
  NL.fml, cluster, na.rm = FALSE, useAcc = TRUE, NL.start, lower,
  upper, env, NL.start.init, offset, nl.gradient, linear.start = 0,
  jacobian.method = c("simple", "Richardson"), useHessian = TRUE,
  opt.control = list(), cores = 1, verbose = 0, theta.init,
  precision.cluster, itermax.cluster = 10000, itermax.deriv = 5000,
  showWarning = TRUE, ...)

Arguments

fml	A formula. This formula gives the linear formula to be estimated (it is similar to a lm formula), for example: fml = z~x+y. To include cluster variables, you can 1) either insert them in this formula using a pipe (e.g. fml = z~x+y|cluster1+cluster2), or 2) either use the argment cluster. To include a non-linear in parameters element, you must use the argment NL.fml.
data	A data.frame containing the necessary variables to run the model. The variables of the non-linear right hand side of the formula are identified with this data.frame names. Note that no NA is allowed in the variables to be used in the estimation. Can also be a matrix.
family	Character scalar. It should provide the family. The possible values are "poisson" (Poisson model with log-link, the default), "negbin" (Negative Binomial model with log-link), "logit" (LOGIT model with log-link), "gaussian" (Gaussian model).
NL.fml	A formula. If provided, this formula represents the non-linear part of the right hand side (RHS). Note that contrary to the fml argument, the coefficients must explicitely appear in this formula. For instance, it can be ~a*log(b*x + c*x^3), where a, b, and c are the coefficients to be estimated. Note that only the RHS of the formula is to be provided, and NOT the left hand side.
cluster	Character vector. The name/s of a/some variable/s within the dataset to be used as clusters. These variables should contain the identifier of each observation (e.g., think of it as a panel identifier).
na.rm	Logical, default is FALSE. If the variables necessary for the estimation contain NAs and na.rm = TRUE, then all observations containing NA are removed prior to estimation and a warning message is raised detailing the number of observations removed.
useAcc	Default is TRUE. Whether an acceleration algorithm (Irons and Tuck iterations) should be used to otbain the cluster coefficients when there are two or more clusters.
NL.start	(For NL models only) A list of starting values for the non-linear parameters. ALL the parameters are to be named and given a staring value. Example: NL.start=list(a=1,b=5,c=0). Though, there is an exception: if all parameters are to be given the same starting value, you can use the argument NL.start.init.
lower	(For NL models only) A list. The lower bound for each of the non-linear parameters that requires one. Example: lower=list(b=0,c=0). Beware, if the estimated parameter is at his lower bound, then asymptotic theory cannot be applied and the standard-error of the parameter cannot be estimated because the gradient will not be null. In other words, when at its upper/lower bound, the parameter is considered as 'fixed'.
upper	(For NL models only) A list. The upper bound for each of the non-linear parameters that requires one. Example: upper=list(a=10,c=50). Beware, if the estimated parameter is at his upper bound, then asymptotic theory cannot be applied and the standard-error of the parameter cannot be estimated because the gradient will not be null. In other words, when at its upper/lower bound, the parameter is considered as 'fixed'.
env	(For NL models only) An environment. You can provide an environement in which the non-linear part will be evaluated. (May be useful for some particular non-linear functions.)
NL.start.init	(For NL models only) Numeric scalar. If the argument NL.start is not provided, or only partially filled (i.e. there remain non-linear parameters with no starting value), then the starting value of all remaining non-linear parameters is set to NL.start.init.
offset	A formula. An offset can be added to the estimation. It should be a formula of the form (for example) ~0.5*x**2. This offset is linearily added to the elements of the main formula 'fml'. Note that when using the argument 'NL.fml', you can directly add the offset there.
nl.gradient	(For NL models only) A formula. The user can prodide a function that computes the gradient of the non-linear part. The formula should be of the form ~f0(a1,x1,a2,a2). The important point is that it should be able to be evaluated by: eval(nl.gradient[[2]], env) where env is the working environment of the algorithm (which contains all variables and parameters). The function should return a list or a data.frame whose names are the non-linear parameters.
linear.start	Numeric named vector. The starting values of the linear part.
jacobian.method	Character scalar. Provides the method used to numerically compute the jacobian of the non-linear part. Can be either "simple" or "Richardson". Default is "simple". See the help of jacobian for more information.
useHessian	Logical. Should the Hessian be computed in the optimization stage? Default is TRUE.
opt.control	List of elements to be passed to the optimization method nlminb. See the help page of nlminb for more information.
cores	Integer, default is 1. Number of threads to be used (accelerates the algorithm via the use of openMP routines). This is particularly efficient for the negative binomial and logit models, less so for the Gaussian and Poisson likelihoods (unless for very large datasets).
verbose	Integer, default is 0. It represents the level of information that should be reported during the optimisation process. If verbose=0: nothing is reported. If verbose=1: the value of the coefficients and the likelihood are reported. If verbose=2: 1 + information on the computing tiime of the null model, the cluster coefficients and the hessian are reported.
theta.init	Positive numeric scalar. The starting value of the dispersion parameter if family="negbin". By default, the algorithm uses as a starting value the theta obtained from the model with only the intercept.
precision.cluster	Precision used to obtain the fixed-effects (ie cluster coefficients). Defaults to 1e-5. It corresponds to the maximum absolute difference allowed between two iterations. Argument precision.cluster cannot be lower than 10000*.Machine$double.eps.
itermax.cluster	Maximum number of iterations in the step obtaining the fixed-effects (only in use for 2+ clusters). Default is 10000.
itermax.deriv	Maximum number of iterations in the step obtaining the derivative of the fixed-effects (only in use for 2+ clusters). Default is 5000.
showWarning	Logical, default is TRUE. Whether warnings should be displayed (concerns warnings relating to: convergence state, collinearity issues and observation removal due to only 0/1 outcomes or presence of NA values).
...	Not currently used.

Details

This function estimates maximum likelihood models where the conditional expectations are as follows:

Gaussian likelihood:

E(Y|X)=X\beta

Poisson and Negative Binomial likelihoods:

E(Y|X)=\exp(X\beta)

where in the Negative Binomial there is the parameter \theta used to model the variance as \mu+\mu^2/\theta, with \mu the conditional expectation. Logit likelihood:

E(Y|X)=\frac{\exp(X\beta)}{1+\exp(X\beta)}

When there are one or more clusters, the conditional expectation can be written as:

E(Y|X) = h(X\beta+\sum_{k}\sum_{m}\gamma_{m}^{k}\times C_{im}^{k}),

where h(.) is the function corresponding to the likelihood function as shown before. C^k is the matrix associated to cluster k such that C^k_{im} is equal to 1 if observation i is of category m in cluster k and 0 otherwise.

When there are non linear in parameters functions, we can schematically split the set of regressors in two:

f(X,\beta)=X^1\beta^1 + g(X^2,\beta^2)

with first a linear term and then a non linear part expressed by the function g. That is, we add a non-linear term to the linear terms (which are X*beta and the cluster coefficients). It is always better (more efficient) to put into the argument NL.fml only the non-linear in parameter terms, and add all linear terms in the fml argument.

To estimate only a non-linear formula without even the intercept, you must exclude the intercept from the linear formula by using, e.g., fml = z~0.

The over-dispersion parameter of the Negative Binomial family, theta, is capped at 10,000. If theta reaches this high value, it means that there is no overdispersion.

Value

An femlm object.

coefficients	The named vector of coefficients.
coeftable	The table of the coefficients with their standard errors, z-values and p-values.
loglik	The loglikelihood.
iterations	Number of iterations of the algorithm.
n	The number of observations.
nparams	The number of parameters of the model.
call	The call.
fml	The linear formula of the call.
ll_null	Log-likelihood of the null model (i.e. with the intercept only).
pseudo_r2	The adjusted pseudo R2.
message	The convergence message from the optimization procedures.
sq.cor	Squared correlation between the dependent variable and the expected predictor (i.e. fitted.values) obtained by the estimation.
hessian	The Hessian of the parameters.
fitted.values	The fitted values are the expected value of the dependent variable for the fitted model: that is E(Y|X).
cov.unscaled	The variance-covariance matrix of the parameters.
se	The standard-error of the parameters.
scores	The matrix of the scores (first derivative for each observation).
family	The ML family that was used for the estimation.
residuals	The difference between the dependent variable and the expected predictor.
sumFE	The sum of the fixed-effects for each observation.
offset	The offset formula.
NL.fml	The nonlinear formula of the call.
bounds	Whether the coefficients were upper or lower bounded. – This can only be the case when a non-linear formula is included and the arguments 'lower' or 'upper' are provided.
isBounded	The logical vector that gives for each coefficient whether it was bounded or not. This can only be the case when a non-linear formula is included and the arguments 'lower' or 'upper' are provided.
clusterNames	The names of each cluster.
id_dummies	The list (of length the number of clusters) of the cluster identifiers for each observation.
clusterSize	The size of each cluster.
obsRemoved	In the case there were clusters and some observations were removed because of only 0/1 outcome within a cluster, it gives the row numbers of the observations that were removed.
clusterRemoved	In the case there were clusters and some observations were removed because of only 0/1 outcome within a cluster, it gives the list (for each cluster) of the clustr identifiers that were removed.
theta	In the case of a negative binomial estimation: the overdispersion parameter.

Author(s)

Laurent Berge

References

Berge, Laurent, 2018, "Efficient estimation of maximum likelihood models with multiple fixed-effects: the R package FENmlm." CREA Discussion Papers, 13 (https://github.com/lrberge/fixest/blob/master/_DOCS/FENmlm_paper.pdf).

For models with multiple fixed-effects:

Gaure, Simen, 2013, "OLS with multiple high dimensional category variables", Computational Statistics & Data Analysis 66 pp. 8–18

On the unconditionnal Negative Binomial model:

Allison, Paul D and Waterman, Richard P, 2002, "Fixed-Effects Negative Binomial Regression Models", Sociological Methodology 32(1) pp. 247–265

See Also

See also summary.femlm to see the results with the appropriate standard-errors, getFE to extract the cluster coefficients, and the functions res2table and res2tex to visualize the results of multiple estimations.

Examples

#
# Linear examples
#

# Load trade data
data(trade)

# We estimate the effect of distance on trade => we account for 3 cluster effects
# 1) Poisson estimation
est_pois = femlm(Euros ~ log(dist_km)|Origin+Destination+Product, trade)
# alternative formulation giving the same results:
# est_pois = femlm(Euros ~ log(dist_km), trade, cluster = c("Origin", "Destination", "Product"))

# 2) Log-Log Gaussian estimation (with same clusters)
est_gaus = update(est_pois, log(Euros+1) ~ ., family="gaussian")

# 3) Negative Binomial estimation
est_nb = update(est_pois, family="negbin")

# Comparison of the results using the function res2table
res2table(est_pois, est_gaus, est_nb)
# Now using two way clustered standard-errors
res2table(est_pois, est_gaus, est_nb, se = "twoway")

# Comparing different types of standard errors
sum_white = summary(est_pois, se = "white")
sum_oneway = summary(est_pois, se = "cluster")
sum_twoway = summary(est_pois, se = "twoway")
sum_threeway = summary(est_pois, se = "threeway")

res2table(sum_white, sum_oneway, sum_twoway, sum_threeway)


#
# Example of Equivalences
#
## Not run: 
# equivalence with glm poisson
est_glm <- glm(Euros ~ log(dist_km) + factor(Origin) +
            factor(Destination) + factor(Product), trade, family = poisson)

# coefficient estimates + Standard-error
summary(est_glm)$coefficients["log(dist_km)", ]
est_pois$coeftable

# equivalence with lm
est_lm <- lm(log(Euros+1) ~ log(dist_km) + factor(Origin) +
            factor(Destination) + factor(Product), trade)

# coefficient estimates + Standard-error
summary(est_lm)$coefficients["log(dist_km)", ]
summary(est_gaus, dof_correction = TRUE)$coeftable

## End(Not run)


#
# Non-linear examples
#

# Generating data for a simple example
n = 100
x = rnorm(n, 1, 5)**2
y = rnorm(n, -1, 5)**2
z1 = rpois(n, x*y) + rpois(n, 2)
base = data.frame(x, y, z1)

# Estimating a 'linear' relation:
est1_L = femlm(z1 ~ log(x) + log(y), base)
# Estimating the same 'linear' relation using a 'non-linear' call
est1_NL = femlm(z1 ~ 1, base, NL.fml = ~a*log(x)+b*log(y), NL.start = list(a=0, b=0))
# we compare the estimates with the function res2table (they are identical)
res2table(est1_L, est1_NL)

# Now generating a non-linear relation (E(z2) = x + y + 1):
z2 = rpois(n, x + y) + rpois(n, 1)
base$z2 = z2

# Estimation using this non-linear form
est2_NL = femlm(z2~0, base, NL.fml = ~log(a*x + b*y),
               NL.start = list(a=1, b=2), lower = list(a=0, b=0))
# we can't estimate this relation linearily
# => closest we can do:
est2_L = femlm(z2~log(x)+log(y), base)

# Difference between the two models:
res2table(est2_L, est2_NL)

# Plotting the fits:
plot(x, z2, pch = 18)
points(x, fitted(est2_L), col = 2, pch = 1)
points(x, fitted(est2_NL), col = 4, pch = 2)


# Using a custom Jacobian for the function log(a*x + b*y)
myGrad = function(a,x,b,y){
	s = a*x+b*y
	data.frame(a = x/s, b = y/s)
}

est2_NL_grad = femlm(z2~0, base, NL.fml = ~log(a*x + b*y),
                     NL.start = list(a=1,b=2), nl.gradient = ~myGrad(a,x,b,y))

FENmlm documentation built on Aug. 22, 2023, 5:11 p.m.