Homework_1130_1

Exercise 1

MNREAD Acuity Charts (Optelec, US Inc.) were used to measure reading speed as a function of print size. MNREAD data were collected from 42 adult observers with normal vision. Reading speed was measured in log word per minute. Print size was measured in in logMAR (related to degrees in visual angle).

Source: Cheung, C., Kallie, S., Legge, G.E., & Cheong, A.M.Y. (2008). Nonlinear Mixed-Effects Modeling of MNREAD Data. Investigative Ophthalmology & Visual Science, 49(2), 828-835.

• Column 1: Subject ID
• Column 2: Reading time in logWPM
• Column 3: Print size measured in logMAR

# set digit and stars options
options(digits = 4, show.signif.stars = FALSE)
# load packages for later use
pacman::p_load(tidyverse, nlme)

# regular contrast polarity (i.e., black on white)
dta1 <- read.table("mnread0.txt", header = T)
str(dta1)

## 'data.frame':    670 obs. of  3 variables:
##  $ ID   : chr  "S101" "S101" "S101" "S101" ...
##  $ pSize: num  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 ...
##  $ logRS: num  0.951 2.028 2.092 2.219 2.292 ...

# set plot theme to black-n-white
ot <- theme_set(theme_bw())

# individual profiles on separate panels
ggplot(dta1, aes(pSize, logRS)) +
 geom_point(size = rel(.8), pch = 20)+
 geom_line()+
 facet_wrap( ~ ID) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

#  individual profiles on one panel
ggplot(dta1, aes(pSize, logRS, group = ID)) +
 geom_point(pch = 20)+
 geom_line(alpha = .3) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

# coerce data into groupedData object
##groupedData=formula+Data
dta_g <- groupedData(logRS ~ pSize | ID, data = as.data.frame(dta1),
                     labels = list(x = "Print size", y = "Reading speed"),
                     units = list(x = "(logMAR)", y = "(logWPM)") )

# lattice plot by subject, use groupedData
plot(dta_g, pch = 20, aspect = .89)

Modeling

# try nls for all subjects
##nlsList= List of nls Objects with a Common Model, without random effect
m0_g <- nlsList(SSasympOff, data = dta_g)
# results of parameter estimates
coef(m0_g)

##       Asym   lrc       c0
## S128 2.131 1.590 -0.30739
## S132 2.186 2.388 -0.40812
## S102 2.154 2.350 -0.04337
## S131 2.193 1.551 -0.60031
## S122 2.166 2.270 -0.34818
## S109 2.202 3.196 -0.01809
## S108 2.171 2.353 -0.23285
## S130 2.187 2.332 -0.31713
## S106 2.184 2.372 -0.38127
## S117 2.212 1.819 -0.22989
## S136 2.275 2.769 -0.08696
## S103    NA    NA       NA
## S113 2.240 2.177 -0.21904
## S101 2.251 2.758 -0.03504
## S123 2.274 1.853 -0.30819
## S133 2.225 2.362 -0.44646
## S112 2.270 2.075 -0.23090
## S142 2.258 2.169 -0.33552
## S118 2.267 2.607 -0.32027
## S104 2.236 2.431 -0.24719
## S120 2.268 1.521 -0.31138
## S105 2.289 2.420 -0.25797
## S124 2.290 1.788 -0.16797
## S129 2.258 2.130 -0.32087
## S134 2.290 2.236 -0.37175
## S135 2.287 2.071 -0.29440
## S125 2.272 1.961 -0.28334
## S126 2.317 1.685 -0.47288
## S141 2.323 2.082 -0.29600
## S111 2.329 2.233 -0.42614
## S115 2.306 1.623 -0.56310
## S110 2.288 2.401 -0.14867
## S107 2.321 1.720 -0.37909
## S114 2.334 2.406 -0.24746
## S139    NA    NA       NA
## S137 2.334 2.622 -0.26045
## S140 2.352 2.190 -0.22763
## S138 2.325 2.358 -0.37132
## S127 2.316 1.982 -0.32574
## S119 2.386 1.902 -0.30892
## S116 2.412 1.710 -0.41221
## S121 2.456 2.476 -0.22895

# random effects model equivalent
m1_g <- nlme(m0_g)
# show summary
summary(m1_g)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: logRS ~ SSasympOff(pSize, Asym, lrc, c0) 
##  Data: dta_g 
##     AIC   BIC logLik
##   -1369 -1324  694.7
## 
## Random effects:
##  Formula: list(Asym ~ 1, lrc ~ 1, c0 ~ 1)
##  Level: ID
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev  Corr         
## Asym     0.06717 Asym   lrc   
## lrc      0.33175 -0.032       
## c0       0.10966 -0.103  0.492
## Residual 0.06186              
## 
## Fixed effects: list(Asym ~ 1, lrc ~ 1, c0 ~ 1) 
##        Value Std.Error  DF t-value p-value
## Asym  2.2700   0.01079 626  210.47       0
## lrc   2.2028   0.05630 626   39.13       0
## c0   -0.2856   0.01739 626  -16.43       0
##  Correlation: 
##     Asym   lrc   
## lrc -0.072       
## c0  -0.111  0.514
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -6.62871 -0.39007  0.08947  0.56091  2.61865 
## 
## Number of Observations: 670
## Number of Groups: 42

