In biomechanics, rotation matrices, axes of rotation, and the center of rotation of a moving body segment or a joint are estimated using conventional ordinary least squares (Kim, 2011, 2013). One can approach the least squares method with a general linear model (GLM). The response variable in this model is position vectors of the markers attached to the points of interest of human body obtained from high-speed camera image. It is a longitudinal data measured continuously over time. The values measured repeatedly over time are not independent of each other, and the variance of the repeated measurements is not usually constant over time (Fitzmaurice, Laird, & Ware, 2004; Hedeker & Gibbons, 2006). A series of position vectors measured over time in the same marker during exercise may have positive correlations. Furthermore, as a joint moves, the variance of position vectors may be changed by soft tissue artefact (STA) (Cereatti, Croce, & Cappozzo, 2006), which is an error due to movement of the soft tissue such as skin.

In the general linear model, the mean of the marker position vectors is expressed as a function of the parameter like the center of rotation of a joint (the joint center), and the variance of position vectors, variance-covariance matrix (covariance matrix), comprise only one parameter assuming independence and homoscedasticity for convenience of estimation (Kim, 2017). The more the data, the more the covariance matrix parameters to be estimated, and in some circumstances, this assumption may be valid. This may not be a problem if it is limited to estimating the center of rotation. However, it may become a problem if the interest of the research is related to the variance of the center of the rotation estimator, the correlation between the position vectors, and the interference about the center of rotation. If the appropriate variance and covariance of the measured values are not taken into account, the efficiency of the estimated parameters associated with the mean is reduced (Diggle et al., 2002; Fitzmaurice et al., 2004) and the validity of the estimated variance of the parameters decreases (Cnaan, Laird, & Slasor, 1997). The transformation method used in biomechanics (Ehrig, Taylor, Duda, & Heller, 2006; Piazza, Erdemir, Okita, & Cavanagh, 2004; Siston & Delp, 2006) can estimate the joint center based on the general linear model (Kim, 2017). Although the values are repeatedly measured over time, they still focus on the modeling of the means assuming independence and homoscedasticity. To estimate an efficient joint center, attention should be paid to modeling of variance.

Given the repeatedly measured values of the markers, the experimental unit, an easy way to model variance is to approach it through marginal analysis (Diggle, Heagerty, Liang, & Zeger, 2002). It can model the mean and variance separately, and it is distinguished from the mixed-effects model (Fitzmaurice et al., 2004; Hedeker & Gibbons, 2006). The covariance pattern model is proposed by Jennrich and Schluchter (1986), and it is distinguished from the unstructured covariance matrix, as it models the covariance matrices of various patterns to explain the correlation and heteroscedasticity of the errors of the general linear model by using a small number of parameters. The covariance pattern model is mainly used in time series analysis (Fitzmaurice et al., 2004; Ware, 1985).

The purpose of this study was to estimate efficient joint centers by applying a general linear model to the position vector of the marker, which was continuously measured over time in an image analysis. In this model, the joint center is difined as the parameter related to the mean, and the variance is the covariance pattern model considering the correlation between the continuously measured vectors and heteroscedasticity. By applying various covariance patterns to the covariance matrix of the marker position vectors, the appropriate models were identified, and the joint center vectors were deduced on the basis of the models.

1. Model

As shown in Figure 1, the three-dimensional position vector of arbitrary i marker (i=1,⋯,m), which is repeatedly measured (j=1,⋯,n) in a frame of images continuously taken by a high-speed camera during human body motion can be expressed as follows:

(1)

Figure 1. A hip joint model. Markers on the thigh (distal segment) are observed with respect to the pelvic (proximal segment) coordinate system.

Modify Delete

where ο (omicron) is the center of rotation (CoR) vector of the joint, is the position vector from the joint center to the marker, and is the error vector. , ο, and are represented by the proximal segment coordinate system , and is expressed by the distal segment coordinate system . is a 3×3 transformation matrix that transforms the distal segment component into the proximal segment component. and are continuously measured and vary in n frames. It should be noted that is not , which includes the subscript i of the marker, as equally acts as a transformation matrix for all markers. Both the proximal and distal coordinate systems are rectangular coordinate systems that are fixed to the segments and have the same motion as the segments. Thus, the distal segment coordinate system serves as the transformation matrix . Assuming the body is a rigid body, ο and are constant vectors.

Equation 1 can be expressed by the equation used in the general linear model.

(2)

is a response variable, and is an error. Consecutive n position vectors and error vectors of arbitrary i marker are arrayed in columns, respectively.

(3)

In addition, and are

(4)

where is the covariate, and is the regression parameter. is a vector with the value of 1 for every n components, is an identity matrix of size 3, and ⊗ is a Kronecker product. The component of is a matrix that arrays the transformation matrix in a column, as in .

