  Marginally specified models have recently become a popular tool for discrete
longitudinal data analysis. Nonetheless, they introduce complex constraint
equations and model fitting algorithms. Moreover, there is a lack of available
software to fit these models. In this paper, we propose a three-level
marginally specified model for analysis of multivariate longitudinal binary
response data. The implicit function theorem is introduced to approximately
solve the marginal constraint equations explicitly. Furthermore, the use of
\textit{probit} link enables direct solutions to the convolution equations. We
propose an R package \textbf{pnmtrem} to fit the model. A simulation study is
conducted to examine the properties of the estimator. We illustrate the model
on the Iowa Youth and Families Project data set.
