Appendix C: Data, documentation and analysis scripts for Accuracy Assessment of Two Electromagnetic Articulographs: NDI WAVE and NDI VOX
Teja Rebernik (t.rebernik@rug.nl), Jidde Jacobi, Mark Tiede, Martijn Wieling
Generated on 28/01/2021 10:59:23
1 Abstract
Purpose: This study compares two electromagnetic articulographs (EMA) manufactured by Northern Digital, Inc.: the NDI Wave System (2008) and the NDI Vox-EMA System (2020).
Method: Four experiments were completed: (a) comparison of statically positioned sensors; (b) tracking dynamic movements of sensors manipulated using a motor-driven LEGO apparatus; (c) tracking small and large movements of sensors mounted in a rigid bar manipulated by hand; and (d) tracking movements of sensors rotated on a circular disc. We assessed spatial variability for statically positioned sensors, variability in the transduced Euclidean distances (EDs) between sensor pairs, and missing data rates. For sensors tracking circular movements, we compared the fit between fitted ideal circles and actual trajectories.
Results: The average sensor pair tracking error (i.e., the standard deviation of the EDs) was 1.37 mm for the WAVE and 0.12 mm for the VOX during automated trials at the fastest speed, and 0.35 mm for the WAVE and 0.14 mm for the VOX during the tracking of large manual movements. The average standard deviation of the fitted circle radii charted by manual circular disc movements was 0.72 mm for the WAVE sensors and 0.14 mm for the VOX sensors. There was no significant difference between the WAVE and the VOX in the number of missing frames.
Conclusions: In general, the VOX system significantly outperformed the WAVE on measures of both static precision and dynamic accuracy (automated and manual). For both systems, positional precision and spatial variability were influenced by the sensors’ position relative to the field generator unit (FGU; worse when further away).
2 Libraries and functions
# load required packages
library(mgcv)## Loading required package: nlme## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.packageVersion("mgcv")## [1] '1.8.33'library(itsadug)## Loading required package: plotfunctions## Loaded package itsadug 2.4 (see 'help("itsadug")' ).packageVersion("itsadug")## [1] '2.4'library(beanplot)
packageVersion("beanplot")## [1] '1.2'# R version information
R.version.string## [1] "R version 4.0.2 (2020-06-22)"3 Static trials
3.1 Dataset
tab2 <- read.csv("static-coordinates.csv")
tab2$articulograph = as.factor(tab2$articulograph)
tab2$position = as.factor(tab2$position)
tab2$trial = as.factor(tab2$trial)
tab2$sensorpos = as.factor(tab2$sensor)
tab2$sensorset = (tab2$trial %in% c('t1','t2'))
tab2$sensorset = paste(tab2$sensorpos,tab2$sensorset,sep='')
tab2$sensorset = as.factor(tab2$sensorset)
tab2$sensor = interaction(tab2$articulograph,tab2$sensorpos)
tab2$coordinate = as.factor(tab2$coordinate)
levels(tab2$coordinate) = c('Pitch', 'X','Y','Z','Yaw')
tab2 = droplevels(tab2[tab2$position != 'p16',]) # remove p16 as it was not used in our data analysisstr(tab2)## 'data.frame': 4800 obs. of 8 variables:
## $ articulograph: Factor w/ 2 levels "VOX","WAVE": 2 2 2 2 2 2 2 2 1 1 ...
## $ position : Factor w/ 15 levels "p1","p10","p11",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ trial : Factor w/ 4 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sensor : Factor w/ 16 levels "VOX.s1","WAVE.s1",..: 2 4 6 8 10 12 14 16 1 3 ...
## $ coordinate : Factor w/ 5 levels "Pitch","X","Y",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ SD : num 0.0151 0.0261 0.0104 0.0051 0.0084 0.0133 0.0604 0.0374 0.0345 0.0031 ...
## $ sensorpos : Factor w/ 8 levels "s1","s2","s3",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ sensorset : Factor w/ 16 levels "s1FALSE","s1TRUE",..: 2 4 6 8 10 12 14 16 2 4 ...head(tab2)## articulograph position trial sensor coordinate SD sensorpos sensorset
## 1 WAVE p1 t1 WAVE.s1 X 0.0151 s1 s1TRUE
## 2 WAVE p1 t1 WAVE.s2 X 0.0261 s2 s2TRUE
## 3 WAVE p1 t1 WAVE.s3 X 0.0104 s3 s3TRUE
## 4 WAVE p1 t1 WAVE.s4 X 0.0051 s4 s4TRUE
## 5 WAVE p1 t1 WAVE.s5 X 0.0084 s5 s5TRUE
## 6 WAVE p1 t1 WAVE.s6 X 0.0133 s6 s6TRUE3.2 Trial variability
We compared the four trials to see if there is a difference between them. All trials show similar patterns, which led us to averaging our results in the tables found in the paper.
# there is no differential effect of trial
summary(m<-gam(SD+1e-04 ~ trial*articulograph+ s(sensorset,articulograph,bs='re'), data = tab2[tab2$coordinate %in% c('X','Y','Z','Pitch','Yaw'),],method='REML',family=Gamma(link=log))) ##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ trial * articulograph + s(sensorset, articulograph,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.512776 0.477279 -7.360 2.15e-13 ***
## trialt2 -0.013647 0.136130 -0.100 0.9202
## trialt3 -0.225549 0.674974 -0.334 0.7383
## trialt4 -0.228505 0.674974 -0.339 0.7350
## articulographWAVE 1.535972 0.674974 2.276 0.0229 *
## trialt2:articulographWAVE 0.003894 0.192517 0.020 0.9839
## trialt3:articulographWAVE 0.388446 0.954558 0.407 0.6841
## trialt4:articulographWAVE 0.411314 0.954558 0.431 0.6666
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sensorset,articulograph) 27.42 28 16.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.00906 Deviance explained = 43.3%
## -REML = -9293.9 Scale est. = 5.5594 n = 4800Which is also clear from the visualization:
par(mfrow=c(2,2)) # note that for the visualization a small value is added as taking the logarithm of 0 is not possible
beanplot(log(SD+1e-04)~articulograph,data=tab2[tab2$trial=='t1',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='Trial 1',ylim=c(-10,6))
beanplot(log(SD+1e-04)~articulograph,data=tab2[tab2$trial=='t2',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='Trial 2',ylim=c(-10,6))
beanplot(log(SD+1e-04)~articulograph,data=tab2[tab2$trial=='t3',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='Trial 3',ylim=c(-10,6))
beanplot(log(SD+1e-04)~articulograph,data=tab2[tab2$trial=='t4',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='Trial 4',ylim=c(-10,6))
3.3 Static precision comparison
This part assesses the static precision of the WAVE and the VOX. It corresponds to Table 2 and Figure 10 in the paper.
3.3.1 GAM model
We assessed the static precision of the VOX and the WAVE. The VOX was significantly better (p < 0.001)
# Gamma model due to non-Gaussian distribution of residuals (due to large outliers of Wave)
summary(m<-gam(SD+1e-04 ~ articulograph+ s(sensorset,articulograph,bs='re'), data = tab2[tab2$coordinate %in% c('X','Y','Z','Pitch','Yaw'),],method='REML',family=Gamma(link=log))) ##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ articulograph + s(sensorset, articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.6297 0.3262 -11.129 < 2e-16 ***
## articulographWAVE 1.7369 0.4613 3.766 0.000168 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sensorset,articulograph) 29.35 30 15.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.00989 Deviance explained = 43.3%
## -REML = -9300 Scale est. = 5.5724 n = 4800# Quantification of the effect on the measurement scale
summary(m<-gam(SD ~ articulograph + s(sensorset,articulograph,bs='re'), data = tab2[tab2$coordinate %in% c('X','Y','Z','Pitch','Yaw')&tab2$SD < 100,],method='REML')) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## SD ~ articulograph + s(sensorset, articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.03277 0.01499 2.186 0.0289 *
## articulographWAVE 0.11541 0.02120 5.443 5.51e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sensorset,articulograph) 22.69 30 3.095 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0426 Deviance explained = 4.73%
## -REML = 1967.2 Scale est. = 0.13139 n = 47983.3.2 Visualization
# Figure 10 in paper
par(mfrow=c(1,1))
beanplot(log(SD+1e-04)~articulograph,data=tab2[tab2$coordinate %in% c('X','Y','Z','Pitch','Yaw'),],bw='bcv',main='',ylab='SD (log-scale)',log='',method='jitter')
3.4 Effect of FGU distance
This part assesses the effect of actual distance from FGU (along the z-axis) on the SDs.
# Defining the distance from the FGU for each position (movement along the z-axis)
tab2$distance = NA
tab2[tab2$position=='p1',]$distance = 9
tab2[tab2$position=='p2',]$distance = 13
tab2[tab2$position=='p3',]$distance = 17.8
tab2[tab2$position=='p4',]$distance = 22.6
tab2[tab2$position=='p5',]$distance = 27.4
tab2[tab2$position=='p6',]$distance = 33.8
tab2p = droplevels(tab2[tab2$position %in% c('p1','p2','p3','p4','p5','p6'),])3.4.1 GAM model
We assessed the effect of distance on the static precision. Both articulographs were less accurate when sensors were further away from the FGU.
summary(m<-bam(SD+1e-4 ~ s(distance,by=articulograph,k=6) + articulograph + s(distance,sensorset,by=articulograph, bs='fs',m=1,k=6), data = tab2p,method='fREML',family=Gamma(link=log), discrete=T)) ##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ s(distance, by = articulograph, k = 6) + articulograph +
## s(distance, sensorset, by = articulograph, bs = "fs", m = 1,
## k = 6)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.7112 0.2152 -17.243 < 2e-16 ***
## articulographWAVE 1.3561 0.2216 6.119 1.15e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(distance):articulographVOX 2.855 3.3 52.506 < 2e-16 ***
## s(distance):articulographWAVE 1.000 1.0 1286.060 < 2e-16 ***
## s(distance,sensorset):articulographVOX 45.850 94.0 6.656 < 2e-16 ***
## s(distance,sensorset):articulographWAVE 9.521 94.0 0.277 0.000348 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.538 Deviance explained = 63.6%
## fREML = 2784.4 Scale est. = 0.98052 n = 19203.4.2 GAM model: WAVE vs. VOX
We additionally assessed statistically whether the WAVE is more affected by distance than the VOX. The results show the WAVE is affected more than the VOX (significantly). This is visualized in Figure 11 in the paper.
tab2p$articulographO = as.ordered(tab2p$articulograph)
contrasts(tab2p$articulographO) = 'contr.treatment'
summary(m<-bam(SD+1e-4 ~ s(distance,k=6) + s(distance,by=articulographO,k=6) + articulographO + s(distance,sensorset,by=articulograph, bs='fs',m=1,k=6), data = tab2p,method='fREML',family=Gamma(link=log), discrete=T))##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ s(distance, k = 6) + s(distance, by = articulographO,
## k = 6) + articulographO + s(distance, sensorset, by = articulograph,
## bs = "fs", m = 1, k = 6)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.6774 0.2143 -17.159 < 2e-16 ***
## articulographOWAVE 1.3455 0.2208 6.093 1.34e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(distance) 3.187 3.757 46.561 < 2e-16 ***
## s(distance):articulographOWAVE 1.000 1.000 11.619 0.000667 ***
## s(distance,sensorset):articulographVOX 46.310 94.000 6.701 < 2e-16 ***
## s(distance,sensorset):articulographWAVE 9.593 94.000 0.283 0.000289 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.539 Deviance explained = 63.6%
## fREML = 2784 Scale est. = 0.97754 n = 1920# As s(distance):articulographOWAVE has 1 edf, it is linear, so we fit a linear interaction directly to assess the direction of the effect
summary(m<-gam(SD+1e-4 ~ articulograph*distance + s(sensorset,articulograph,bs='re') + s(sensorset,distance,bs='re'), data = tab2p,method='REML',family=Gamma(link=log))) ##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ articulograph * distance + s(sensorset, articulograph,
## bs = "re") + s(sensorset, distance, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.815512 0.270890 -21.468 < 2e-16 ***
## articulographWAVE 0.609554 0.383096 1.591 0.112
## distance 0.103157 0.007066 14.600 < 2e-16 ***
## articulographWAVE:distance 0.033707 0.006060 5.562 3.05e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sensorset,articulograph) 28.17 30 520.4 < 2e-16 ***
## s(sensorset,distance) 11.90 15 594.2 7.06e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.528 Deviance explained = 62%
## -REML = -3818.3 Scale est. = 1.2477 n = 19203.4.3 Visualization
# Figure 11 in paper
tab2p = droplevels(tab2[tab2$position %in% c('p1','p2','p3','p4','p5','p6'),])
levels(tab2p$position) = c('1','2','3','4','5','6')
summary(m<-bam(SD+1e-4 ~ s(distance,by=articulograph,k=6) + articulograph + s(distance,sensorset,by=articulograph, bs='fs',m=1,k=6), data = tab2p,method='fREML',family=Gamma(link=log), discrete=T))##
## Family: Gamma
## Link function: log
##
## Formula:
## SD + 1e-04 ~ s(distance, by = articulograph, k = 6) + articulograph +
## s(distance, sensorset, by = articulograph, bs = "fs", m = 1,
## k = 6)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.7112 0.2152 -17.243 < 2e-16 ***
## articulographWAVE 1.3561 0.2216 6.119 1.15e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(distance):articulographVOX 2.855 3.3 52.506 < 2e-16 ***
## s(distance):articulographWAVE 1.000 1.0 1286.060 < 2e-16 ***
## s(distance,sensorset):articulographVOX 45.850 94.0 6.656 < 2e-16 ***
## s(distance,sensorset):articulographWAVE 9.521 94.0 0.277 0.000348 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.538 Deviance explained = 63.6%
## fREML = 2784.4 Scale est. = 0.98052 n = 1920par(mfrow=c(1,3))
plot_smooth(m,view='distance',plot_all='articulograph', shade=T,ylab='SD (log-scale)',xlab='Distance from FGU',rug=F, col=c('black','darkgray'), rm.ranef=T)## Summary:
## * articulograph : factor; set to the value(s): VOX, WAVE.
## * distance : numeric predictor; with 30 values ranging from 9.000000 to 33.800000.
## * sensorset : factor; set to the value(s): s1FALSE. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(distance,sensorset):articulographVOX,s(distance,sensorset):articulographWAVE
## beanplot(log(SD+1e-04)~position,data=tab2p[tab2p$articulograph=='VOX',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='VOX',ylim=c(-10,2),xlab='Position')
beanplot(log(SD+1e-04)~position,data=tab2p[tab2p$articulograph=='WAVE',],bw='bcv',ylab='SD (log-scale)',log='',method='jitter',main='WAVE',ylim=c(-10,2),xlab='Position')
3.5 Distance measurements
This part assesses the error between the actual (hand-measured) distance and the distance determined using the reported sensor coordinates. It corresponds to Table 3 and Figure 12 in the paper. The VOX and the WAVE do not significantly differ regarding the inferred distances (p = 0.5)
3.5.1 Dataset
tab3 <- read.csv("distance-measures.csv")
tab3$articulograph = as.factor(tab3$articulograph)
tab3$trial = as.factor(tab3$trial)
tab3$sensor = as.factor(tab3$sensor)
tab3$sensorSet = as.factor(tab3$sensorSet)
tab3$sens = interaction(tab3$sensor,tab3$sensorSet)
tab3$positionPair = as.factor(tab3$positionPair)str(tab3)## 'data.frame': 176 obs. of 7 variables:
## $ articulograph: Factor w/ 2 levels "VOX","WAVE": 2 2 2 2 2 2 2 2 2 2 ...
## $ trial : Factor w/ 4 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sensor : Factor w/ 2 levels "s1","s2": 1 1 1 1 1 1 1 1 1 1 ...
## $ sensorSet : Factor w/ 4 levels "sv1","sv2","sw1",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ positionPair : Factor w/ 11 levels "eightnine","eleventwelve",..: 6 11 10 4 3 7 1 8 2 9 ...
## $ error : num 0.065 0.279 0.521 0.732 1.158 ...
## $ sens : Factor w/ 8 levels "s1.sv1","s2.sv1",..: 5 5 5 5 5 5 5 5 5 5 ...head(tab3)## articulograph trial sensor sensorSet positionPair error sens
## 1 WAVE t1 s1 sw1 onetwo 0.0650 s1.sw1
## 2 WAVE t1 s1 sw1 twothree 0.2786 s1.sw1
## 3 WAVE t1 s1 sw1 threefour 0.5211 s1.sw1
## 4 WAVE t1 s1 sw1 fourfive 0.7316 s1.sw1
## 5 WAVE t1 s1 sw1 fivesix 1.1579 s1.sw1
## 6 WAVE t1 s1 sw1 seveneight 1.0948 s1.sw13.5.2 GAM model
# Gamma-distributed model with log-link: not significant
summary(m<-gam(error ~ articulograph + s(sens,articulograph,bs='re') + s(positionPair,articulograph,bs='re'), data = tab3,method='REML',family=Gamma(link=log)))##
## Family: Gamma
## Link function: log
##
## Formula:
## error ~ articulograph + s(sens, articulograph, bs = "re") + s(positionPair,
## articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7289 0.2701 -6.401 1.84e-09 ***
## articulographWAVE 0.2630 0.3820 0.689 0.492
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sens,articulograph) 4.246 6 3.549 0.000103 ***
## s(positionPair,articulograph) 18.844 20 24.619 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.709 Deviance explained = 70.7%
## -REML = -85.986 Scale est. = 0.32557 n = 176# Gaussian model (using the original measurement scale): also not significant
summary(m<-gam(error ~ articulograph + s(sens,articulograph,bs='re') + s(positionPair,articulograph,bs='re'), data = tab3,method='REML')) ##
## Family: gaussian
## Link function: identity
##
## Formula:
## error ~ articulograph + s(sens, articulograph, bs = "re") + s(positionPair,
## articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2332 0.1426 1.635 0.104
## articulographWAVE 0.2391 0.2017 1.185 0.238
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sens,articulograph) 4.984 6 4.907 1.51e-05 ***
## s(positionPair,articulograph) 18.923 20 17.577 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.702 Deviance explained = 74.5%
## -REML = 62.003 Scale est. = 0.076248 n = 1763.5.3 Visualization
# Figure 12 in paper
par(mfrow=c(1,2)); beanplot(log(error)~articulograph,data=tab3,bw='bcv',main='',ylab='error (log-scale)',log='')
beanplot(error~articulograph,data=tab3,bw='bcv',main='',ylab='error (mm)',log='')
4 Dynamic trials
This part assesses dynamic accuracy during automated dynamic trials. It corresponds to Table 4 and Figures 17 and 18 in the paper.
4.1 Dataset
datAut <- read.csv("vox-wave-automatic.csv")
datAut$articulograph <- as.factor(datAut$articulograph)
levels(datAut$articulograph) = c('VOX','WAVE') # capitalize levels
datAut$rigidBody <- as.factor(datAut$rigidBody)
datAut$trial <- as.factor(datAut$trial)
datAut$sensorSet <- as.factor(datAut$sensorSet)str(datAut)## 'data.frame': 160 obs. of 6 variables:
## $ articulograph: Factor w/ 2 levels "VOX","WAVE": 2 2 2 2 2 2 2 2 2 2 ...
## $ rigidBody : Factor w/ 5 levels "circle","eccentricLarge",..: 1 1 1 1 3 3 3 3 2 2 ...
## $ speed : int 0 1 2 3 0 1 2 3 0 1 ...
## $ trial : Factor w/ 4 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sensorSet : Factor w/ 4 levels "v1","v2","w1",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ SD : num 0.76 3.74 5.16 6.27 0.09 0.23 0.33 0.44 0.03 0.09 ...head(datAut)## articulograph rigidBody speed trial sensorSet SD
## 1 WAVE circle 0 t1 w1 0.76
## 2 WAVE circle 1 t1 w1 3.74
## 3 WAVE circle 2 t1 w1 5.16
## 4 WAVE circle 3 t1 w1 6.27
## 5 WAVE eccentricSmall 0 t1 w1 0.09
## 6 WAVE eccentricSmall 1 t1 w1 0.234.2 GAM model
Here we assess the articulographs’ accuracy during dynamic trials. Only automated data is used, as the manual dynamic movements are likely not consistent across trials and articulographs. The circularly rotating bar showed the largest standard deviations.
summary(gam(SD ~ articulograph + s(rigidBody,articulograph,bs='re'), data = datAut,family=Gamma(link=log), method='REML'))##
## Family: Gamma
## Link function: log
##
## Formula:
## SD ~ articulograph + s(rigidBody, articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.9420 0.4858 -6.056 1.07e-08 ***
## articulographWAVE 1.6603 0.6870 2.417 0.0169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(rigidBody,articulograph) 7.856 8 74.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.702 Deviance explained = 85.2%
## -REML = -179.75 Scale est. = 0.33917 n = 160# Quantification of the effect on the measurement scale: non-circle
summary(gam(SD ~ articulograph + s(rigidBody,articulograph,bs='re'), data = datAut[datAut$rigidBody != 'circle',], method='REML')) # normal scale: 0.1 mm difference##
## Family: gaussian
## Link function: identity
##
## Formula:
## SD ~ articulograph + s(rigidBody, articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.04469 0.02487 1.797 0.07486 .
## articulographWAVE 0.11078 0.03517 3.150 0.00206 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(rigidBody,articulograph) 5.225 6 6.738 8.18e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.478 Deviance explained = 50.3%
## -REML = -143.27 Scale est. = 0.0051156 n = 128# Quantification of the effect on the measurement scale: circle
summary(gam(SD ~ articulograph, data = datAut[datAut$rigidBody == 'circle',], method='REML')) # normal scale for the circle: 3.4 mm difference##
## Family: gaussian
## Link function: identity
##
## Formula:
## SD ~ articulograph
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2025 0.3789 0.534 0.597
## articulographWAVE 3.4256 0.5359 6.392 4.68e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## R-sq.(adj) = 0.563 Deviance explained = 57.7%
## -REML = 57.819 Scale est. = 2.2976 n = 32# Figure 17 in paper
par(mfrow=c(2,3))
beanplot(log(SD)~articulograph,data=droplevels(datAut),bw='bcv',main='All rigid bodies',ylab='SD (log-scale)',log='',method='jitter',ylim=c(-7,3.5))
beanplot(log(SD)~articulograph,data=droplevels(datAut[datAut$rigidBody=='static',]),bw='bcv',main='Static',ylab='SD (log-scale)',ylim=c(-7,3.5),method='jitter',log='')
beanplot(log(SD)~articulograph,data=droplevels(datAut[datAut$rigidBody=='eccentricSmall',]),bw='bcv',main='Eccentric (small)',ylab='SD (log-scale)',ylim=c(-7,3.5),method='jitter',log='')
beanplot(log(SD)~articulograph,data=droplevels(datAut[datAut$rigidBody=='eccentricLarge',]),bw='bcv',main='Eccentric (large)',ylab='SD (log-scale)',ylim=c(-7,3.5),method='jitter',log='')
beanplot(log(SD)~articulograph,data=droplevels(datAut[datAut$rigidBody=='piston',]),bw='bcv',main='Piston',ylab='SD (log-scale)',ylim=c(-7,3.5),method='jitter',log='')
beanplot(log(SD)~articulograph,data=droplevels(datAut[datAut$rigidBody=='circle',]),bw='bcv',main='Circle',ylab='SD (log-scale)',ylim=c(-7,3.5),method='jitter',log='')
4.3 Effect of speed on accuracy
4.3.1 GAM model
summary(m<-gam(SD ~ s(speed,by=articulograph,k=4) + articulograph + s(speed,rigidBody,by=articulograph,bs='fs',m=1,k=4), data = datAut,family=Gamma(link=log), method='REML'))##
## Family: Gamma
## Link function: log
##
## Formula:
## SD ~ s(speed, by = articulograph, k = 4) + articulograph + s(speed,
## rigidBody, by = articulograph, bs = "fs", m = 1, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.1494 0.3769 -8.356 1.12e-13 ***
## articulographWAVE 1.6062 0.7301 2.200 0.0297 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speed):articulographVOX 2.389 2.480 16.84 3.92e-07 ***
## s(speed):articulographWAVE 2.438 2.487 12.09 2.92e-05 ***
## s(speed,rigidBody):articulographVOX 14.304 18.000 60.09 < 2e-16 ***
## s(speed,rigidBody):articulographWAVE 15.312 18.000 140.81 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.968 Deviance explained = 98.5%
## -REML = -306.26 Scale est. = 0.058838 n = 1604.3.2 Visualization
# Figure 18 in paper
par(mfrow=c(1,3))
plot_smooth(m,view='speed',plot_all='articulograph', shade=T,ylab='SD (log-scale)',xlab='Movement speed',rug=F, col=c('black','darkgray'), rm.ranef=T)## Summary:
## * articulograph : factor; set to the value(s): VOX, WAVE.
## * speed : numeric predictor; with 30 values ranging from 0.000000 to 3.000000.
## * rigidBody : factor; set to the value(s): circle. (Might be canceled as random effect, check below.)
## * NOTE : The following random effects columns are canceled: s(speed,rigidBody):articulographVOX,s(speed,rigidBody):articulographWAVE
## beanplot(log(SD)~speed,data=droplevels(datAut[datAut$articulograph=='VOX',]),bw='bcv',main='VOX',log="",ylab='SD (log-scale)',xlab='Speed',ylim=c(-8,4),method='jitter')
beanplot(log(SD)~speed,data=droplevels(datAut[datAut$articulograph=='WAVE',]),bw='bcv',main='WAVE',log="",ylab='SD (log-scale)',xlab='Speed',ylim=c(-8,4),method='jitter')
4.3.3 GAM model: WAVE vs. VOX
datAut$articulographO = as.ordered(datAut$articulograph)
contrasts(datAut$articulographO) = 'contr.treatment'summary(m<-gam(SD ~ s(speed,k=4) + s(speed, by=articulographO,k=4) + articulographO + s(speed,rigidBody,by=articulograph,bs='fs',m=1,k=4), data = datAut[datAut$rigidBody!='circle',], method='REML', family=Gamma(link=log))) #when log-transformed, WAVE is not significantly worse than VOX at faster speeds (p = 0.9)##
## Family: Gamma
## Link function: log
##
## Formula:
## SD ~ s(speed, k = 4) + s(speed, by = articulographO, k = 4) +
## articulographO + s(speed, rigidBody, by = articulograph,
## bs = "fs", m = 1, k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.4679 0.2736 -12.676 < 2e-16 ***
## articulographOWAVE 1.2657 0.3683 3.437 0.000863 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(speed) 2.409 2.483 13.285 9.55e-06 ***
## s(speed):articulographOWAVE 1.000 1.000 0.012 0.913
## s(speed,rigidBody):articulographVOX 11.199 14.000 30.364 < 2e-16 ***
## s(speed,rigidBody):articulographWAVE 12.420 14.000 25.481 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.84 Deviance explained = 95.2%
## -REML = -308.17 Scale est. = 0.062998 n = 1285 Missing data
This part assesses missing data using logistic regression. It corresponds to Table 5 and Figure 19 in the paper.
5.1 Dataset
mis <- read.csv("missing-data.csv")
mis$articulograph = as.factor(mis$articulograph)
mis$trial = as.factor(mis$trial)
mis$sensorSet = as.factor(mis$sensorSet)
mis$dynamicTask = as.factor(mis$dynamicTask)
mis$sensor = as.factor(mis$sensor)
mis$sens = interaction(mis$sensor,mis$sensorSet)
mis$missingFrames = round(8000 * mis$missing / 100)
mis$nonMissingFrames = 8000 - mis$missingFramesstr(mis)## 'data.frame': 360 obs. of 9 variables:
## $ articulograph : Factor w/ 2 levels "VOX","WAVE": 2 2 2 2 2 2 2 2 2 2 ...
## $ trial : Factor w/ 4 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ sensorSet : Factor w/ 8 levels "sv1a","sv1b",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ dynamicTask : Factor w/ 6 levels "automatedFast",..: 3 3 3 3 3 3 3 3 2 2 ...
## $ sensor : Factor w/ 21 levels "bs1","bs2","bs3",..: 14 15 16 17 18 19 20 21 14 15 ...
## $ missing : num 1.6 4.5 0 0 0 0 0 0 6 6 ...
## $ sens : Factor w/ 168 levels "bs1.sv1a","bs2.sv1a",..: 98 99 100 101 102 103 104 105 98 99 ...
## $ missingFrames : num 128 360 0 0 0 0 0 0 480 480 ...
## $ nonMissingFrames: num 7872 7640 8000 8000 8000 ...head(mis)## articulograph trial sensorSet dynamicTask sensor missing sens
## 1 WAVE t1 sw1a automatedSlow sa1 1.6 sa1.sw1a
## 2 WAVE t1 sw1a automatedSlow sa2 4.5 sa2.sw1a
## 3 WAVE t1 sw1a automatedSlow sa3 0.0 sa3.sw1a
## 4 WAVE t1 sw1a automatedSlow sa4 0.0 sa4.sw1a
## 5 WAVE t1 sw1a automatedSlow sa5 0.0 sa5.sw1a
## 6 WAVE t1 sw1a automatedSlow sa6 0.0 sa6.sw1a
## missingFrames nonMissingFrames
## 1 128 7872
## 2 360 7640
## 3 0 8000
## 4 0 8000
## 5 0 8000
## 6 0 80005.2 GAM model: WAVE vs. VOX
summary(m<-gam(cbind(missingFrames,nonMissingFrames) ~ articulograph + s(sens,articulograph,bs='re'), data = mis,family=binomial, method='REML'))##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(missingFrames, nonMissingFrames) ~ articulograph + s(sens,
## articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.3845 0.9245 -12.314 <2e-16 ***
## articulographWAVE 1.7694 1.2564 1.408 0.159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(sens,articulograph) 65.37 82 8325 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.674 Deviance explained = 90.8%
## -REML = 3333.7 Scale est. = 1 n = 3605.2.1 Visualization
# Figure 19 in paper
par(mfrow=c(1,1))
beanplot(missing~articulograph,data=mis,bw='bcv',log="",ylab='Missing (%)', cutmin = 0)
5.2.2 Effect of test type
# Assessing whether 'type' (manual vs. automated) plays a role
mis$type = 'Manual'
mis[mis$dynamicTask %in% c('automatedFast', 'automatedMedium','automatedSlow'),]$type = 'Automatic'
mis$type = as.factor(mis$type)
mis$type = relevel(mis$type,'Manual')summary(m<-gam(cbind(missingFrames,nonMissingFrames) ~ articulograph*type + s(sens,articulograph,bs='re'), data = mis,family=binomial, method='REML'))##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(missingFrames, nonMissingFrames) ~ articulograph * type +
## s(sens, articulograph, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.1861 1.1041 -9.226 <2e-16 ***
## articulographWAVE 0.5238 1.5373 0.341 0.7333
## typeAutomatic -3.9421 2.3196 -1.700 0.0892 .
## articulographWAVE:typeAutomatic 4.0860 2.9008 1.409 0.1590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(sens,articulograph) 60.65 80 3318 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.677 Deviance explained = 90.8%
## -REML = 3325.9 Scale est. = 1 n = 360