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)
The final dark
object with 15 elements:
print(tmp,2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.