(5)

It should be noted and in Equation 2. Since matrix consists of A, it is independent of the subscript i, like A.

(6)

This is referred to as a completely balanced design (Ware, 1985), and i is not needed, but considering entire markers below, the model matrix X will be different from , so i is remained for distinction. On the other hand, in terms of , i is necessary because the position of the component differs from marker to marker. In addition, β in the model considering entire markers varies in size depending on the number of markers used. In other words, if all the markers (m) are used, the size of β is . The position of a marker can be a vector of the sum of the joint center ο and its own position . The error is assumed to come from a normal distribution with E()=0 and Var()=. Therefore, the mean and the variance of are as follows:

(7)

The mean vector is a function of and , and the covariance matrix is a conditional variance () when is given. In other words, it is the covariance matrix given the joint center ο and the transformed , which is a within-markers variance. It is completely different from the actual covariance matrix of the observed in the proximal segment coordinate system, and it is difficult to determine the form of using the observed covariance matrix. is a positive definite matrix of size of . If it is assumed an unstructured pattern, has no specific pattern.

. (8)

The diagonal is a 3×3 variance-covariance matrix of j frame of i marker, and the off diagonal is a 3×3 matrix of covariance between j and j' (j≠j ') frames of the i marker.

The positions of the markers obtained from the images correspond to the balanced data measured at the same time points and time intervals without missing values. The model in Equation 2 can be applied to both balanced and unbalanced data, and the reason for having subscript i is to distinguish a different number of observed values for each experimental unit (Hedeker & Gibbons, 2006). In the case of balanced data, it is convenient to handle the marker position vectors corresponding to all number of m in Equation (2) by stacking them in columns.

, (9)

Here, , , and are

, (10)

(N = mn). The mean of y is E(y)=μ =Xβ, and the variance of y is Var(y)=Σ, a matrix with a size of . The maximum likelihood estimator (MLE) of β that substitutes the estimated matrix () of Σ is as follows:

(11)

This is the best linear unbiased estimator (BLUE), and the expected value and variance of the estimated are

,

(12)

(Fitzmaurice et al., 2004; Gałecki & Burzykowski, 2013). It is important to determine the appropriate variance structure because the estimated value (Equation 11) and the variance of (Equation 12) vary according to the form of the covariance matrix . The number of parameters to be estimated in the covariance matrix in Equation 8 is . In general, estimates of the covariance matrix of the unstructured pattern should be considered as almost unidentifiable because the image data has a very large value of n. In conclusion, the research problem is to find a suitable covariance pattern model with fewer parameters.

2. Assumptions and covariance pattern models

Assuming independence and homoscedasticity of all measured vectors including components of (x,y,z), covariance matrix becomes Σ =. This is a classical general linear model and corresponds to the conventional transformation method for estimating the joint center, and in this case the problem becomes considerably simple. However, this study involves various models of the covariance matrix Σ, which shows the correlation and heteroscedasticity of the longitudinal data measured continuously over time.

First, the errors between the different markers were assumed to be independent of each other. Therefore, the measured values between the different markers were also independent. Second, it was assumed that all resulted from a common covariance matrix. This means that the size and component of the covariance matrices of all markers were the same (= ⋯ =) in the case of balanced data as in this study. These two assumptions are general assumptions in the covariance pattern models. In addition, two assumptions about the vector were added in this study. It was assumed that the components (x,y,z) of the vector were independent of each other and that the covariance matrices of the components (x,y,z) were the same. Therefore, the diagonal matrix and off diagonal matrix in Equation 8 are

(13)

where , , and the components (x,y,z) of the error vector are independent of each other and have the same variance. The covariance matrix of the components (x,y,z) is , and the covariance matrix of any component (x,y,z) is

(14)

Therefore, the covariance matrix of arbitrary marker i is

(15)

The covariance matrix Σ of all the markers is a block-diagonal matrix consisting of a single marker covariance matrix as follows:

(16)

The estimator () of the covariance matrix used in Equations 11 and 12 is mainly estimated by a restricted maximum likelihood (REML) (Fitzmaurice et al., 2004; Hedeker & Gibbons, 2006), and the REML function excluding the constant term by taking the logarithm is

(17)

Modelling of covariance matrix of this study uses the method by Pinheiro and Bates (2000). The covariance matrix was modeled according to variance structure and correlation structure. The 'gls' function (Pinheiro et al., 2017) of the free statistical program R(3.3.2) of 'nlme' package has various variance and correlation structure functions as shown in Table 1. In Equation 14, the variance structure is related to the diagonal components, and the correlation structure is related to the off-diagonal components. A combination of these two structures can produce more various covariance patterns.

