A study of Hip angle and Knee angle with Functional Principal Component Analysis
1 Objectives:
- Demonstrate how to model Biomechanic Data with Functional Principal Component Analysis
2 Biomechanics_dataset : Gait Dataset
- Gait Dataset is gait cycle for 39 boys
- Observation Time Normalized to \([0,1]\)
- The angles in the sagittal plane measured by the hip and knee when 39 boys walking gait cycle.
- Measured Objects : Hip Angle and Knee Angle Per Boy in each completed gait cycle.
3 Literature Reviews and Methodology : Skipped.
Details can refer to
- Schreiber and Moissenet (2019)
- Phinyomark et al. (2018)
- Sang, Wang, and Cao (2017)
- Olshen et al. (1989)
- https://cran.r-project.org/web/views/FunctionalData.html
4 Load Gatis Dataset
5 Construct Normalized Observation Time Domain
TimeRange = c(0,1)
LinearCoeff = c(0,(2*pi/diff(TimeRange))^2,0)
ObservationTimeNormalized = seq(from=0.025,to=0.975,by=0.05)6 Case 1: Hip Angle
6.1 Non-Parametric Fourier Basis Function Approach
\[\begin{align*} \phi_{k,0}(t) &:= 1 \\ \phi_{k,1}(t) &:= \cos(2\pi kt)\\ \phi_{k,2}(t) &:= \sin(2\pi kt)\\ f(t)&:=a_0 \phi_{k,0}(t) +\sum_{k=1}^{\infty}a_k \phi_{k,1}(t) + b_k \phi_{k,2}(t) \\ f(t)&:=a_0+\sum_{k=1}^{\infty}a_k\cos(2\pi kt)+b_k\sin(2\pi kt) \\ \end{align*}\]
HipAngle = gait[,,"Hip Angle"]
FourierBasisTimeFunction = create.fourier.basis(c(0, 1), nbasis=10, period=1)
LinearDifferentialOperator = vec2Lfd(LinearCoeff, TimeRange)6.1.1 Plot Hip Angle with FourierBasisTimeFunction
6.1.2 Find Optimal Non-Parametric \(\lambda\)
by minimizing generalized cross-validation error: \(gcv\)
\[gcv = \mathrm{\frac{n*SSE}{(n-df)^2}}\]
PowerArray = -10: +10
PowerArraySize = length(PowerArray)
SmothedBasisFunList = rep(NA,PowerArraySize)6.1.3 Plot Generalized Cross-Validation Error
for (hh in 1:PowerArraySize) {
lambda = 10^PowerArray[hh]
fdParobj = fdPar(FourierBasisTimeFunction, LinearDifferentialOperator, lambda)
SmothedList = smooth.basis(ObservationTimeNormalized,HipAngle,fdParobj,dfscale=1.2)
SmothedBasisFunList[hh] = sum(SmothedList$gcv)
}
LambdaMin_Hip = 10^(PowerArray[which.min(SmothedBasisFunList)])
LambdaMin_Hip## [1] 1e-07HipAngle.fit = smooth.basisPar(argvals=ObservationTimeNormalized,
y=HipAngle, fdobj=FourierBasisTimeFunction,
Lfdobj=LinearDifferentialOperator, lambda=LambdaMin_Hip)
HipAngle.fd = HipAngle.fit$fd6.1.4 Contour Map generated by Fitted Fourier Basis
TimeGrid = ObservationTimeNormalized
HipAngle.var = var.fd(HipAngle.fd)
tvvals = eval.bifd(TimeGrid, TimeGrid, HipAngle.var)
contour(TimeGrid, TimeGrid,tvvals,xlab='NormalizedTimeGrid',ylab='NormalizedTimeGrid',cex.lab=1.2,cex.axis=1.2, col = "red", lwd = 2)
title("NormalizedTimeGrid Contour Map: Hip Angle")
image.plot(TimeGrid,TimeGrid,tvvals,xlab='NormalizedTimeGrid',ylab='NormalizedTimeGrid',cex.lab=1.2,cex.axis=1.2)
title("NormalizedTimeGrid Contour Density Map: Hip Angle")
6.1.5 Compute Non-parametric Fouriers Basis FPCA Components
HipAngle.pcalist.fourier = pca.fd(HipAngle.fd, nharm=5)6.1.6 Plot PCA Components: Fourier Basis
plot(HipAngle.pcalist.fourier, col = "red", xlab = "Normalized Obervation Time of Hip Angle(Fourier Basis)")
CoeffFourierTimeBasisFun = HipAngle.pcalist.fourier$harmonics$coefs
CoeffFourierTimeBasisFun## PC1 PC2 PC3 PC4 PC5
## const 0.9780343233 -0.120024033 0.167580215 0.010586168 -0.020757521
## sin1 0.0423949210 0.883682906 0.363168243 0.270661962 -0.036633921
## cos1 0.2010314630 0.375500843 -0.859344066 -0.066835435 0.255618419
## sin2 -0.0304712561 -0.175501249 0.130329216 0.469518258 0.842818016
## cos2 -0.0065657211 0.168530492 0.271350021 -0.822387686 0.460498365
## sin3 0.0045666094 0.051233989 0.033170595 0.155439643 0.051775035
## cos3 -0.0062410975 -0.033968793 0.095376050 -0.009214563 -0.017802206
## sin4 -0.0135102772 0.011195221 0.018322643 0.024356108 0.086246348
## cos4 -0.0045931319 0.021740607 0.017217805 -0.020225279 0.002399875
## sin5 -0.0008410917 0.009293970 0.003316611 0.011100455 -0.006053436
## cos5 -0.0014631497 -0.005611583 0.009842979 0.004718771 0.002635010randfd1.fourier = fd(CoeffFourierTimeBasisFun[,1],FourierBasisTimeFunction)
randfd2.fourier = fd(CoeffFourierTimeBasisFun[,2],FourierBasisTimeFunction)
randfd3.fourier = fd(CoeffFourierTimeBasisFun[,3],FourierBasisTimeFunction)
randfd4.fourier = fd(CoeffFourierTimeBasisFun[,4],FourierBasisTimeFunction)
randfd5.fourier = fd(CoeffFourierTimeBasisFun[,5],FourierBasisTimeFunction)
fdvals1.fourier = eval.fd(ObservationTimeNormalized,randfd1.fourier)
fdvals2.fourier = eval.fd(ObservationTimeNormalized,randfd2.fourier)
fdvals3.fourier = eval.fd(ObservationTimeNormalized,randfd3.fourier)
fdvals4.fourier = eval.fd(ObservationTimeNormalized,randfd4.fourier)
fdvals5.fourier = eval.fd(ObservationTimeNormalized,randfd5.fourier)## [1] 0.72666812 0.12292826 0.08497140 0.03687697 0.01688092## [1] 0.9883257
6.2 Parametric Monomial Basis Function Approach
\[\begin{align*} \phi_{k}(t^{k}) &:= t^{k}\\ f(t) &:= \sum_{k= 0}^{ \infty} c_{t}\phi_{k}(t^{k}) \\ f(t) &:= c_{0} 1 +c_{1} t^{1} +c_{2} t^{2} + c_{3}t^{3} + c_{4}t^{4} + \dotsm \end{align*}\]
Notes:
- Parametric FPCA with B-Spline : set \(\lambda = 0\)
- int2Lfd(0) : penalizes the squared difference from 0 of the function
6.2.1 Find Optimal Degree of Order for Monomial Basis Function
DegreeOfOrder = 2:11
diferences = numeric(length(DegreeOfOrder))
for(kk in DegreeOfOrder){
timebasis.monomial = create.monomial.basis(c(0,1), kk+1)
fdParobj.monomial= fdPar(timebasis.monomial,int2Lfd(0),lambda=0)
HipAngle.pcalist.monomial = pca.fd(HipAngle.fd, nharm=3, fdParobj.monomial)
coef.monomial = HipAngle.pcalist.monomial$harmonics$coefs
randfd1.monomial = fd(coef.monomial[,1],timebasis.monomial)
randfd2.monomial = fd(coef.monomial[,2],timebasis.monomial)
randfd3.monomial = fd(coef.monomial[,3],timebasis.monomial)
fdvals1.monomial = eval.fd(ObservationTimeNormalized,randfd1.monomial)
fdvals2.monomial = eval.fd(ObservationTimeNormalized,randfd2.monomial)
fdvals3.monomial = eval.fd(ObservationTimeNormalized,randfd3.monomial)
diferences[kk-1] =
HipAngle.pcalist.fourier$varprop[1]*trapz(x=ObservationTimeNormalized, y=abs(fdvals1.monomial - fdvals1.fourier)) +
HipAngle.pcalist.fourier$varprop[2]*trapz(x=ObservationTimeNormalized, y=abs(fdvals2.monomial - fdvals2.fourier)) +
HipAngle.pcalist.fourier$varprop[3]*trapz(x=ObservationTimeNormalized, y=abs(fdvals3.monomial - fdvals3.fourier))
}diferences## [1] 0.2511946 0.2031234 0.1884396 0.1667319 0.3573146 0.1553945 0.1529367
## [8] 0.1497931 0.1484710 0.1447140
6.2.2 Compute Parametric Monomial Basis FPCA Componets
timebasis.monomial = create.monomial.basis(c(0,1), 5)
fdParobj.monomial= fdPar(timebasis.monomial,int2Lfd(0),lambda=0)
HipAngle.pcalist.monomial = pca.fd(HipAngle.fd, nharm=3, fdParobj.monomial)
coef.monomial = HipAngle.pcalist.monomial$harmonics$coefs
randfd1.monomial = fd(coef.monomial[,1],timebasis.monomial)
randfd2.monomial = fd(coef.monomial[,2],timebasis.monomial)
randfd3.monomial = fd(coef.monomial[,3],timebasis.monomial)- Therefore the optimal Degree of Order for Monomial Basis Function should be \(5\)
## [1] 0.7578349 0.1204087 0.0871829## [1] 0.9654264
6.2.3 Plot PCA Components: Monomial Basis
plot(HipAngle.pcalist.monomial, col = "blue", xlab = "Normalized Obervation Time of Hip Angle(Monomial Basis)")
6.3 FPCA Comparsion of Fourier Basis Vs Monomial Basis : Hip Angle
6.4 Summary : Hip Angle
10 Fourier Basis Functions: Non-Parametric Fourier Basis FPCA Components: \[NPFPCA_1 + NPFPCA_2 + NPFPCA_3 = 0.9345678\]
\(5\) Degree of Order of Monomial Basis Functions: Parametric Monomial Basis FPCA Components:
\[PFPCA_1 + PFPCA_2 + PFPCA_3 = 0.9654264\]
7 Case 2: Knee Angle
7.1 Non-Parametric Fourier Basis Function Approach
\[\begin{align*} \phi_{k,0}(t) &:= 1 \\ \phi_{k,1}(t) &:= \cos(2\pi kt)\\ \phi_{k,2}(t) &:= \sin(2\pi kt)\\ f(t)&:=a_0 \phi_{k,0}(t) +\sum_{k=1}^{\infty}a_k \phi_{k,1}(t) + b_k \phi_{k,2}(t) \\ f(t)&:=a_0+\sum_{k=1}^{\infty}a_k\cos(2\pi kt)+b_k\sin(2\pi kt) \\ \end{align*}\]
KneeAngle = gait[,,"Knee Angle"]
FourierBasisTimeFunction = create.fourier.basis(c(0, 1), nbasis=10, period=1)
LinearDifferentialOperator = vec2Lfd(LinearCoeff, TimeRange)7.1.1 Plot Knee Angle with FourierBasisTimeFunction
7.1.2 Find Optimal Non-Parametric \(\lambda\)
by minimizing generalized cross-validation error: \(gcv\)
\[gcv = \mathrm{\frac{n*SSE}{(n-df)^2}}\]
PowerArray = -10: +10
PowerArraySize = length(PowerArray)
SmothedBasisFunList = rep(NA,PowerArraySize)7.1.3 Plot Generalized Cross-Validation Error
for (hh in 1:PowerArraySize) {
lambda = 10^PowerArray[hh]
fdParobj = fdPar(FourierBasisTimeFunction, LinearDifferentialOperator, lambda)
SmothedList = smooth.basis(ObservationTimeNormalized,KneeAngle,fdParobj,dfscale=1.2)
SmothedBasisFunList[hh] = sum(SmothedList$gcv)
}
LambdaMin_Knee = 10^(PowerArray[which.min(SmothedBasisFunList)])
LambdaMin_Knee## [1] 1e-08KneeAngle.fit = smooth.basisPar(argvals=ObservationTimeNormalized,
y=KneeAngle, fdobj=FourierBasisTimeFunction,
Lfdobj=LinearDifferentialOperator, lambda=LambdaMin_Knee)
KneeAngle.fd = KneeAngle.fit$fd7.1.4 Contour Map generated by Fitted Fourier Basis
TimeGrid = ObservationTimeNormalized
KneeAngle.var = var.fd(KneeAngle.fd)
tvvals = eval.bifd(TimeGrid, TimeGrid, KneeAngle.var)
contour(TimeGrid, TimeGrid,tvvals,xlab='NormalizedTimeGrid',ylab='NormalizedTimeGrid',cex.lab=1.2,cex.axis=1.2, col = "red", lwd = 2)
title("NormalizedTimeGrid Contour Map: Knee Angle")
image.plot(TimeGrid,TimeGrid,tvvals,xlab='NormalizedTimeGrid',ylab='NormalizedTimeGrid',cex.lab=1.2,cex.axis=1.2)
title("NormalizedTimeGrid Contour Density Map: Knee Angle")
7.1.5 Compute Non-parametric Fouriers Basis FPCA Components
KneeAngle.pcalist.fourier = pca.fd(KneeAngle.fd, nharm=5)7.1.6 Plot PCA Components: Fourier Basis
plot(KneeAngle.pcalist.fourier, col = "red", xlab = "Normalized Obervation Time of Knee Angle(Fourier Basis)")
CoeffFourierTimeBasisFun = KneeAngle.pcalist.fourier$harmonics$coefs
CoeffFourierTimeBasisFun## PC1 PC2 PC3 PC4 PC5
## const 0.50723064 0.797212882 -0.067563421 -0.24870268 -0.18720279
## sin1 0.11438970 0.062043518 0.930967807 -0.02214217 0.24467106
## cos1 0.69317040 -0.240568524 -0.092487319 0.35476738 0.44986618
## sin2 -0.37898751 0.533570533 -0.067092536 0.64970943 0.34668347
## cos2 -0.23534786 0.106644463 -0.088730303 -0.61485282 0.63615406
## sin3 -0.04450655 0.022296451 0.232570144 0.04402370 -0.35854920
## cos3 -0.18979547 0.071343003 0.186657701 0.02891098 -0.05198628
## sin4 -0.04445262 0.005795555 -0.013481204 0.03602307 -0.02663664
## cos4 -0.06275549 0.003662455 0.136261256 0.07343932 -0.19127249
## sin5 0.03677769 -0.023532451 0.002109702 0.02201160 -0.09570999
## cos5 -0.06968987 0.021696057 0.011065911 0.04100031 -0.01072461randfd1.fourier = fd(CoeffFourierTimeBasisFun[,1],FourierBasisTimeFunction)
randfd2.fourier = fd(CoeffFourierTimeBasisFun[,2],FourierBasisTimeFunction)
randfd3.fourier = fd(CoeffFourierTimeBasisFun[,3],FourierBasisTimeFunction)
randfd4.fourier = fd(CoeffFourierTimeBasisFun[,4],FourierBasisTimeFunction)
randfd5.fourier = fd(CoeffFourierTimeBasisFun[,5],FourierBasisTimeFunction)
fdvals1.fourier = eval.fd(ObservationTimeNormalized,randfd1.fourier)
fdvals2.fourier = eval.fd(ObservationTimeNormalized,randfd2.fourier)
fdvals3.fourier = eval.fd(ObservationTimeNormalized,randfd3.fourier)
fdvals4.fourier = eval.fd(ObservationTimeNormalized,randfd4.fourier)
fdvals5.fourier = eval.fd(ObservationTimeNormalized,randfd5.fourier)## [1] 0.42816992 0.24992303 0.15592858 0.08602767 0.03903258## [1] 0.9590818
7.2 Parametric Monomial Basis Function Approach
\[\begin{align*} \phi_{k}(t^{k}) &:= t^{k}\\ f(t) &:= \sum_{k= 0}^{ \infty} c_{t}\phi_{k}(t^{k}) \\ f(t) &:= c_{0} 1 +c_{1} t^{1} +c_{2} t^{2} + c_{3}t^{3} + c_{4}t^{4} + \dotsm \end{align*}\]
Notes:
- Parametric FPCA with B-Spline : set \(\lambda = 0\)
- int2Lfd(0) : penalizes the squared difference from 0 of the function
7.2.1 Find Optimal Degree of Order for Monomial Basis Function
DegreeOfOrder = 2:11
diferences = numeric(length(DegreeOfOrder))
for(kk in DegreeOfOrder){
timebasis.monomial = create.monomial.basis(c(0,1), kk+1)
fdParobj.monomial= fdPar(timebasis.monomial,int2Lfd(0),lambda=0)
KneeAngle.pcalist.monomial = pca.fd(KneeAngle.fd, nharm=3, fdParobj.monomial)
coef.monomial = KneeAngle.pcalist.monomial$harmonics$coefs
randfd1.monomial = fd(coef.monomial[,1],timebasis.monomial)
randfd2.monomial = fd(coef.monomial[,2],timebasis.monomial)
randfd3.monomial = fd(coef.monomial[,3],timebasis.monomial)
fdvals1.monomial = eval.fd(ObservationTimeNormalized,randfd1.monomial)
fdvals2.monomial = eval.fd(ObservationTimeNormalized,randfd2.monomial)
fdvals3.monomial = eval.fd(ObservationTimeNormalized,randfd3.monomial)
diferences[kk-1] =
KneeAngle.pcalist.fourier$varprop[1]*trapz(x=ObservationTimeNormalized, y=abs(fdvals1.monomial - fdvals1.fourier)) +
KneeAngle.pcalist.fourier$varprop[2]*trapz(x=ObservationTimeNormalized, y=abs(fdvals2.monomial - fdvals2.fourier)) +
KneeAngle.pcalist.fourier$varprop[3]*trapz(x=ObservationTimeNormalized, y=abs(fdvals3.monomial - fdvals3.fourier))
}diferences## [1] 0.8823307 0.8935286 0.5174817 0.3408332 0.3409829 0.3209010 0.3070720
## [8] 0.3078678 0.3072810 0.3058688
7.2.2 Compute Parametric Monomial Basis FPCA Componets
optimaldfo = 10
timebasis.monomial = create.monomial.basis(c(0,1), optimaldfo)
fdParobj.monomial= fdPar(timebasis.monomial,int2Lfd(0),lambda=0)
KneeAngle.pcalist.monomial = pca.fd(KneeAngle.fd, nharm=3, fdParobj.monomial)
coef.monomial = KneeAngle.pcalist.monomial$harmonics$coefs
randfd1.monomial = fd(coef.monomial[,1],timebasis.monomial)
randfd2.monomial = fd(coef.monomial[,2],timebasis.monomial)
randfd3.monomial = fd(coef.monomial[,3],timebasis.monomial)- Therefore the optimal Degree of Order for Monomial Basis Function should be \(10\)
## [1] 0.4311328 0.2525375 0.1567733## [1] 0.8404436
7.2.3 Plot PCA Components: Monomial Basis
plot(KneeAngle.pcalist.monomial, col = "blue", xlab = "Normalized Obervation Time of Knee Angle(Monomial Basis)")
7.3 FPCA Comparsion of Fourier Basis Vs Monomial Basis : Knee Angle
7.4 Summary : Knee Angle
- 10 Fourier Basis Functions: Non-Parametric Fourier Basis FPCA Components:
\[NPFPCA_1 + NPFPCA_2 + NPFPCA_3 = 0.8340215\]
- \(5\) Degree of Order of Monomial Basis Functions: Parametric Monomial Basis FPCA Components:
\[PFPCA_1 + PFPCA_2 + PFPCA_3 = 0.8404436\]
8 Co-movement: Hip and Knee
8.1 Bivariate Plot Per Each Boy Trajectory
8.2 Bivariate plot
8.3 Harmonics Cycle Plots For FPCA
Olshen, Richard A., Edmund N. Biden, Marilynn P. Wyatt, and David H. Sutherland. 1989. “Gait Analysis and the Bootstrap.” The Annals of Statistics 17 (4): 1419–40. https://doi.org/10.1214/aos/1176347372.
Phinyomark, Angkoon, Giovanni Petri, Esther Ibáñez-Marcelo, Sean T. Osis, and Reed Ferber. 2018. “Analysis of Big Data in Gait Biomechanics: Current Trends and Future Directions.” Journal of Medical and Biological Engineering 38: 244–60.
Sang, Peijun, Liangliang Wang, and Jiguo Cao. 2017. “Parametric Functional Principal Component Analysis.” Biometrics 73 (September): 802–10.
Schreiber, Céline, and Florent Moissenet. 2019. “A Multimodal Dataset of Human Gait at Different Walking Speeds Established on Injury-Free Adult Participants.” Scientific Data 6 (1): 111. https://doi.org/10.1038/s41597-019-0124-4.