In emkayoh/Dark: The Analysis of Dark Adaptation Data

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. \emph{The Journal of Physiology,} 1975.

A 'dark' object can have up to 16 elements and no fewer than 2.

library(Dark)
tmp <- dark
names(tmp)

load("dark.rda")
tmp<-dark
Names <- c("time", "thrs", "fit", "resid", "R2", "Bootstrap", "weight", 
"valid", "Mod", "Pn", "AIC", "data", "opt", "call", "val")
print(Names)

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

Pd <- structure(c(-1.93048394296102, -0.966243347355156, -0.990918795667905, 
0.253616641335044, 0.420301711722179, 1.32959014056482, 3.97646484343339, 
0.763150661132819, 0.905683974690625, -0.389894947473737, -0.4064032016008, 
-0.182491245622811, 2.52208104052026, 6.96883941970985, 3.12198099049525, 
0.130565282641321, 0.207403701850925, 0.0885442721360681, 32.5607121101566, 
32.9919673191487, 29.4211741886944), .Dim = c(3L, 7L), .Dimnames = list(
    NULL, c("CT", "CC", "Tau", "S2", "Alph", "S3", "Beta")))
print.table(Pd,3)

The 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

MSCd <- structure(list(AIC = c(0, 0, -347.688152541027, 0, -353.337417397606, 
-556.296993756504, -561.687551083725), param = structure(c(-7.03677859813905, 
-4.84533369033486, -5.52767312215107, -1.97221068357342, 6.23084047079413, 
4.01859727139177, 1.43155445018767, 1.39834034016783, 28.1212523680337, 
17.2254594347471, 0.557811726205019, 1.62283784201578, 0, -0.0365810342284519, 
9.61920222958987, -0.227896945728703, 0, 9.40371211032976, 3.50297041829524, 
8.73599874809917, 0, 0, 0.0972021474126348, 0.180809685149208, 
0, 0, 0, 19.4250494832971), .Dim = c(4L, 7L))), .Names = c("AIC", 
"param"))
print(MSCd, 3)

Then 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;

source("../R/H.R")
source("../R/P7c.R")
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)
X <- tmp$time
lines(X, tmp$fit)
lines(X, P7c(tmp$Bootstrap[,1], X), col = 2)
lines(X, P7c(tmp$Bootstrap[,3], X), col = 2)