The following describes the work flow for the analysis of dark adaptation data. The data were extracted from Figure 1 in Rushton’s paradox: rod dark adaptation after flash photolysis, E.N.Pugh Jr. 1975.
A ‘dark’ object can have up to 16 elements and no fewer than 2.
library(Dark)
tmp <- dark
names(tmp)## [1] "time" "thrs" "fit" "resid" "R2"
## [6] "Bootstrap" "weight" "valid" "Mod" "Pn"
## [11] "AIC" "data" "opt" "call" "val"The raw data are shown below
par(las=1, bty='n',mfrow=c(1,1))
XL <- expression(bold(Time~(min)))
YL <- expression(bold(Threshold~(LU)))
plot(tmp$time, tmp$thrs, xlab = XL, ylab=YL)
The function Start generates an array of possible parameter values, P, calculated from the data, which are then passed to ModelSelect.
set.seed(1234)
P <-Start(tmp, 2000)
head(P,3)The first three rows are shown below
## CT CC Tau S2 Alph S3 Beta
## [1,] -1.9305 0.2536 3.9765 -0.3899 2.5221 0.1306 32.5607
## [2,] -0.9662 0.4203 0.7632 -0.4064 6.9688 0.2074 32.9920
## [3,] -0.9909 1.3296 0.9057 -0.1825 3.1220 0.0885 29.4212The data frame P is processed by ModelSelect, this finds the row of parameters that give the lowest sum of squared errors and calculates an AIC score. This is not calculated again.
MSC<-ModelSelect(tmp,P)
MSC## $AIC
## [1] 0 0 -348 0 -353 -556 -562
##
## $param
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -7.04 6.23 28.121 0.0000 0.00 0.0000 0.0
## [2,] -4.85 4.02 17.225 -0.0366 9.40 0.0000 0.0
## [3,] -5.53 1.43 0.558 9.6192 3.50 0.0972 0.0
## [4,] -1.97 1.40 1.623 -0.2279 8.74 0.1808 19.4Then the function BestFit optimises the model using the lowest AIC score to select the model and the best initial estimate of the parameters for that model.
tmp<-BestFit(tmp,MSC)A further function MultiStart ensures that the estimated parameter values are optimal by repeatedly optimising the model with starting locations in the vicinity of the initial estimate. Finally BootDark uses bootstrap methods to calculate confidence intervals for the parameter estimates.
tmp<-MultiStart(tmp,repeats = 25)
tmp<-BootDark(tmp,R = 500)The data and model are shown in the figure below; 
The final dark object with 15 elements:
print(tmp,2)## $time
## [1] -0.0494 0.0014 0.0521 0.1503 0.2518 0.2518 0.4022 0.4040
## [9] 0.6788 0.7541 0.8574 0.9423 1.2075 1.2109 1.4087 1.4202
## [17] 1.5907 1.6476 1.8917 1.9293 2.1700 2.3712 2.4448 2.5970
## [25] 3.2045 3.5825 3.6315 4.0226 4.3369 4.4499 5.0504 5.5658
## [33] 5.6409 6.0925 6.4575 6.7097 7.1385 7.5279 7.9567 8.1581
## [41] 8.2332 8.6490 9.0288 9.0550 9.3430 9.3577 9.8603 10.0694
## [49] 10.1447 10.1560 10.3329 10.4229 10.8483 11.0006 11.0006 11.1398
## [57] 11.3017 11.4558 11.5931 11.7436 12.0412 12.1724 12.2836 12.4508
## [65] 12.6389 12.8401 12.9775 13.1544 13.5062 13.5683 13.8187 13.8466
## [73] 14.4584 14.5092 14.6729 14.7596 14.8137 15.0607 15.3880 15.4370
## [81] 15.8789 16.0408 16.2700 16.7234 17.3011 17.3747 17.8655 17.8788
## [89] 18.1833 18.3715 18.6708 18.6821 19.1750 19.5118 19.6248 20.1551
## [97] 20.3958 20.7100 20.8735 21.7427 21.7427 22.5593 22.5967 22.8979
## [105] 24.5538 25.2689 26.2867 26.5879 26.7891 27.1033 27.8822 27.9199
## [113] 28.3618 29.3911 29.4664 29.9950 30.8361 31.3875 31.6263 32.0159
## [121] 32.4691 32.9322 34.0630 34.3019 34.5785 35.4178 36.5502 36.6517
## [129] 37.5339 37.9743 38.9938 39.7500 41.4565 41.8986 42.9541 44.4235
## [137] 46.4952 49.4832
##
## $thrs
## [1] -0.57 -0.25 -0.74 -0.85 -0.98 -1.03 -0.33 -1.15 -0.82 -0.92 -1.06
## [12] -1.27 -1.37 -1.55 -1.49 -1.22 -1.59 -1.34 -1.52 -1.64 -1.67 -1.72
## [23] -1.63 -1.75 -1.74 -1.75 -1.84 -1.77 -1.79 -1.86 -1.92 -1.88 -1.83
## [34] -1.96 -1.92 -1.99 -1.95 -2.00 -1.97 -2.03 -1.95 -2.01 -2.06 -2.14
## [45] -2.25 -2.15 -2.08 -2.33 -2.25 -2.18 -2.41 -2.51 -2.34 -2.14 -2.66
## [56] -2.52 -2.59 -2.73 -2.65 -2.60 -2.38 -2.73 -2.69 -2.72 -2.78 -2.92
## [67] -2.98 -3.03 -3.02 -3.20 -3.10 -3.27 -3.34 -3.28 -3.08 -3.47 -3.21
## [78] -3.56 -3.64 -3.41 -3.75 -3.60 -3.83 -3.94 -3.95 -4.06 -4.17 -4.06
## [89] -4.02 -4.07 -4.19 -4.26 -4.11 -4.29 -4.36 -4.25 -4.38 -4.29 -4.42
## [100] -4.48 -4.55 -4.59 -4.41 -4.56 -4.67 -4.67 -4.74 -4.80 -4.72 -4.87
## [111] -4.95 -4.85 -4.89 -4.90 -4.99 -4.97 -5.04 -4.96 -5.11 -5.09 -5.08
## [122] -5.14 -5.10 -5.22 -5.26 -5.15 -5.19 -5.29 -5.33 -5.33 -5.36 -5.40
## [133] -5.35 -5.40 -5.52 -5.46 -5.55 -5.59
##
## $fit
## [1] -0.53 -0.58 -0.62 -0.70 -0.78 -0.78 -0.88 -0.88 -1.05 -1.10 -1.15
## [12] -1.19 -1.31 -1.31 -1.39 -1.39 -1.45 -1.47 -1.54 -1.55 -1.61 -1.65
## [23] -1.66 -1.69 -1.78 -1.82 -1.82 -1.85 -1.87 -1.88 -1.91 -1.92 -1.92
## [34] -1.94 -1.94 -1.95 -1.95 -1.95 -1.96 -1.96 -1.96 -1.96 -2.04 -2.04
## [45] -2.11 -2.11 -2.23 -2.27 -2.29 -2.29 -2.34 -2.36 -2.45 -2.49 -2.49
## [56] -2.52 -2.56 -2.59 -2.62 -2.66 -2.73 -2.76 -2.78 -2.82 -2.86 -2.91
## [67] -2.94 -2.98 -3.06 -3.07 -3.13 -3.14 -3.28 -3.29 -3.33 -3.35 -3.36
## [78] -3.41 -3.49 -3.50 -3.60 -3.64 -3.69 -3.79 -3.92 -3.94 -4.05 -4.06
## [89] -4.13 -4.17 -4.24 -4.24 -4.35 -4.41 -4.42 -4.44 -4.45 -4.47 -4.48
## [100] -4.52 -4.52 -4.56 -4.56 -4.57 -4.65 -4.68 -4.73 -4.75 -4.76 -4.77
## [111] -4.81 -4.81 -4.83 -4.88 -4.88 -4.91 -4.95 -4.97 -4.98 -5.00 -5.02
## [122] -5.04 -5.10 -5.11 -5.12 -5.16 -5.22 -5.22 -5.26 -5.28 -5.33 -5.37
## [133] -5.45 -5.47 -5.52 -5.59 -5.68 -5.82
##
## $resid
## [1] -0.04057 0.32852 -0.12241 -0.14992 -0.20206 -0.25384 0.55760
## [8] -0.26303 0.23789 0.17606 0.08952 -0.07589 -0.05625 -0.23738
## [15] -0.10452 0.17509 -0.13863 0.13047 0.01824 -0.09030 -0.06555
## [22] -0.06909 0.03569 -0.06254 0.03775 0.06644 -0.01961 0.08213
## [29] 0.08383 0.02182 -0.01814 0.04462 0.09890 -0.02458 0.02277
## [36] -0.04558 -0.00152 -0.04950 -0.00939 -0.07276 0.00544 -0.04887
## [43] -0.02306 -0.09480 -0.14492 -0.04355 0.14787 -0.05592 0.04561
## [50] 0.11877 -0.07163 -0.15046 0.11309 0.34460 -0.17331 -0.00472
## [57] -0.03247 -0.13649 -0.03069 0.06100 0.34684 0.02158 0.09151
## [64] 0.09810 0.07984 -0.01080 -0.03870 -0.05378 0.03726 -0.12941
## [71] 0.03119 -0.12933 -0.06577 0.01225 0.24825 -0.12068 0.14472
## [78] -0.14236 -0.14750 0.08516 -0.15200 0.03274 -0.14264 -0.15206
## [85] -0.02241 -0.12363 -0.11706 -0.00689 0.11040 0.09954 0.04974
## [92] -0.02161 0.24014 0.12145 0.05849 0.19113 0.07073 0.17995
## [99] 0.05282 0.04200 -0.03578 -0.03387 0.15066 0.01532 -0.01933
## [106] 0.01437 -0.01158 -0.04940 0.03231 -0.09544 -0.13820 -0.04202
## [113] -0.05997 -0.02563 -0.10566 -0.06220 -0.09647 0.00728 -0.13073
## [120] -0.08466 -0.05224 -0.09880 0.00074 -0.10790 -0.13918 0.00946
## [127] 0.02047 -0.07469 -0.06827 -0.04751 -0.03269 -0.03416 0.09975
## [134] 0.06880 0.00200 0.13024 0.13013 0.23554
##
## $R2
## [1] 0.99
##
## $Bootstrap
## 2.5% 50% 97.5%
## CT -2.04 -1.97 -1.89
## CC 1.29 1.39 1.50
## Tau 1.32 1.61 1.99
## S2 -0.24 -0.23 -0.22
## Alpha 8.32 8.71 9.06
## *S3* 0.17 0.18 0.19
## Beta 18.95 19.45 19.98
##
## $weight
## CT CC Tau S2 Alpha *S3* Beta
## 6.80 4.67 1.49 50.00 1.35 47.62 0.97
##
## $valid
## [1] 1 1 1 1 1 1 1
##
## $Mod
## [1] "P7c"
##
## $Pn
## [1] 7
##
## $AIC
## [1] 0 0 -348 0 -356 0 -562
##
## $opt
## [1] -1.97 1.39 1.61 -0.23 8.70 0.18 19.43
##
## $val
## [1] 2.1
##
## $call
## BootDark(obj = tmp, R = 500, graph = T)
##
## $data
## [1] "Pugh_1975"
##
## attr(,"class")
## [1] "dark"