W12 in-class exercise2：cochlear implant eyetracking data example
Ching-Fang Wu
Mon Dec 07 07:34:36 2020
1 Introduction
Twenty-eight postlingually-deafened cochlear implant (CI) users and 26 normal hearing (NH) adults perform a spoken-word recognition task. On each trial, 4 pictures appear on a monitor. The subject hears the name of one of them (e.g., wizard as the target) and clicks on the referent. The competitors, in this case, are: lizard (rhyme), whistle (cohort), baggage (unrelated).
Source: Farris-Trimble, A., McMurray, B., Cigrand, N., & Tomblin, J.B. (2014). The process of spoken word recognition in the face of signal degradation: Cochlear implant users and normal-hearing listeners. Journal of Experimental Psychology: Human Perception and Performance, 40(1), 308-327.
Data: ci{bdots}
2 Data Management
# set options
options(digits = 4, show.signif.stars = FALSE)# load packages for later use
pacman::p_load(bdots, tidyverse, nlme)# input data
data(ci)head(ci)## protocol Subject Time Fixations LookType
## 1 NH 36 0 0 Cohort
## 2 NH 36 4 0 Cohort
## 3 NH 36 8 0 Cohort
## 4 NH 36 12 0 Cohort
## 5 NH 36 16 0 Cohort
## 6 NH 36 20 0 Cohortstr(ci)## 'data.frame': 108216 obs. of 5 variables:
## $ protocol : Factor w/ 2 levels "CI","NH": 2 2 2 2 2 2 2 2 2 2 ...
## $ Subject : int 36 36 36 36 36 36 36 36 36 36 ...
## $ Time : int 0 4 8 12 16 20 24 28 32 36 ...
## $ Fixations: num 0 0 0 0 0 0 0 0 0 0 ...
## $ LookType : Factor w/ 4 levels "Cohort","Rhyme",..: 1 1 1 1 1 1 1 1 1 1 ...# set plot themes to black-n-white
ot <- theme_set(theme_bw())#ci$Subject <- as.factor(ci$Subject)
ggplot(ci, aes(Time, Fixations)) +
geom_point(size=rel(.1), pch = 1, alpha = .2) +
facet_grid(protocol ~ LookType)
不同的Looktype。比較CI與NH兩組的Fixations，pattern都差不多。 比較四種不同LookType，CI組的Rhyme和Target的Fixations變異比NH大。
3 plot
# select target looking time course data
dta_t <- subset(ci, LookType == "Target") #只擷取LookType其中一個level作為subset# order data by subject , time
dta_t <- dta_t[order(dta_t$Subject, dta_t$Time, decreasing = F), ]# have a look
head(dta_t)## protocol Subject Time Fixations LookType
## 39079 CI 2 0 0 Target
## 39080 CI 2 4 0 Target
## 39081 CI 2 8 0 Target
## 39082 CI 2 12 0 Target
## 39083 CI 2 16 0 Target
## 39084 CI 2 20 0 Target# convert to groupedData object for nlme processing
dtag <- groupedData(Fixations ~ Time | Subject, outer = ~ protocol,data = as.data.frame(dta_t),
labels = list(x = "Time", y = "Fixations"),
units = list(x = "(ms)", y = "(proportion)") )# lattice plot by subject
plot(dtag, pch = 20, aspect = .89)
# plot by subject given group = protocol
p <- plot(dtag, outer = TRUE, cex = .1, layout = c(2,1), aspect = .89)
update(p, legend = NULL)
分開畫每一個人的curve來看，大部分的人都滿相似的，有部分如68,64,24的curve不太相同。
把所有人的curve畫一起，並分CI與NH來看，兩組曲線圖變異程度不同，CI組個體間的curve圖看起來比較鬆散，NH比較緊密
4 Function form
##
# define the 3-parameter curve fitting function
fx <- function(x, P, S, C) {
P / (1 + exp(S*(C - x)/P))
}5 NLME
5.1 m0
#
# fit nls separately for all subjects
#
m0 <- nlsList(Fixations ~ fx(Time, P, S, C),
start = c(P = 0.87, S = 0.004, C = 757.8),
data = dtag)# results of parameter estimates
coef(m0)## P S C
## 24 0.5571 0.004108 708.9
## 81 0.7501 0.004582 804.5
## 6 0.7480 0.004294 801.2
## 90 0.7795 0.004376 780.4
## 2 0.7770 0.003588 811.9
## 41 0.8007 0.006199 728.2
## 84 0.8190 0.004794 792.4
## 30 0.8259 0.005005 819.4
## 16 0.8409 0.004715 804.8
## 10 0.8438 0.004999 788.8
## 23 0.8913 0.005303 924.2
## 75 0.8929 0.005409 805.3
## 79 0.9074 0.006026 698.8
## 48 0.8881 0.004349 792.9
## 76 0.8908 0.005846 761.5
## 4 0.9082 0.006417 736.7
## 74 0.9020 0.004945 762.7
## 72 0.8887 0.007259 822.3
## 17 0.9086 0.005659 788.8
## 29 0.9401 0.006854 777.9
## 19 0.9224 0.007384 709.1
## 82 0.9483 0.007750 734.1
## 83 0.9332 0.009380 761.1
## 94 0.9353 0.008019 735.2
## 22 0.9474 0.005805 814.3
## 86 0.9674 0.007940 690.5
## 99 0.9681 0.007187 817.4
## 3 0.9716 0.007568 677.7
## 68 0.6759 0.004333 745.9
## 64 0.8652 0.007813 761.6
## 55 0.8759 0.007716 711.8
## 98 0.8861 0.005607 744.8
## 100 0.8969 0.007230 795.8
## 51 0.9269 0.006000 705.6
## 105 0.9254 0.007172 685.0
## 38 0.9375 0.008724 696.7
## 40 0.9477 0.007433 721.8
## 65 0.9450 0.009253 671.1
## 42 0.9622 0.008466 622.0
## 71 0.9539 0.008714 658.7
## 63 0.9633 0.007947 671.7
## 95 0.9625 0.007309 725.4
## 66 0.9669 0.008559 777.7
## 106 0.9630 0.005560 743.0
## 56 0.9674 0.007671 681.7
## 60 0.9673 0.008303 680.0
## 85 0.9713 0.007033 690.3
## 102 0.9636 0.008226 722.6
## 36 0.9795 0.007957 698.8
## 61 0.9684 0.007889 683.7
## 73 0.9773 0.008403 673.0
## 47 0.9812 0.010171 704.2
## 53 0.9811 0.009324 722.4
## 58 0.9768 0.007508 682.55.2 m1
# covert fitting separately by subject to subjects as random effects
m1 <- nlme(m0)## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: false
## convergence (8)## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 2, LME step: nlminb() did not converge (code = 1). PORT message: false
## convergence (8)出現錯誤訊息： Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)Iteration 2, LME step: nlminb() did not converge (code = 1). PORT message: false convergence (8)
# extract mean parameters values to use as initial values
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)
prm0## P S C
## 0.903 0.007 741.2716 m2
#
# update to let parameters dependent on protocol
# and simplify the covariance structure to a simple separate
# variances for each one of the parameters
#
m2 <- update(m1,
fixed = list(P + S + C ~ protocol), #let parameters dependent on protocol
random = pdDiag(P + S + C ~ 1),
groups = ~ Subject,
start = list(fixed = c(P = c(prm0[1], 0),
S = c(prm0[2], 0),
C = c(prm0[3], 0))))
summary(m2)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Fixations ~ fx(Time, P, S, C)
## Data: dtag
## AIC BIC logLik
## -142249 -142167 71134
##
## Random effects:
## Formula: list(P ~ 1, S ~ 1, C ~ 1)
## Level: Subject
## Structure: Diagonal
## P.(Intercept) S.(Intercept) C.(Intercept) Residual
## StdDev: 0.07686 0.001344 45.51 0.01707
##
## Fixed effects: list(P + S + C ~ protocol)
## Value Std.Error DF t-value p-value
## P.(Intercept) 0.9 0.015 26995 59.87 0.0000
## P.protocolNH 0.1 0.021 26995 3.26 0.0011
## S.(Intercept) 0.0 0.000 26995 23.30 0.0000
## S.protocolNH 0.0 0.000 26995 4.87 0.0000
## C.(Intercept) 773.2 8.606 26995 89.85 0.0000
## C.protocolNH -66.4 12.401 26995 -5.35 0.0000
## Correlation:
## P.(In) P.prNH S.(In) S.prNH C.(In)
## P.protocolNH -0.694
## S.(Intercept) 0.000 0.000
## S.protocolNH 0.000 0.000 -0.694
## C.(Intercept) 0.000 0.000 0.000 0.000
## C.protocolNH 0.000 0.000 0.000 0.000 -0.694
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.2733 -0.7499 -0.1167 0.3908 4.9355
##
## Number of Observations: 27054
## Number of Groups: 547 m3 update model 2
#
# update model 2 to include auto-regressive dependence over time
# results have convergence issues - takes a few minutes
# set lmeControl options - change optimizer, etc.
#
m3 <- update(m2, corr = corAR1(0.95, form = ~ 1 | Subject),
control = lmeControl(opt = 'optim',
niterEM = 500,
maxIter = 1e6,
msMaxIter = 1e3))## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value
## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value
## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value
## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value
## Warning in nlm(f = function(nlmePars) -logLik(nlmeSt, nlmePars), p =
## c(coef(nlmeSt)), : NA/Inf replaced by maximum positive value## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 1, LME step: nlm() did not converge (code = 4). Do increase
## 'msMaxIter'!## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 2, LME step: nlm() did not converge (code = 4). Do increase
## 'msMaxIter'!## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 4, LME step: nlm() did not converge (code = 3).## Warning in (function (model, data = sys.frame(sys.parent()), fixed, random, :
## Iteration 5, LME step: nlm() did not converge (code = 3).NA/Inf replaced by maximum positive valueNA/Inf replaced by maximum positive valueNA/Inf replaced by maximum positive valueNA/Inf replaced by maximum positive valueNA/Inf replaced by maximum positive valueIteration 1, LME step: nlm() did not converge (code = 4). Do increase ‘msMaxIter’!Iteration 2, LME step: nlm() did not converge (code = 4). Do increase ‘msMaxIter’!Iteration 4, LME step: nlm() did not converge (code = 3). Iteration 5, LME step: nlm() did not converge (code = 3).
# show summary
summary(m3)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: Fixations ~ fx(Time, P, S, C)
## Data: dtag
## AIC BIC logLik
## -244799 -244709 122411
##
## Random effects:
## Formula: list(P ~ 1, S ~ 1, C ~ 1)
## Level: Subject
## Structure: Diagonal
## P.(Intercept) S.(Intercept) C.(Intercept) Residual
## StdDev: 0.07106 0.001399 46.43 0.02976
##
## Correlation Structure: AR(1)
## Formula: ~1 | Subject
## Parameter estimate(s):
## Phi
## 0.9962
## Fixed effects: list(P + S + C ~ protocol)
## Value Std.Error DF t-value p-value
## P.(Intercept) 0.9 0.014 26995 61.37 0.0000
## P.protocolNH 0.1 0.021 26995 3.17 0.0015
## S.(Intercept) 0.0 0.000 26995 21.23 0.0000
## S.protocolNH 0.0 0.000 26995 4.49 0.0000
## C.(Intercept) 754.6 9.312 26995 81.04 0.0000
## C.protocolNH -69.0 13.205 26995 -5.23 0.0000
## Correlation:
## P.(In) P.prNH S.(In) S.prNH C.(In)
## P.protocolNH -0.696
## S.(Intercept) 0.007 -0.005
## S.protocolNH -0.005 0.008 -0.692
## C.(Intercept) 0.012 -0.008 0.001 0.000
## C.protocolNH -0.008 0.008 -0.001 0.001 -0.705
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.30299 -1.01280 -0.47437 -0.06979 1.58923
##
## Number of Observations: 27054
## Number of Groups: 54#
# check the auto-correlation function of residuals
#
plot(ACF(m3, resType = "n"), alpha = 0.05)
# residual plot to identify outlier
plot(m3, id = 0.05, pch = '.', adj = -0.1, cex = .4)
# residual plot by protocol
plot(m3, resid(., type = "p") ~ fitted(.) | protocol, abline = 0, pch = '.')
# residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject, abline = 0, pch = '.')
# residual plot by subject
plot(m3, resid(., type = "p") ~ fitted(.) | Subject*protocol,
abline = 0, pch = '.')
# plot fitted lines over observations by subject
plot(augPred(m3, level = 0:1))
# fortify data with fitted values
dta_t$yhat <- fitted(m3)# Normal group
ggplot(subset(dta_t, protocol == "NH"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "skyblue") +
geom_line(aes(y = yhat, group = Subject)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ Subject, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")
# implant group
ggplot(subset(dta_t, protocol == "CI"), aes(Time, Fixations)) +
geom_point(size = rel(.5), pch = 1, col = "pink") +
geom_line(aes(y = yhat, group = Subject)) +
coord_cartesian(ylim = c(0, 1)) +
facet_wrap(~ Subject, ncol = 7) +
labs(x = "Time (ms)", y = "Proportion of fixations") +
theme(legend.position = "NONE")