# fix c0
m1a_g <- update(m1_g, random = Asym + lrc ~ 1)

# model without random effects for the covariance terms involving Asym
m2_g <- update(m1_g, random = pdBlocked(list(Asym ~ 1, lrc + c0 ~ 1)))

Model comparison

anova(m1_g, m2_g, m1a_g) #choose m1_g

##       Model df     AIC     BIC logLik   Test L.Ratio p-value
## m1_g      1 10 -1369.5 -1324.4  694.7                       
## m2_g      2  8 -1373.0 -1337.0  694.5 1 vs 2     0.4  0.8178
## m1a_g     3  7  -551.2  -519.7  282.6 2 vs 3   823.8  <.0001

Plot Fitted values and Residuals

# residuals
plot(m1a_g, id = .05, cex = .6, adj=-.05)

# random effects
plot(ranef(m1a_g, augFrame = T))

# parameter estimates CIs
intervals(m1a_g)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##        lower    est.   upper
## Asym  2.2820  2.3055  2.3289
## lrc   1.4610  1.5918  1.7226
## c0   -0.4176 -0.4096 -0.4016
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: ID 
##                  lower     est.    upper
## sd(Asym)       0.04262  0.06013  0.08484
## sd(lrc)        0.31109  0.39712  0.50693
## cor(Asym,lrc) -0.69702 -0.42048 -0.03506
## 
##  Within-group standard error:
##  lower   est.  upper 
## 0.1336 0.1415 0.1498

# plot model fits
dta_m1a <- data.frame(dta1, fit = fitted(m1a_g), fit1 = fitted(m1_g))

ggplot(dta_m1a, aes(pSize, fit)) +
 geom_point(aes(pSize, logRS), pch = 1)+
 geom_line(col ="steelblue", lwd = 1)+
 geom_line(aes(pSize, fit1), col = "indianred", lwd = 1) +
 facet_wrap( ~ ID) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

#  individual profiles on one panel
ggplot(dta_m1a, aes(pSize, fit, group = ID)) +
 geom_point(aes(pSize, logRS), pch = 1, alpha = .6)+
 geom_line(alpha = .3, col="steelblue")+
 geom_line(aes(pSize, fit1), col = "indianred", alpha = .3) +
 labs(x = "Print size (in logMAR)", y = "Reading speed (in logWPM)")

Exercise 2

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}

# load packages for later use
pacman::p_load(bdots, tidyverse, nlme)

# load 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   Cohort

str(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)

# target looks different, select target looking time course data
dta_t <- subset(ci, LookType == "Target")
# 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)

# define the 3-parameter curve fitting function
fx <- function(x, P, S, C) {
               P / (1 + exp(S*(C - x)/P))
}

lm(Fixations~ Time, data= ci)

## 
## Call:
## lm(formula = Fixations ~ Time, data = ci)
## 
## Coefficients:
## (Intercept)         Time  
##    0.039572     0.000134

slope is about 0.0001()

# fit nls separately for all subjects
## use groupedDate because want to discuss individual differences
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.5

# covert fitting separately by subject to subjects as random effects
m1 <- nlme(m0)

# extract mean parameters values to use as initial values
prm0 <- round(unlist(lapply(m1$coef$fixed, mean)), digits=3)

# 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),
             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))))

# 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))

# 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.30290 -1.01277 -0.47429 -0.06974  1.58938 
## 
## Number of Observations: 27054
## Number of Groups: 54

# check the auto-correlation function(ACF) 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*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")

Exercise 3

The data are for 200 patients who had suffered from traumatic brain injuries resulting in comas of varying duration. After awakening from their comas, patients were administered a standard IQ test, on average, twice per patient.

Source: Wong, P.P., Monette, G., & Weiner, N.I. (2001). Mathematical models of cognitive recovery. Brain Injury, 15, 519-530.

Data: Wong{car}

pacman::p_load(tidyverse, nlme, car)
dta3 <- Wong
# check data structure
str(dta3)

## 'data.frame':    331 obs. of  7 variables:
##  $ id      : int  3358 3535 3547 3592 3728 3790 3807 3808 4253 4356 ...
##  $ days    : num  30 16 40 13 19 13 37 31 40 31 ...
##  $ duration: int  4 17 1 10 6 3 5 7 3 7 ...
##  $ sex     : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
##  $ age     : num  20.7 55.3 55.9 61.7 30.1 ...
##  $ piq     : int  87 95 95 59 67 76 74 91 115 86 ...
##  $ viq     : int  89 77 116 73 73 69 77 110 110 83 ...

# see the first 6 lines
head(dta3)

##     id days duration  sex   age piq viq
## 1 3358   30        4 Male 20.67  87  89
## 2 3535   16       17 Male 55.29  95  77
## 3 3547   40        1 Male 55.92  95 116
## 4 3592   13       10 Male 61.66  59  73
## 5 3728   19        6 Male 30.13  67  73
## 6 3790   13        3 Male 57.06  76  69