Among these structures, the varPower of the variance structure and corAR1, corExp, corLin, and corSpher of the correlation structure are, respectively

(18)

and δ in the variance structure λ(⋅) are covariate and parameter, respectively. The s and ρ in the correlation structure are distance and parameter, respectively. Other functions are found in Pinheiro et al. (2017).

3. Experiment

An experiment of obtaining efficient estimation vectors of the hip joint center was conducted. The subjects were five males who were 179±2 cm in height, 73±3 kg in weight and 20.6±.89 years of age. All of them did not have any history of musculoskeletal diseases and agreed to participate in the experiment. The experiment was conducted at K University gymnasium with the approval of K University Institutional Review Board. The subjects attached nine markers (m=9) to the right thigh for the position vectors to be used in the model (Equation 9), and the position vectors were converted to the component of the pelvis coordinate system. To define the pelvic coordinate system (Wu et al., 2002), markers were attached to both the anterior superior iliac spine (ASIS) and posterior superior iliac spine (PSIS) of the pelvis. To define the femur coordinate system (Kim, 2017), markers were attached to greater trochanter (GT), lateral epicondyle (LE), and medial epicondyle (ME). The subjects performed various functional motions (De Rosario, Page, Besa, & Valera, 2013; Ojeda, Martínez-Reina, & Mayo, 2014) to estimate the hip joint center. The motion was a combination of three types of hip motions that allowed flexion/extension, internal/external rotation, and adduction/abduction. The first motion consisted of 30º of flexion, neutral position, and 30º of extension, and the second motion was composed of 30º of abduction and neutral position. The third motion consisted of 30º of flexion, and circumduction to 30º of extension and neutral position. The subjects performed three motions in order maintaining the upper body upright, both arms raised to the side, and the knee and ankle joint fixed (Figure 2).

Figure 2. A subject performing hip joint functional motions to estimate the CoR.

Modify Delete

The three-dimensional position vectors of the markers were obtained in centimeters by using the OptiTrack Motion Capture System of NaturalPoint (USA). Eight Prime13 cameras were used, and sampling was conducted at 120 frames per second. The raw three-dimensional position vectors did not undergo any filtering. All coordinate systems were as follows: +x for anterior, +y for superior, and +z for right. MATLAB ver 6.1 (MathWorks) was used for data processing, and 'gls' function of 'nlme' package of 'R(ver 3.3.2)' was used for inferences. The various models in Table 1 were applied to covariance matrices, and the GLM assuming was also applied. The appropriate covariance pattern model was selected on the basis of the likelihood ratio test, Akaike information criterion (AIC), and Bayesian information criterion (BIC). Significance level for the likelihood ratio test was α=.1 (Fitzmaurice et al., 2004; Hedeker & Gibbons, 2006).

Table 2 shows AICs and BICs of several selected models (Equation 18) applied to the covariance matrix and the results of the likelihood ratio tests of the models with GLM. The selected models had relatively small AIC and BIC values, and all the subjects showed the same results. All of the covariance pattern models showed better fit compared with of GLM, and all subjects showed the same results (p<.0001). The degree of freedom (df) means the number of parameters to be estimated. Even though the functions of variance related to time, varExp and varPower, and varExp, showed better fit than GLM, it was not shown in the table because AIC and BIC had larger values than varPower. Correlation structures also underwent the same procedure, and four functions (corAR1, corExp, corLin, and corSpher) were selected. In accordance with AIC and BIC, the best fit model was corSpher for subjects 1, 2, and 5, and corAR1 and corExp for subject 3. Subject 4 showed similar results of the four correlation structures.

Table 3 shows each subject's estimated joint center using the estimated covariance matrix of Table 2 and its standard error () of diagonal component in Equation 12). The variance structure function, varPower, showed a subtle difference in the estimated joint center and a slight reduction in standard error compared with GLM. However, in the case of the correlation structure function, the changes in the components (x,y,z) of the estimated joint center were observed, and estimated standard error was large. Compared with the GLM, the components (x,y,z) of the joint center were slightly shifted to the forward, upward, and medial sides. The results were the same for all the subjects except z for subject 2 and y for subject 3. In particular, in the case of subject 1, the change of the z value of the joint center was remarkable from approximately 11 to 5.5 cm. In all the subjects, the standard error was increased at a similar rate to the GLM, with a maximum increase of 30 times (subject 1 corAR1 and corExp).

Since the repeatedly measured position vectors are correlated with each other, and the variance of the position vectors is not constant, the joint center should be estimated by considering these characteristics.

