gbmt: Estimation of a group-based multivariate trajectory model
In gbmt: Group-Based Multivariate Trajectory Modeling
gbmt R Documentation
Estimation of a group-based multivariate trajectory model
Description
Estimation of a group-based multivariate trajectory model through the Expectation-Maximization (EM) algorithm. Missing values are allowed and the panel may be unbalanced.
Usage
gbmt(x.names, unit, time, ng=1, d=2, data, scaling=2, pruning=TRUE, nstart=NULL,
tol=1e-4, maxit=1000, quiet=FALSE)
Arguments
x.names
Character vector including the names of the indicators.
unit
Character indicating the name of the variable identifying the units.
time
Character indicating the name of the variable identifying the time points.
ng
Positive integer value indicating the number of groups to create. Default is 1.
d
Positive integer value indicating the polynomial degree of group trajectories. Default is 2.
data
Object of class data.frame containing the variables indicated in arguments x.names, unit and time. The variable indicated in argument unit must be of type 'character' or 'factor' and cannot contain missing values. The variable indicated in argument time must be of type 'numeric' or 'Date' and cannot contain missing values. Variables indicated in argument x.names must be of type 'numeric' and may contain missing values. Variables indicated in argument x.names which are completely missing or not present in data will be ignored. The time points may differ across units (unbalanced panel) but they must be unique within units.
scaling
Normalisation method, that should be indicated as: 0 (no normalisation), 1 (centering), 2 (standardization), 3 (ratio to the mean) and 4 (logarithmic ratio to the mean). Default is 2 (standardization). See 'Details'.
pruning
Logical value indicating whether non-significant polynomial terms should be dropped. Default is TRUE. See 'Details'.
nstart
Positive integer value indicating the number of random restarts of the EM algorithm. If NULL (the default), the EM algorithm is started from the solution of a hierarchical cluster with Ward's linkage.
tol
Positive value indicating the tolerance of the EM algorithm. Default is 1e-4.
maxit
Positive integer value indicating the maximum number of iterations of the EM algorithm. Default is 1000.
quiet
Logical value indicating whether prompt messages should be suppressed. Default is FALSE.
Details
Let Y_1,…,Y_k,…,Y_K be the considered indicators and \mbox{y}_{i,t}=(y_{i,t,1},…,y_{i,t,k},…,y_{i,t,K})' denote their observation on unit i (i=1,…,n) at time t (t=1,…,T).
Also, let \bar{y}_{i,k} and s_{i,k} be, respectively, sample mean and sample standard deviation of indicator Y_k for unit i across the whole period of observation.
Each indicator is normalized within units according to one among the following normalisation methods:
0) no normalisation:
y^*_{i,t,k}=y_{i,t,k}
1) centering:
y^*_{i,t,k}=y_{i,t,k}-\bar{y}_{i,k}
2) standardization:
y^*_{i,t,k}=\frac{y_{i,t,k}-\bar{y}_{i,k}}{s_{i,k}}
3) ratio to the mean:
y^*_{i,t,k}=\frac{y_{i,t,k}}{\bar{y}_{i,k}}
4) logarithmic ratio to the mean:
y^*_{i,t,k}=\log≤ft(\frac{y_{i,t,k}}{\bar{y}_{i,k}}\right)\approx\frac{y_{i,t,k}-\bar{y}_{i,k}}{\bar{y}_{i,k}}
Normalisation is required if the trajectories have different levels across units. When indicators have different scales of measurement, standardization is needed to compare the measurements of different indicators. Ratio to the mean and logaritmic ratio to the mean allow comparisons among different indicators as well, but they can be applied only in case of strictly positive measurements.
Denote the hypothesized groups as j=1,…,J and let G_i be a latent variable taking value j if unit i belongs to group j.
A group-based multivariate trajectory model with polynomial degree d is defined as:
\mbox{y}^*_{i,t}\mid G_i=j\sim\mbox{MVN}≤ft(μ_j,Σ_j\right)\hspace{.9cm}j=1,…,J
μ_j=\mbox{B}_j'≤ft(1,t,t^2,…,t^d\right)'
where \mbox{B}_j is the (d+1)\times K matrix of regression coefficients in group j, and Σ_j is the K \times K covariance matrix of the indicators in group j.
The likelihood of the model is:
\mathcal{L}(\mbox{B}_1,…,\mbox{B}_J,Σ_1,…,Σ_J,π_1,…,π_J)=∏_{i=1}^n≤ft[∑_{j=1}^Jπ_j ∏_{t=1}^Tφ(\mbox{y}^*_{i,t}\mid \mbox{B}_j,Σ_j)\right]
where φ(\mbox{y}^*_{i,t}\mid \mbox{B}_j,Σ_j) is the multivariate Normal density of \mbox{y}^*_{i,t} in group j, and π_j is the prior probability of group j.
The posterior probability of group j for unit i is computed as:
\mbox{Pr}(G_i=j\mid \mbox{y}^*_i)\equivπ_{i,j}=\frac{\widehat{π}_j ∏_{t=1}^{T}φ(\mbox{y}^*_{i,t}\mid \widehat{\mbox{B}}_j,\widehat{Σ}_j)}{∑_{j=1}^J\widehat{π}_j ∏_{t=1}^{T}φ(\mbox{y}^*_{i,t}\mid \widehat{\mbox{B}}_j,\widehat{Σ}_j)}
where the hat symbol above a parameter denotes the estimated value for that parameter.
See the vignette of the package and Magrini (2022) for details on maximum likelihood estimation through the EM algorithm.
S3 methods available for class gbmt include:
print: to see the estimated regression coefficients for each group;
summary: to obtain the summary of the linear regressions (a list with one component for each group and each indicator);
plot: to display estimated and predicted trajectories. See plot.gbmt for details;
coef: to see the estimated coefficients (a list with one component for each group);
fitted: to obtain the fitted values, equating to the estimated group trajectories (a list with one component for each group);
residuals: to obtain the residuals (a list with one component for each group);
predict: to perform prediction on trajectories. See predict.gbmt for details.
logLik: to get the log likelihood;
AIC, extractAIC: to get the Akaike information criterion;
BIC: to get the Bayesian information criterion.
Value
An object of class gbmt, including the following components:
call: list including details on the call.
prior: vector including the estimated prior probabilities.
beta: list of matrices, one for each group, including the estimated regression coefficients.
Sigma: list of matrices, one for each group, including the estimated covariance matrix of the indicators.
posterior: matrix including posterior probabilities.
Xmat: the model matrix employed for each indicator in each group.
fitted: list of matrices, one for each group, including the estimated group trajectories.
reg: list of objects of class lm, one for each group and each indicator, including the fitted regressions.
assign: vector indicating the assignement of the units to the groups.
assign.list: list indicating the assignement of the units to the groups.
logLik: log-likelihood of the model.
npar: total number of free parameters in the model.
ic: information criteria for the model (see Magrini, 2022 for details.
appa: average posterior probability of assignments (APPA) for the model.
data.orig: data provided to argument data.
data.scaled: data after normalization.
data.imputed: data after imputation of missing values, equal to data.orig if there are no missing data.
em: matrix with one row for each run of the EM algorithm, including log-likelihood, number of iterations and convergence status (1=yes, 0=no).
References
A. Magrini (2022). Assessment of agricultural sustainability in European Union countries: A group-based multivariate trajectory approach. Advances in Statistical Analysis, published online: March 2022. DOI: 10.1007/s10182-022-00437-9
See Also
plot.gbmt; predict.gbmt.
Examples
data(agrisus2)
# names of indicators (just a subset for illustration)
varNames <- c("TFP_2005", "NetCapital_GVA",
"Income_rur", "Unempl_rur", "GHG_UAA", "GNB_N_UAA")
# model with 2 degrees and 3 groups using the imputed dataset
# - log ratio to the mean is used as normalisation (scaling=4), thus values
# represent relative changes with respect to country averages (see Magrini, 2022)
# - by default, standardization (scaling=2) is used.
m3_2 <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=3, data=agrisus2, scaling=4)
## NOT RUN: same model with multiple random restarts
#m3_2r <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=3, data=agrisus2,
# scaling=4, nstart=10)
# resulting groups
m3_2$assign.list
# estimated group trajectories
m3_2$fitted
# summary of regressions by group
summary(m3_2)
# fit a model with 4 groups
m4_2 <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=4, data=agrisus2,
scaling=4)
rbind(m3_2$ic, m4_2$ic) ## comparison
gbmt documentation built on March 18, 2022, 6:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| gbmt | R Documentation |
Estimation of a group-based multivariate trajectory model
Description
Estimation of a group-based multivariate trajectory model through the Expectation-Maximization (EM) algorithm. Missing values are allowed and the panel may be unbalanced.
Usage
gbmt(x.names, unit, time, ng=1, d=2, data, scaling=2, pruning=TRUE, nstart=NULL, tol=1e-4, maxit=1000, quiet=FALSE)
Arguments
x.names |
Character vector including the names of the indicators. |
unit |
Character indicating the name of the variable identifying the units. |
time |
Character indicating the name of the variable identifying the time points. |
ng |
Positive integer value indicating the number of groups to create. Default is 1. |
d |
Positive integer value indicating the polynomial degree of group trajectories. Default is 2. |
data |
Object of class |
scaling |
Normalisation method, that should be indicated as: 0 (no normalisation), 1 (centering), 2 (standardization), 3 (ratio to the mean) and 4 (logarithmic ratio to the mean). Default is 2 (standardization). See 'Details'. |
pruning |
Logical value indicating whether non-significant polynomial terms should be dropped. Default is |
nstart |
Positive integer value indicating the number of random restarts of the EM algorithm. If |
tol |
Positive value indicating the tolerance of the EM algorithm. Default is 1e-4. |
maxit |
Positive integer value indicating the maximum number of iterations of the EM algorithm. Default is 1000. |
quiet |
Logical value indicating whether prompt messages should be suppressed. Default is |
Details
Let Y_1,…,Y_k,…,Y_K be the considered indicators and \mbox{y}_{i,t}=(y_{i,t,1},…,y_{i,t,k},…,y_{i,t,K})' denote their observation on unit i (i=1,…,n) at time t (t=1,…,T). Also, let \bar{y}_{i,k} and s_{i,k} be, respectively, sample mean and sample standard deviation of indicator Y_k for unit i across the whole period of observation. Each indicator is normalized within units according to one among the following normalisation methods:
0) no normalisation:
y^*_{i,t,k}=y_{i,t,k}
1) centering:
y^*_{i,t,k}=y_{i,t,k}-\bar{y}_{i,k}
2) standardization:
y^*_{i,t,k}=\frac{y_{i,t,k}-\bar{y}_{i,k}}{s_{i,k}}
3) ratio to the mean:
y^*_{i,t,k}=\frac{y_{i,t,k}}{\bar{y}_{i,k}}
4) logarithmic ratio to the mean:
y^*_{i,t,k}=\log≤ft(\frac{y_{i,t,k}}{\bar{y}_{i,k}}\right)\approx\frac{y_{i,t,k}-\bar{y}_{i,k}}{\bar{y}_{i,k}}
Normalisation is required if the trajectories have different levels across units. When indicators have different scales of measurement, standardization is needed to compare the measurements of different indicators. Ratio to the mean and logaritmic ratio to the mean allow comparisons among different indicators as well, but they can be applied only in case of strictly positive measurements.
Denote the hypothesized groups as j=1,…,J and let G_i be a latent variable taking value j if unit i belongs to group j. A group-based multivariate trajectory model with polynomial degree d is defined as:
\mbox{y}^*_{i,t}\mid G_i=j\sim\mbox{MVN}≤ft(μ_j,Σ_j\right)\hspace{.9cm}j=1,…,J
μ_j=\mbox{B}_j'≤ft(1,t,t^2,…,t^d\right)'
where \mbox{B}_j is the (d+1)\times K matrix of regression coefficients in group j, and Σ_j is the K \times K covariance matrix of the indicators in group j. The likelihood of the model is:
\mathcal{L}(\mbox{B}_1,…,\mbox{B}_J,Σ_1,…,Σ_J,π_1,…,π_J)=∏_{i=1}^n≤ft[∑_{j=1}^Jπ_j ∏_{t=1}^Tφ(\mbox{y}^*_{i,t}\mid \mbox{B}_j,Σ_j)\right]
where φ(\mbox{y}^*_{i,t}\mid \mbox{B}_j,Σ_j) is the multivariate Normal density of \mbox{y}^*_{i,t} in group j, and π_j is the prior probability of group j. The posterior probability of group j for unit i is computed as:
\mbox{Pr}(G_i=j\mid \mbox{y}^*_i)\equivπ_{i,j}=\frac{\widehat{π}_j ∏_{t=1}^{T}φ(\mbox{y}^*_{i,t}\mid \widehat{\mbox{B}}_j,\widehat{Σ}_j)}{∑_{j=1}^J\widehat{π}_j ∏_{t=1}^{T}φ(\mbox{y}^*_{i,t}\mid \widehat{\mbox{B}}_j,\widehat{Σ}_j)}
where the hat symbol above a parameter denotes the estimated value for that parameter. See the vignette of the package and Magrini (2022) for details on maximum likelihood estimation through the EM algorithm.
S3 methods available for class gbmt include:
print: to see the estimated regression coefficients for each group;summary: to obtain the summary of the linear regressions (a list with one component for each group and each indicator);plot: to display estimated and predicted trajectories. See plot.gbmt for details;coef: to see the estimated coefficients (a list with one component for each group);fitted: to obtain the fitted values, equating to the estimated group trajectories (a list with one component for each group);residuals: to obtain the residuals (a list with one component for each group);predict: to perform prediction on trajectories. See predict.gbmt for details.logLik: to get the log likelihood;AIC,extractAIC: to get the Akaike information criterion;BIC: to get the Bayesian information criterion.
Value
An object of class gbmt, including the following components:
call: list including details on the call.prior: vector including the estimated prior probabilities.beta: list of matrices, one for each group, including the estimated regression coefficients.Sigma: list of matrices, one for each group, including the estimated covariance matrix of the indicators.posterior: matrix including posterior probabilities.Xmat: the model matrix employed for each indicator in each group.fitted: list of matrices, one for each group, including the estimated group trajectories.reg: list of objects of classlm, one for each group and each indicator, including the fitted regressions.assign: vector indicating the assignement of the units to the groups.assign.list: list indicating the assignement of the units to the groups.logLik: log-likelihood of the model.npar: total number of free parameters in the model.ic: information criteria for the model (see Magrini, 2022 for details.appa: average posterior probability of assignments (APPA) for the model.data.orig: data provided to argumentdata.data.scaled: data after normalization.data.imputed: data after imputation of missing values, equal todata.origif there are no missing data.em: matrix with one row for each run of the EM algorithm, including log-likelihood, number of iterations and convergence status (1=yes, 0=no).
References
A. Magrini (2022). Assessment of agricultural sustainability in European Union countries: A group-based multivariate trajectory approach. Advances in Statistical Analysis, published online: March 2022. DOI: 10.1007/s10182-022-00437-9
See Also
plot.gbmt; predict.gbmt.
Examples
data(agrisus2)
# names of indicators (just a subset for illustration)
varNames <- c("TFP_2005", "NetCapital_GVA",
"Income_rur", "Unempl_rur", "GHG_UAA", "GNB_N_UAA")
# model with 2 degrees and 3 groups using the imputed dataset
# - log ratio to the mean is used as normalisation (scaling=4), thus values
# represent relative changes with respect to country averages (see Magrini, 2022)
# - by default, standardization (scaling=2) is used.
m3_2 <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=3, data=agrisus2, scaling=4)
## NOT RUN: same model with multiple random restarts
#m3_2r <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=3, data=agrisus2,
# scaling=4, nstart=10)
# resulting groups
m3_2$assign.list
# estimated group trajectories
m3_2$fitted
# summary of regressions by group
summary(m3_2)
# fit a model with 4 groups
m4_2 <- gbmt(x.names=varNames, unit="Country", time="Year", d=2, ng=4, data=agrisus2,
scaling=4)
rbind(m3_2$ic, m4_2$ic) ## comparison
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.