Description Usage Arguments Details Value References Examples
Maximum likelihood estimation of the parameters of a correlated probit model via EM algorithm. The function works with wide format of the response data. The function allows for NA values in the outcome.
1 2 3 | emcorrprobit(model, y, xfixed, xrand, start.values.beta,
start.values.delta = NULL, start.values.sigma.rand, exact, montecarlo,
epsilon, ...)
|
model |
specifies the response model; |
y |
2-way array for the the response variable with dimensions: individuals and multiple observations. The ordinal data should be represented by numeric values in the following way: the first level is denoted by the number 1, second by the number 2 and so on. |
xfixed |
3-way array for the predictors for the fixed effects with dimensions: individuals, dimension of the fixed effects and multiple observations. The intercept should be included as well. |
xrand |
3-way array for the predictors for the random effects with dimensions: individuals, dimension of the random effects and multiple observations. |
start.values.beta |
start values for the regression parameters β |
start.values.delta |
start values for the differences in the consecutive thresholds δ. Default NULL for binary data. Otherwise it should be specified. |
start.values.sigma.rand |
a matrix with the start values for the covariance matrix of the random effects Σ |
exact |
logical; if TRUE analytical calculation of the moments of
truncated normal distribution is obtained ( |
montecarlo |
numeric; the number of generated values used for the estimation of the first two moments of truncated normal distribution. If exact=TRUE this parameter is not needed. |
epsilon |
a value for the stopping criterion. |
... |
not used. |
The function fits the latent class probit model:
y_{ij} = x'_{ij}β+ z'_{ij}b_i+ε_{ij},
where y*_{ij} = k is observed, if y_{ij} < α_k and y*_{ij} = m, if y_{ij} > α_{m-1}, the response variable y*_ij may take a value from 1 to m. We assume b_i ~ N(0,Σ) and ε_{ij} ~ N(0,1).
The model is fitted using re-parametrisation where new parameters are defined as: δ_k=α_k-α_{k-1}, k=2,...,m-1.
The stopping criterion of the algorithm is when the differences between the
estimates from two successive iterations of the algorithm are less than
epsilon
for each parameter.
One should choose carefully the starting values for the parameters (especially
for the covariance matrix of the random effects) and the value of epsilon
.
It is possible that the algorithm stops before convergence and over- or
underestimate the parameters. We recommend using different starting values for
the parameters and only after getting similar results, it can be assumed that
obtained estimates are the MLEs.
When the data consists of 2 or 3 observations per subject it is recommended using the
analytical calculation of the moments of truncated normal distribution
(exact=TRUE
).
An object of class emcorrprobit
. List with following components:
Sigma.rand.effects |
The estimated covariance matrix of the random effects Σ. |
regression.coefficients |
The estimated regression coefficients β. |
differences.in.thresholds |
The estimated differences in the consecutive thresholds δ. |
thresholds |
Estimated thresholds α. By definition the first threshold α_1 is zero. |
random.effects |
The estimated random effects for each individual b_i. |
loglikelihood |
Log-likelihood of the model. |
AIC |
Akaike information criterion. |
BIC |
Bayesian information criterion. |
number.iterations |
The number of iterations. |
R. V. Gueorguieva. Correlated probit model. In Encyclopedia of Biopharmaceutical Statistics, chapter 59, pages 355-362. 2006. doi: 10.3109/9781439822463.057. URL http://informahealthcare.com/doi/abs/10.3109/9781439822463.057
D. Grigorova and R. Gueorguieva. Implementation of the EM algorithm for maximum likelihood estimation of a random effects model for one longitudinal ordinal outcome. Pliska Stud. Math. Bulgar., 22:41-56, 2013
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 | ### data simulation
############################################################
### Random intercept model for 3-level ordinal variable ####
### 750 individuals with 2 observations per subject ########
### Predcitors - intercept and time ########################
############################################################
random.int=rnorm(750,0,0.1)
l=length(random.int)
int=-0.5
mult.obs=2
y1=sapply(1:mult.obs,function(i) random.int+int+i+rnorm(l), simplify="array")
data.ordinal=ifelse(y1<=0,1,ifelse(y1<=1.5,2,3))
table(data.ordinal)
head(data.ordinal)
time=sapply(1:mult.obs, function(i) rep(i,l), simplify="array")
predictors.fixed=sapply(1:mult.obs, function(i) cbind(1,time[,i]), simplify="array")
predictors.random=sapply(1:mult.obs, function(i) matrix(rep(1,l),ncol=1), simplify="array")
sigma.rand=matrix(.01)
beta=c(-0.55,0.95)
delta=c(1.5)
mc=500
e=TRUE
### estimation
###should work
example1=emcorrprobit(model = "oneord", y=data.ordinal,xfixed=predictors.fixed,
xrand=predictors.random,
start.values.beta=beta,start.values.delta=delta,
start.values.sigma.rand=sigma.rand,
exact=e,montecarlo=mc,epsilon=.0002)
###doesn't work
#example2=emcorrprobit(model = "1ord", y=data.ordinal,xfixed=predictors.fixed,
# xrand=predictors.random,
# start.values.beta=beta,start.values.delta=delta,
# start.values.sigma.rand=sigma.rand,
# exact=e,montecarlo=mc,epsilon=.0002)
#example2=emcorrprobit(y=data.ordinal,xfixed=predictors.fixed,
# xrand=predictors.random,
# start.values.beta=beta,start.values.delta=delta,
# start.values.sigma.rand=sigma.rand,
# exact=e,montecarlo=mc,epsilon=.0002)
###to see the estimates
example1
## the same as
print(example1)
### Monte Carlo approach for estimation of moments of truncated normal distribution - slower in this
### case
example2=emcorrprobit(model = "oneord", y=data.ordinal,xfixed=predictors.fixed,
xrand=predictors.random,
start.values.beta=beta,start.values.delta=delta,
start.values.sigma.rand=sigma.rand,
exact=FALSE,montecarlo=mc,epsilon=.0002)
### example with missing data
data.ordinal[1,2]=NA
head(data.ordinal)
example3=emcorrprobit(model = "oneord", y=data.ordinal,xfixed=predictors.fixed,
xrand=predictors.random,
start.values.beta=beta,start.values.delta=delta,
start.values.sigma.rand=sigma.rand,
exact=TRUE,montecarlo=mc,epsilon=.0002)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.