The covariance pattern model is generated by finding an appropriate covariance matrix describing the correlation and heteroscedasticity of the error vector in the GLM of the position vector. Pinheiro and Bates (2000) referred to this as an extended linear model, and Hedeker and Gibbons (2006) explained it as an extension of the multivariate analysis of variance (MANOVA) that can be applied to an analysis of unbalanced data. In addition, the covariance pattern model may have been derived from a mixed model (Cnaan et al., 1997; Littell, Pendergast, & Natarajan, 2000; Liu, Rovine, & Molenaar, 2012; Wolfinger, 1993). If the markers are set as a random effect, the variance is divided into between- and within-markers variances. The combination of these two covariance matrices enables various covariance structures. Here, the between-markers variance is set to '0', and only the within-markers variance has various structures. Thus, the covariance pattern model is more general than the above-mentioned approaches in that it arbitrarily parametrizes the covariance structure (Jennrich & Schluchter, 1986). In most cases, the covariance pattern model is more appropriate for balanced data measured at the same interval of the longitudinal design (Fitzmaurice et al., 2004). This study focused on the covariance structure of error based on the model and the related joint center estimation.

If in Equation 14 is an unstructured pattern, the number of parameters is . In the image analysis, in the case of a high-speed camera, the number of covariance matrix parameters rapidly increases when the number of frame samples per second is increased or the operation time is long. In this study, we could not estimate the unstructured covariance matrix due to a computer memory problem. However, it is expected that even if the memory problem is resolved, the para- meters of unstructured covariance matrix may not be estimated parameters nor stable values. This is because the estimated unstructured pattern matrix estimation is unstable when the number of covariance matrix parameters is relatively large as compared with the number of markers (an experimental unit) (Fitzmaurice et al., 2004). Simplified modeling of the covariance matrix leads to a considerable improvement in the efficiency of estimated values (Jennrich & Schluchter, 1986). This study demonstrated efficiency by using a valid and simplified covariance structure model.

Several covariance pattern models (Jennrich & Schluchter, 1986; Pinheiro & Bates, 2000; Wolfinger, 1996) are provided in different statistical packages, and this study selected Pinheiro and Bates (2000) method, which can separately model the variance and the correlation structures. It was provided by the statistical program R, and the estimated joint center vector was obtained through various covariance matrix structures. The results of applying the varPower variance structure in Table 3 showed significant differences in the likelihood ratio test with the GLM, but the difference in the center of rotation and the standard errors was not significant. However, the correlation functions, corExp, corLin, and corSpher, showed markedly different results from the GLM. The estimated joint center value showed differences, and the standard errors also varied significantly. As the correlated values contained less information than the independent values, the standard error of the mean structure parameter increased when considering the correlation in the model. However, the estimated standard error at this time was less biased than the independent ones, which made it easier to control for the type I error (Gałecki & Burzykowski, 2013). Therefore, the greatly increased standard error in this study is considered to be an appropriate value considering the characteristics of repeatedly measurements. The differences in AIC and BIC between corAR1, corExp, corLin, and corSpher are so small that it is difficult to determine which is better fit. However, we can conclude that they are all better fits than the conventional GLM. In addition, the correlation structure rather than the variance structure better describes the characteristics of the longitudinal data for estimating the joint center.

It is common to identify the effect of response variables on changes in covariate when measurements are repeated. However, this study assumed that the position of the marker measured over time was a constant, as we presumed the general rigid body used in the mechanical analysis of the human body. The covariance structure is relatively simple as compared with the model for handling general longitudinal data. The covariance structure aims to obtain an efficient estimate β, but it has important implications as a covariance structure itself. This is because the covariance structure completely determines how the individual experimental units vary with deviations over time with respect to population mean (Schluchter, 1988). Therefore, if the appropriate covariance structure is determined, the STA or translational motion of the joints over time may be identified through this structure.

In this study, we did not obtain the result of applying the variance and correlation function at the same time. This may have been due to overparameterization (Cnaan et al., 1997) for making various structures of with limited data. This may have resulted from the limits of the capacity needed to calculate the matrix of massive data. Finding the optimal estimated parameters by applying various covariance matrix models is of great importance. Considering the characteristics of the image analysis data, modeling the covariance matrix into a complex structure is inefficient and impossible, so a simple and appropriate model should be selected. The data applied to the model in this study were vectors. A method that can have different correlation structures for each component of the vector is needed, as well as an estimation method considering a multilevel or hierarchical structure, and this should be further studied.

This study aimed to make appropriate inferences of the hip joint center through the GLM, considering the correlation and heteroscedasticity of the marker position vectors continuously measured from the high-speed camera images of the functional motion for estimating the hip center, and the following conclusions were obtained. The covariance pattern model of the GLM error considering the correlation and heteroscedasticity was a better fit than the GLM assuming classical independence and homoscedasticity. In addition, the correlation structure rather than the variance structure better explained the characteristics of the longitudinal data for estimating the joint center.