Analysis of CE
Yishi Wang
1/10/2021
1 Preprocessing
Read the original data from excel.
Organize the data into new format with column names id, age, body mass, height, nc1, nc7, week, group, and response variables.
## id age bm h nc1 nc7 resp week group
## 1 1 22 72.4 183.2 36.4 39.0 6.0 1 1
## 2 1 22 72.4 183.2 36.4 39.0 5.2 1 1
## 3 1 22 72.4 183.2 36.4 39.0 5.8 1 1
## 4 2 21 81.7 181.9 35.6 40.7 4.8 1 0
## 5 2 21 81.7 181.9 35.6 40.7 4.4 1 0
## 6 2 21 81.7 181.9 35.6 40.7 7.0 1 0- Creat functions to show the averages of a reapsonse variables in each week over time for every subject.
2 Data Analysis of CE
2.1 Discriptive and Visual Analysis
Here is an overview of the top six records for CE data, and individual subject’s CE measurement over the experiment period.
## id age bm h nc1 nc7 resp week group
## 1 1 22 72.4 183.2 36.4 39.0 6.0 1 1
## 2 1 22 72.4 183.2 36.4 39.0 5.2 1 1
## 3 1 22 72.4 183.2 36.4 39.0 5.8 1 1
## 4 2 21 81.7 181.9 35.6 40.7 4.8 1 0
## 5 2 21 81.7 181.9 35.6 40.7 4.4 1 0
## 6 2 21 81.7 181.9 35.6 40.7 7.0 1 0
We may notice that some subject only have records of week1 (Sub 8,10,23,27,29,32).
2.2 Visual of changes in average over time.
Since we have three repetitive measurements for each of the 3 weeks, the average of three measuremens with each week of each subject are obtained, and plotted below. As shown in the above figure, the Treatment group (red) tends to have more upward moving trend than the Control group (Blue).
##Statistical assessment of the different patterns in CE between the two groups. ### Anslysis with simple linear regressions For each subject, we first conduct a simple linear regression of CE with week as the covariate. A pair of slope and intercept is then obtained for each subject. The following figure is based on those pairs of slopes and intercepts. It is observed that the Topspin group may have greater slopes, which is consistent with the figure in 2.2. As to the intercept, it is hard to tell a difference between the two groups.
##
## Welch Two Sample t-test
##
## data: slope.ml.ce[uniq.group == 1] and slope.ml.ce[uniq.group == 0]
## t = -2.6911, df = 23.538, p-value = 0.01289
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.41095913 -0.05399814
## sample estimates:
## mean of x mean of y
## 0.06031746 0.29279609##
## Welch Two Sample t-test
##
## data: inter.ml.ce[uniq.group == 1] and inter.ml.ce[uniq.group == 0]
## t = -0.13572, df = 29.875, p-value = 0.893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.10663 1.84413
## sample estimates:
## mean of x mean of y
## 4.547321 4.678571The above two two-sample tests confirms the impression we have from the figure. The first test is comparing the slopes of the two groups, by using a two-sample t-test. With a p-value of 0.01289, we can conclude that there is a statistically significant difference in the linear slopes of increase for Cervical Extension, between the two groups. On the other hand, the second test suggests that there is no significant difference in the intercepts for the two groups. This also indicates that the two groups are similar from week zero in terms of the initial CE values.
###Modeling with Linear Mixed Model for longitudinal data. The above analysis is simple and straightfroward. However, it is based on heavy assumptions that each subject’s CE grows linearly, and variations of observations around each subject’s weekly average measurements are also ignored. Both of the first two figures in section 2 suggests that these two assumptions are not valid. Also, other variables measured in the study besides week and interested response variables are not involved in the modeling. This leads us to consider Liner Mixed Model, a popular approach in modeling longitudinal data. (REference here). Given the background of the study, the full model we consider is the following:
\[\begin{align} CE_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_a Age_i+\beta_b BM_{ij}+\beta_h Height_{ij}\nonumber\\ &+\beta_{nc1} NC1_{ijk}+\beta_{nc7} NC7_{ijk}+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\] subject \(i\), week \(j\) and repetition \(k\)
2.2.0.1 Full model
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + age + bm + h + nc1 + nc7 + (week | id)
## Data: dat.ce
##
## REML criterion at convergence: 840.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7527 -0.5198 0.0071 0.4745 2.9816
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.50133 1.5816
## week 0.03229 0.1797 -0.36
## Residual 0.87792 0.9370
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -27.12723 4.86682 -5.574
## group.faTopspin -1.17581 0.62941 -1.868
## week 0.10429 0.05561 1.875
## age 0.12619 0.08293 1.522
## bm 0.06010 0.03464 1.735
## h 0.15577 0.03497 4.455
## nc1 0.14598 0.15991 0.913
## nc7 -0.17895 0.10605 -1.687
## group.faTopspin:week 0.18084 0.07782 2.324
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week age bm h nc1 nc7
## grop.fTpspn -0.011
## week -0.172 0.259
## age -0.247 0.010 -0.075
## bm 0.508 -0.077 -0.130 -0.145
## h -0.766 0.056 0.118 0.146 -0.304
## nc1 -0.192 -0.197 0.142 -0.323 -0.448 -0.160
## nc7 -0.006 0.171 -0.091 0.164 -0.097 -0.224 -0.481
## grp.fTpspn: 0.124 -0.393 -0.713 0.058 0.077 -0.117 -0.054 0.056Based on the output, besides Week and Height, another possibly significant variable in the interaction effect between Group and Week. Let us examine the following model without Age, NC1 and NC7.
2.2.0.2 Reduced model
\[\begin{align} CE_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_b BM_{ij}+\beta_h Height_{ij}\nonumber\\ &+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\]
## F-test with Kenward-Roger approximation; computing time: 0.36 sec.
## large : resp ~ group.fa + week + age + bm + h + nc1 + nc7 + (week | id) +
## group.fa:week
## small : resp ~ group.fa * week + bm + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 2.0147 3.0000 49.3780 0.99388 0.124In the above output, to test the fixed effect terms, the Kenward-Roger adjusted F-test is used. With a p-value of 0.124, we can conclude that the full model has a similar fit with the data as the 1st full model does. As a common practice in stats, when there is no significant difference between a full model and a simplified model, the simplified model is prefered due to a smaller number of parameters needed.
The following output shows the summary results of linear mixed model based on the equation in 3.2.2.
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + bm + h + (week | id)
## Data: dat.ce
##
## REML criterion at convergence: 839.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6142 -0.5060 0.0134 0.4937 2.9954
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.82027 1.6794
## week 0.03502 0.1871 -0.41
## Residual 0.87689 0.9364
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -24.27062 4.61298 -5.261
## group.faTopspin -0.91997 0.64637 -1.423
## week 0.09934 0.05673 1.751
## bm 0.07759 0.02726 2.846
## h 0.13815 0.03314 4.168
## group.faTopspin:week 0.17953 0.08005 2.243
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week bm h
## grop.fTpspn -0.065
## week -0.177 0.320
## bm 0.401 -0.201 -0.113
## h -0.948 0.066 0.148 -0.663
## grp.fTpspn: 0.148 -0.447 -0.713 0.105 -0.132Based on the output above, it shows that the Topspin group has less CE to begin with (by -0.92). However, for every week, subjects’ CE values in the Topspin group increase more than the ones from the Control group (by 0.18, which is reasonable to believe that after 10 week of experiment, on average, CE from the Topspin group is 0.875 higher than the CE from the Control group)
The following analysis is conducted for testing the significance of the interaction term between week and group. As the p-value suggests, the interaction term is significant with a p-value of 0.03535, and thus it is indispensable in the model.
## F-test with Kenward-Roger approximation; computing time: 0.18 sec.
## large : resp ~ group.fa + week + bm + h + (week | id) + group.fa:week
## small : resp ~ group.fa + week + bm + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 4.9494 1.0000 25.0153 1 0.03535 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2.0.3 Testing the significance of random terms
## Computing bootstrap confidence intervals ...##
## 4 warning(s): Model failed to converge with max|grad| = 0.00507303 (tol = 0.002, component 1) (and others)## 2.5 % 97.5 %
## .sig01 1.107312131 2.1952322
## .sig02 -0.735888856 0.1167388
## .sig03 0.114724733 0.2505774
## .sigma 0.837986244 1.0308135
## (Intercept) -34.362606659 -15.1164327
## group.faTopspin -2.228711915 0.2831739
## week -0.004925324 0.2151514
## bm 0.025556625 0.1308625
## h 0.071313699 0.2077637
## group.faTopspin:week 0.017663161 0.3436443A parametric bootstrap approach is used here to evaluate the significance of the random terms.The sd of the intercept for each subject is within (1.107, 2.195) with 95% confidence, whereas (0.115, 0.251) for the sd of slope for week. Since both of the CI’s are above zero, the random terms are kept in the model.
2.2.0.4 Diagnosis plot for residuals.
The following analysis is to evaluate the normality assumption in the error term of the 1st reduced model. Based on the output, we can claim that the normality assumption is valid for the anslysis of the dataset.
3 Data Analysis of CF
3.1 Discriptive and Visual Analysis
Here is an overview of the top six records for CF data, and individual subject’s CF measurement over the experiment period. Notice that the resp is now for CF.
## id age bm h nc1 nc7 resp week group
## 1 1 22 72.4 183.2 36.4 39.0 4.6 1 1
## 2 1 22 72.4 183.2 36.4 39.0 5.2 1 1
## 3 1 22 72.4 183.2 36.4 39.0 5.2 1 1
## 4 2 21 81.7 181.9 35.6 40.7 5.6 1 0
## 5 2 21 81.7 181.9 35.6 40.7 4.4 1 0
## 6 2 21 81.7 181.9 35.6 40.7 5.2 1 0
We may notice that some subject only have records of week1 (Sub 8,10,23,27,29,32).
3.2 Visual of changes in average over time.
##Statistical assessment of the different patterns in CF between the two groups.
3.2.1 Anslysis with simple linear regressions
Similar with the CE data, we first conduct a simple linear regression of CF with week as the covariate.
It is observed that the Topspin group may have greater slopes, which is consistent with the figure in 3.2. As to the intercept, it is hard to tell a difference between the two groups.
##
## Welch Two Sample t-test
##
## data: slope.ml.cf[uniq.group == 1] and slope.ml.cf[uniq.group == 0]
## t = -3.2028, df = 23.442, p-value = 0.003889
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.45258427 -0.09759888
## sample estimates:
## mean of x mean of y
## -0.03443223 0.24065934##
## Welch Two Sample t-test
##
## data: inter.ml.cf[uniq.group == 1] and inter.ml.cf[uniq.group == 0]
## t = -0.056766, df = 29.49, p-value = 0.9551
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.048400 1.937686
## sample estimates:
## mean of x mean of y
## 4.133631 4.188988The above two two-sample tests confirms the impression we have from the figure. The first test is comparing the slopes of the two groups, by using a two-sample t-test. With a p-value of 0.0038895, we can conclude that there is a statistically significant difference in the linear slopes of increase for Cervical Extension, between the two groups. On the other hand, the second test suggests that there is no significant difference (p-value of 0.9551147)in the intercepts for the two groups. This also indicates that the two groups are similar from week zero in terms of the initial CF values.
###Modeling with Linear Mixed Model for longitudinal data. In this section, a Liner Mixed Model is considered as following:
\[\begin{align} CF_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_a Age_i+\beta_b BM_{ij}+\beta_h Height_{ij}\nonumber\\ &+\beta_{nc1} NC1_{ijk}+\beta_{nc7} NC7_{ijk}+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\] subject \(i\), week \(j\) and repetition \(k\)
3.2.1.1 Full model
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + age + bm + h + nc1 + nc7 + (week | id)
## Data: dat.cf
##
## REML criterion at convergence: 739.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6785 -0.4640 -0.0140 0.4474 3.6366
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 3.02834 1.7402
## week 0.03543 0.1882 -0.21
## Residual 0.50223 0.7087
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -28.324629 5.226138 -5.420
## group.faTopspin -1.156667 0.667650 -1.732
## week 0.007882 0.056435 0.140
## age 0.110057 0.090992 1.210
## bm 0.032732 0.037844 0.865
## h 0.165865 0.034081 4.867
## nc1 0.151662 0.161758 0.938
## nc7 -0.145946 0.097013 -1.504
## group.faTopspin:week 0.229543 0.078487 2.925
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week age bm h nc1 nc7
## grop.fTpspn 0.017
## week -0.176 0.149
## age -0.267 -0.005 -0.071
## bm 0.502 -0.077 -0.150 -0.167
## h -0.769 0.040 0.123 0.150 -0.347
## nc1 -0.277 -0.175 0.173 -0.283 -0.490 -0.086
## nc7 -0.005 0.152 -0.090 0.153 -0.062 -0.190 -0.473
## grp.fTpspn: 0.119 -0.237 -0.712 0.035 0.087 -0.108 -0.075 0.054Based on the output, let us examine the following model without Age, Body Mass, NC1 and NC7.
3.2.1.2 Reduced model
\[\begin{align} CF_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_h Height_{ij}\nonumber\\ &+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\]
## F-test with Kenward-Roger approximation; computing time: 0.16 sec.
## large : resp ~ group.fa + week + age + bm + h + nc1 + nc7 + (week | id) +
## group.fa:week
## small : resp ~ group.fa * week + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 1.9347 4.0000 53.1702 0.99368 0.1181In the above output, to test the fixed effect terms, the Kenward-Roger adjusted F-test is used. With a p-value of 0.118, we can conclude that the full model has a similar fit with the data as the reduced model does.
The following output shows the summary results of linear mixed model based on the equation in 3.3.2.2.
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + h + (week | id)
## Data: dat.cf
##
## REML criterion at convergence: 733.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5999 -0.4972 -0.0166 0.4176 3.6619
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 3.53654 1.8806
## week 0.03424 0.1850 -0.06
## Residual 0.49836 0.7059
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -27.64165 4.65417 -5.939
## group.faTopspin -0.66373 0.68727 -0.966
## week 0.01093 0.05460 0.200
## h 0.18772 0.02732 6.872
## group.faTopspin:week 0.22906 0.07694 2.977
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week h
## grop.fTpspn 0.020
## week -0.133 0.080
## h -0.995 -0.094 0.121
## grp.fTpspn: 0.095 -0.121 -0.710 -0.086Based on the output above, it shows that the Topspin group has less CF to begin with (by -0.664). However, for every week, subjects’ CF values in the Topspin group increase more than the ones from the Control group (by around 0.229, which is reasonable to believe that after 10 week of experiment, on average, CF from the Topspin group is 1.627 higher than the CF from the Control group)
## F-test with Kenward-Roger approximation; computing time: 0.19 sec.
## large : resp ~ group.fa + week + h + (week | id) + group.fa:week
## small : resp ~ group.fa + week + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 8.7188 1.0000 24.0554 1 0.006931 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The above analysis is conducted for testing the significance of the interaction term between week and group. As the p-value suggests, the interaction term is significant with a p-value of 0.007, and thus it is indispensable in the model.
3.2.1.3 Testing the significance of random terms
## Computing bootstrap confidence intervals ...##
## 9 warning(s): Model failed to converge with max|grad| = 0.00246034 (tol = 0.002, component 1) (and others)## 2.5 % 97.5 %
## .sig01 1.30648688 2.3765516
## .sig02 -0.46783595 0.4014254
## .sig03 0.11885705 0.2451362
## .sigma 0.63173781 0.7778349
## (Intercept) -36.96091455 -19.1233081
## group.faTopspin -2.13470204 0.7391242
## week -0.08863348 0.1119399
## h 0.13762629 0.2415658
## group.faTopspin:week 0.08463467 0.3949762A parametric bootstrap approach is used here to evaluate the significance of the random terms. The sd of the intercept for each subject is within (1.306, 2.377) with 95% confidence, whereas (0.119, 0.245) for the sd of slope for week. Since both of the CI’s are above zero, the random terms are kept in the model.
3.2.1.4 Diagnosis plot for residuals.
The following analysis is to evaluate the normality assumption in the error term of the reduced model. Based on the output, we can claim that the normality assumption is valid for the anslysis of the dataset.
4 Data Analysis of LCLF
4.1 Discriptive and Visual Analysis
Here is an overview of the top six records for LCLF data, and individual subject’s LCLF measurement over the experiment period. Notice that the resp is now for LCLF.
## id age bm h nc1 nc7 resp week group
## 1 1 22 72.4 183.2 36.4 39.0 7.4 1 1
## 2 1 22 72.4 183.2 36.4 39.0 6.4 1 1
## 3 1 22 72.4 183.2 36.4 39.0 7.6 1 1
## 4 2 21 81.7 181.9 35.6 40.7 7.2 1 0
## 5 2 21 81.7 181.9 35.6 40.7 8.4 1 0
## 6 2 21 81.7 181.9 35.6 40.7 7.4 1 0
We may notice that some subject only have records of week1 (Sub 8,10,23,27,29,32).
4.2 Visual of changes in average over time.
##Statistical assessment of the different patterns in LCLF between the two groups.
4.2.1 Anslysis with simple linear regressions
Similar with the CE data, we first conduct a simple linear regression of LCLF with week as the covariate.
It is observed that the two groups may have similar slopes, which is consistent with the figure in 4.2. As to the intercept, it is also hard to tell a difference between the two groups.
##
## Welch Two Sample t-test
##
## data: slope.ml.lclf[uniq.group == 1] and slope.ml.lclf[uniq.group == 0]
## t = -0.55406, df = 23.144, p-value = 0.5849
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2452845 0.1416215
## sample estimates:
## mean of x mean of y
## 0.1154457 0.1672772##
## Welch Two Sample t-test
##
## data: inter.ml.lclf[uniq.group == 1] and inter.ml.lclf[uniq.group == 0]
## t = -0.6325, df = 29.25, p-value = 0.532
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.429342 1.808806
## sample estimates:
## mean of x mean of y
## 4.950446 5.760714The above two two-sample tests confirms the impression we have from the figure. The first test is comparing the slopes of the two groups, by using a two-sample t-test. With a p-value of 0.5848516, we can conclude that there is no statistically significant difference in the linear slopes of increase for LCLF, between the two groups. Also, the second test suggests that there is no significant difference (p-value of 0.5319727)in the intercepts for the two groups. This also indicates that the two groups are similar from week zero in terms of the initial LCLF values.
###Modeling with Linear Mixed Model for longitudinal data. In this section, a Liner Mixed Model is considered as following:
\[\begin{align} LCLF_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_a Age_i+\beta_b BM_{ij}+\beta_h Height_{ij}\nonumber\\ &+\beta_{nc1} NC1_{ijk}+\beta_{nc7} NC7_{ijk}+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\] subject \(i\), week \(j\) and repetition \(k\)
4.2.1.1 Full model
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + age + bm + h + nc1 + nc7 + (week | id)
## Data: dat.lclf
##
## REML criterion at convergence: 751
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.41205 -0.51483 -0.05491 0.59605 2.98978
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 5.41586 2.3272
## week 0.05327 0.2308 -0.53
## Residual 0.48749 0.6982
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -22.65921 5.73133 -3.954
## group.faTopspin -0.52905 0.87109 -0.607
## week 0.12886 0.06632 1.943
## age 0.11551 0.10850 1.065
## bm 0.12674 0.04257 2.977
## h 0.11778 0.03833 3.073
## nc1 0.06232 0.18497 0.337
## nc7 -0.15714 0.10678 -1.472
## group.faTopspin:week 0.03722 0.09299 0.400
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week age bm h nc1 nc7
## grop.fTpspn 0.003
## week -0.183 0.335
## age -0.231 -0.015 -0.082
## bm 0.452 -0.070 -0.124 -0.146
## h -0.757 0.038 0.116 0.104 -0.291
## nc1 -0.270 -0.154 0.140 -0.283 -0.490 -0.117
## nc7 0.039 0.137 -0.083 0.129 -0.057 -0.207 -0.470
## grp.fTpspn: 0.126 -0.496 -0.713 0.080 0.073 -0.114 -0.057 0.048Based on the output, let us examine the following model without Age, NC1, NC7, Group, and its interaction term with Week.
4.2.1.2 Reduced model
\[\begin{align} LCLF_{ijk}=&\mu+\beta_b BM_{ij} +\beta_h Height_{ij}\nonumber\\ &+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\]
## F-test with Kenward-Roger approximation; computing time: 0.27 sec.
## large : resp ~ group.fa + week + age + bm + h + nc1 + nc7 + (week | id) +
## group.fa:week
## small : resp ~ week + bm + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 0.7132 5.0000 47.3507 0.97206 0.6166In the above output, to test the fixed effect terms, the Kenward-Roger adjusted F-test is used. With a p-value of 0.617, we can conclude that the full model has a similar fit with the data as the reduced model does.
The following output shows the summary results of linear mixed model based on the equation in 4.3.2.2.
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ week + bm + h + (week | id)
## Data: dat.lclf
##
## REML criterion at convergence: 746
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39699 -0.55194 -0.01713 0.59047 2.97239
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 4.96332 2.2279
## week 0.05239 0.2289 -0.50
## Residual 0.48911 0.6994
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -21.19219 5.04938 -4.197
## week 0.14542 0.04583 3.173
## bm 0.12663 0.03088 4.101
## h 0.10127 0.03555 2.848
##
## Correlation of Fixed Effects:
## (Intr) week bm
## week -0.104
## bm 0.289 -0.047
## h -0.933 0.072 -0.608## F-test with Kenward-Roger approximation; computing time: 0.19 sec.
## large : resp ~ group.fa + week + bm + h + (week | id) + group.fa:week
## small : resp ~ group.fa + week + bm + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 0.1351 1.0000 25.3469 1 0.7163The above analysis is conducted for testing the significance of the interaction term between week and group. As the p-value suggests, the interaction term is not significant with a p-value of 0.716, and thus we will keep using the reduced model in 4.3.2.2.
4.2.1.3 Testing the significance of random terms
## Computing bootstrap confidence intervals ...##
## 3 warning(s): Model failed to converge with max|grad| = 0.0028143 (tol = 0.002, component 1) (and others)## 2.5 % 97.5 %
## .sig01 1.60305416 2.7685926
## .sig02 -0.76964464 -0.0853969
## .sig03 0.15573684 0.2965496
## .sigma 0.62636664 0.7707887
## (Intercept) -31.92274773 -11.1378423
## week 0.06497780 0.2397044
## bm 0.06490695 0.1905066
## h 0.03113543 0.1787358A parametric bootstrap approach is used here to evaluate the significance of the random terms. The sd of the intercept for each subject is within (1.603, 2.769) with 95% confidence, whereas (0.156, 0.297) for the sd of slope for week. Since both of the CI’s are above zero, the random terms are kept in the model.
4.2.1.4 Diagnosis plot for residuals.
The following analysis is to evaluate the normality assumption in the error term of the reduced model. Based on the output, we can claim that the normality assumption is valid for the anslysis of the dataset.
5 Data Analysis of RCLF
5.1 Discriptive and Visual Analysis
Here is an overview of the top six records for RCLF data, and individual subject’s RCLF measurement over the experiment period. Notice that the resp is now for RCLF.
## id age bm h nc1 nc7 resp week group
## 1 1 22 72.4 183.2 36.4 39.0 5.0 1 1
## 2 1 22 72.4 183.2 36.4 39.0 7.2 1 1
## 3 1 22 72.4 183.2 36.4 39.0 7.0 1 1
## 4 2 21 81.7 181.9 35.6 40.7 7.6 1 0
## 5 2 21 81.7 181.9 35.6 40.7 9.8 1 0
## 6 2 21 81.7 181.9 35.6 40.7 9.4 1 0
We may notice that some subject only have records of week1 (Sub 8,10,23,27,29,32).
5.2 Visual of changes in average over time.
##Statistical assessment of the different patterns in RCLF between the two groups.
5.2.1 Anslysis with simple linear regressions
Similar with the CE data, we first conduct a simple linear regression of RCLF with week as the covariate.
It is observed that the two groups may have similar slopes, which is consistent with the figure in 5.2. As to the intercept, it is also hard to tell a difference between the two groups.
##
## Welch Two Sample t-test
##
## data: slope.ml.rclf[uniq.group == 1] and slope.ml.rclf[uniq.group == 0]
## t = -1.5534, df = 23.999, p-value = 0.1334
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2733847 0.0385862
## sample estimates:
## mean of x mean of y
## 0.0525641 0.1699634##
## Welch Two Sample t-test
##
## data: inter.ml.rclf[uniq.group == 1] and inter.ml.rclf[uniq.group == 0]
## t = -0.83401, df = 28.885, p-value = 0.4111
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.467619 1.458989
## sample estimates:
## mean of x mean of y
## 4.863542 5.867857The above two two-sample tests confirms the impression we have from the figure. The first test is comparing the slopes of the two groups, by using a two-sample t-test. With a p-value of 0.1334277, we can conclude that there is no statistically significant difference in the linear slopes of increase for RCLF, between the two groups. Also, the second test suggests that there is no significant difference (p-value of 0.4111218)in the intercepts for the two groups. This also indicates that the two groups are similar from week zero in terms of the initial RCLF values.
###Modeling with Linear Mixed Model for longitudinal data. In this section, a Liner Mixed Model is considered as following:
\[\begin{align} RCLF_{ijk}=&\mu+\beta_g Group_i+\beta_w Week_j+\beta_{gw} Group_i\times Week_j +\beta_a Age_i+\beta_b BM_{ij}+\beta_h Height_{ij}\nonumber\\ &+\beta_{nc1} NC1_{ijk}+\beta_{nc7} NC7_{ijk}+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\] subject \(i\), week \(j\) and repetition \(k\)
5.2.1.1 Full model
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ group.fa * week + age + bm + h + nc1 + nc7 + (week | id)
## Data: dat.rclf
##
## REML criterion at convergence: 785.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.74509 -0.50407 0.00361 0.43111 2.49559
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 4.17584 2.0435
## week 0.03021 0.1738 -0.42
## Residual 0.63025 0.7939
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -24.62897 5.48338 -4.492
## group.faTopspin -0.50916 0.77700 -0.655
## week 0.08350 0.05306 1.574
## age 0.02727 0.10092 0.270
## bm 0.10171 0.04046 2.514
## h 0.09975 0.03660 2.725
## nc1 0.27802 0.17676 1.573
## nc7 -0.12085 0.10235 -1.181
## group.faTopspin:week 0.09337 0.07375 1.266
##
## Correlation of Fixed Effects:
## (Intr) grp.fT week age bm h nc1 nc7
## grop.fTpspn 0.011
## week -0.205 0.264
## age -0.235 -0.011 -0.094
## bm 0.466 -0.071 -0.159 -0.140
## h -0.761 0.038 0.138 0.116 -0.302
## nc1 -0.276 -0.164 0.183 -0.289 -0.498 -0.107
## nc7 0.039 0.145 -0.104 0.136 -0.051 -0.212 -0.472
## grp.fTpspn: 0.141 -0.404 -0.715 0.072 0.093 -0.130 -0.077 0.062Based on the output, let us examine the following model without Age, NC1, NC7, Group, and its interaction term with Week.
5.2.1.2 Reduced model
\[\begin{align} RCLF_{ijk}=&\mu+\beta_w Week_j+\beta_b BM_{ij} +\beta_h Height_{ij}\nonumber\\ &+\gamma_{0i}+\gamma_{1i}Week_i+\epsilon_{ijk} \end{align}\]
## F-test with Kenward-Roger approximation; computing time: 0.23 sec.
## large : resp ~ group.fa + week + age + bm + h + nc1 + nc7 + (week | id) +
## group.fa:week
## small : resp ~ week + bm + h + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 0.950 5.000 43.812 0.97179 0.4587In the above output, to test the fixed effect terms, the Kenward-Roger adjusted F-test is used. With a p-value of 0.459, we can conclude that the full model has a similar fit with the data as the reduced model does.
The following output shows the summary results of linear mixed model based on the equation in 4.3.2.2. Since the group covariate is no longer in the model, we may claim that there is no statistical difference in RCLF between the two groups.
## Linear mixed model fit by REML ['lmerMod']
## Formula: resp ~ week + bm + h + (week | id)
## Data: dat.rclf
##
## REML criterion at convergence: 781.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.78645 -0.52548 -0.02848 0.43006 2.50922
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 4.2686 2.0660
## week 0.0318 0.1783 -0.42
## Residual 0.6256 0.7910
## Number of obs: 252, groups: id, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -22.20098 4.90695 -4.524
## week 0.12031 0.03727 3.228
## bm 0.13313 0.02993 4.448
## h 0.10464 0.03455 3.029
##
## Correlation of Fixed Effects:
## (Intr) week bm
## week -0.118
## bm 0.290 -0.061
## h -0.934 0.094 -0.608## F-test with Kenward-Roger approximation; computing time: 0.17 sec.
## large : resp ~ group.fa + week + bm + (week | id) + group.fa:week
## small : resp ~ group.fa + week + bm + (week | id)
## stat ndf ddf F.scaling p.value
## Ftest 3.0041 1.0000 24.5206 1 0.09561 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The above analysis is conducted for testing the significance of the interaction term between week and group. As the p-value suggests, the interaction term is not significant with a p-value of 0.096, and thus we will keep using the reduced model in 5.3.2.2.
5.2.1.3 Testing the significance of random terms
## Computing bootstrap confidence intervals ...##
## 3 warning(s): Model failed to converge with max|grad| = 0.00257739 (tol = 0.002, component 1) (and others)## 2.5 % 97.5 %
## .sig01 1.45935335 2.5725863
## .sig02 -0.73402278 0.0413978
## .sig03 0.11415577 0.2362640
## .sigma 0.70835216 0.8716684
## (Intercept) -32.96559783 -12.5102330
## week 0.05299239 0.1962544
## bm 0.07332558 0.1922281
## h 0.03681866 0.1824160A parametric bootstrap approach is used here to evaluate the significance of the random terms. The sd of the intercept for each subject is within (1.459, 2.573) with 95% confidence, whereas (0.114, 0.236) for the sd of slope for week. Since both of the CI’s are above zero, the random terms are kept in the model.
5.2.1.4 Diagnosis plot for residuals.
The following analysis is to evaluate the normality assumption in the error term of the reduced model. Based on the output, we can claim that the normality assumption is valid for the anslysis of the dataset.
6 Summary and Conclusions:
| Resp.name | Slope.test | Reduced.Model | Sim.pval | Int.term.pval | intercept.group | Imp.Topgroup | X2.5.. | X97.5.. | X2.5…1 | X97.5…1 |
|---|---|---|---|---|---|---|---|---|---|---|
| CE | 0.013 | resp ~ group.fa * week + bm + h + (week | id) | 0.124 | 0.035 | -0.920 | 0.180 | 1.107 | 2.195 | 0.115 | 0.251 |
| CF | 0.004 | resp ~ group.fa * week + h + (week | id) | 0.118 | 0.007 | -0.664 | 0.229 | 1.306 | 2.377 | 0.119 | 0.245 |
| LCLF | 0.585 | resp ~ week + bm + h + (week | id) | 0.617 | 0.716 | NA | NA | 1.603 | 2.769 | 0.156 | 0.297 |
| RCLF | 0.133 | resp ~ week + bm + h + (week | id) | 0.459 | 0.096 | NA | NA | 1.459 | 2.573 | 0.114 | 0.236 |
Resp.name: Response variables considered in the modeling.
Slope.test: We first get slopes of linear regression of each individual subject. Then using a two sample t-test for the difference of the slopes in the two groups. it shows what the p-values are for those tests. Based on the output, it’s evident that the CE and CF growth rates are different between the topspin and control group.
REduced.Model: it is the final model evaluated after eliminating non-significant covariates. The group.fa*week include the group effect, weekly effect, and interaction effect of the two. (week|id) considers the longidunal effect of each individaul over the weeks of experiemnt.
Sim.pval: this might be not as important. it baiscallly telling us that by droping those non-significant variables, such as age NC1, NC7, the reduced model is still as powerful as the model including all covairates.
Int.term.pvalue: This is for whether we should drop the interaction effect between group.fa and week. The p-value shows that for both CE and CF. the interaction effect is significant. Which means the topspin group has significanly different growth rates of CE and CF over the time of the experiment.
intercept.group: This shows the difference between the topspin and control group at the beginning of the study. For example, for CE, topspin group is on average 0.92 less than the control group to begin with.
Imp.topgroup: This is ralated with the Int.term.pvalue column. it shows how much MORE the topspin group gain in CE or CF on average than the control group, in each week of the 10-week period in the experiment.
with the intercept.group and Imp.topgroup, we can calculate that after 10 week of experiment, on average, CE from the Topspin group is -0.92+10\(\times\) 0.18=0.875 higher than the CE from the control group;
The Last four columns are for the random effect (week | id) part of the models. It shows the 95% CIs for the intercept and slope of the random effect terms. Since each of the pairs of the CI’s are above zero, we can conclude that the random effect terms are significant and should be included in the models.