# plot individuals
## gam=generalized additive model
ggplot(dta3, aes(sqrt(days), piq, group = id)) +
 geom_point(pch = 20) +
 geom_line() +
 stat_smooth(aes(group=1), method = "gam", formula = y ~ s(x)) +
 labs(x = "Days after recovering from coma (square-root)",
      y = "Performance IQ score")

# examine the piq against coma duratin (sqrt)
ggplot(dta3, aes(sqrt(duration), piq)) +
 geom_point(pch = 20) +
 stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) +
 labs(x = "Duration of coma (square-root)", y = "Performance IQ score")

# use this to guess the initial values for parameters
summary(lm(piq ~ sqrt(duration), data = dta3))

## 
## Call:
## lm(formula = piq ~ sqrt(duration), data = dta3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -40.30 -10.89  -1.73   9.33  45.34 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)      91.304      1.275   71.61  < 2e-16
## sqrt(duration)   -1.289      0.337   -3.82  0.00016
## 
## Residual standard error: 14.8 on 329 degrees of freedom
## Multiple R-squared:  0.0425, Adjusted R-squared:  0.0396 
## F-statistic: 14.6 on 1 and 329 DF,  p-value: 0.000158

• The asympotic IQ is about 100.
• The average rate of increase, relative to coma duration, is about -1.3(by regression coefficience)
• The recoverable IQ is about 15
• The change in expected eventual PIQ gain with a 1-unit increase in length of the coma on the sqrt scale in about 0
• The half-recovery time is about 100, yielding log(2)/100=0.0069

# fit the model
m1 <- nlme(piq ~ b0 + b1*exp(-b2*days), data = dta3,
           fixed = list(b0 ~ 1 + sqrt(duration),
                        b1 ~ 1 + sqrt(duration),
                        b2 ~ 1),
           random = list(id = list(b0 ~ 1, b1 ~ 1)),
           start = list(fixed = c(100, -1.5, -10, 0, 0.007)))
# show output summary
summary(m1)

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: piq ~ b0 + b1 * exp(-b2 * days) 
##  Data: dta3 
##    AIC  BIC logLik
##   2593 2628  -1288
## 
## Random effects:
##  Formula: list(b0 ~ 1, b1 ~ 1)
##  Level: id
##  Structure: General positive-definite, Log-Cholesky parametrization
##                StdDev Corr  
## b0.(Intercept) 13.769 b0.(I)
## b1.(Intercept)  2.601 -0.996
## Residual        6.736       
## 
## Fixed effects: list(b0 ~ 1 + sqrt(duration), b1 ~ 1 + sqrt(duration), b2 ~ 1) 
##                    Value Std.Error  DF t-value p-value
## b0.(Intercept)     97.09     2.037 127   47.68  0.0000
## b0.sqrt(duration)  -1.25     0.480 127   -2.59  0.0107
## b1.(Intercept)    -11.15     3.208 127   -3.47  0.0007
## b1.sqrt(duration)  -3.25     1.077 127   -3.02  0.0031
## b2                  0.01     0.002 127    5.00  0.0000
##  Correlation: 
##                   b0.(I) b0.s() b1.(I) b1.s()
## b0.sqrt(duration) -0.724                     
## b1.(Intercept)    -0.596  0.463              
## b1.sqrt(duration)  0.463 -0.455 -0.789       
## b2                -0.309  0.013  0.092 -0.380
## 
## Standardized Within-Group Residuals:
##       Min        Q1       Med        Q3       Max 
## -3.332390 -0.365700  0.008978  0.382756  2.303091 
## 
## Number of Observations: 331
## Number of Groups: 200

# augment data with fitted values
dta_m1 <- data.frame(dta3, fit = predict(m1))

# fitted lines
ggplot(dta_m1, aes(sqrt(days), piq, group = id)) +
 geom_point(pch = 20, alpha = .5) +
 geom_line(aes(sqrt(days), fit, group = id)) +
 labs(x = "Days after recovering from coma (square-root)",
      y = "Performance IQ score")

# plot mode fits - basic graphics - detail contol
#
newdata <- expand.grid(duration = c(1, 10, 20, 50, 100, 200),
                       days = seq(0, 1000, 20))
newdata$piq <- predict(m1, newdata, level = 0)

plot(piq ~ days, type="n", data = newdata,
     xlab = "Days Post Coma", ylab = "Average Performance IQ",
     ylim = c(15, 105), xlim = c(-100, 1000), axes = FALSE, frame=TRUE)
axis(2) #right side y axis
axis(4) #left side y axis
axis(1, at = seq(0, 1000, by=100)) # X axis

# choose 1, 10, 20, 50, 100, 200的curve
for (dc in c(1, 10, 20, 50, 100, 200)){
 with(newdata, {
      lines(spline(seq(0, 1000, 20), piq[duration == dc]), lwd = 2)
      text(-25, piq[duration == dc][1], labels = dc, adj = 0.9)
 }
)}
grid()
text(-100, 95, labels="Length of\nComa", adj = 0, cex = .7)