Accroding to your method, I used mean-2sd to produce a data. However it is possible to have negative value. So I modify it as [mean-sd, mean-0.5sd, mean, mean+0.5 mean+sd], in which minimum value is mean-sd, positive.
I compare results using R codes developed by me with FDA, and the results are “same”
rawdata <- read.csv("Z:/My Project/data.csv", header=FALSE)
names(rawdata) <- c('t','m1','s1','m2','s2','m3','s3')
time <- rawdata$t
time1 <- c(rep(time[1],5),rep(time[2],5),rep(time[3],5),rep(time[4],5))
sce1 <- c(rawdata$m1, rawdata$m1 + rawdata$s1, rawdata$m1 + 0.5*rawdata$s1,
rawdata$m1 - rawdata$s1, rawdata$m1 - 0.5*rawdata$s1) ## here, I minus one sd.
sce2 <- c(rawdata$m2, rawdata$m2 + rawdata$s2, rawdata$m2 + 0.5*rawdata$s2,
rawdata$m2 - rawdata$s2, rawdata$m2 - 0.5*rawdata$s2)
sce3 <- c(rawdata$m3, rawdata$m3 + rawdata$s3, rawdata$m3 + 0.5*rawdata$s3,
rawdata$m3 - rawdata$s3, rawdata$m3 - 0.5*rawdata$s3)
sce <- array(c(sce1,sce2,sce3),dim = c(20,3))
data1 <- cbind(rep(time,5),sce)
data1 <- data1[order(data1[,1]),]
t <- data1[,1]
n <- length(t)
ni=rep(5,4) - 1 # degree of freedom
tol <- 0
tolimit <-array(1:80,dim = c(20,4))
tol_fda1 <-array(1:80,dim = c(20,4))
bartlet <- c(1:3)
##--------------------------------------------------------------------------------------##
for (j in 2:4)
{
con <- log(data1[,j])
si <- 0
for (i in 1:4)
si[i] <- var(con[(5*(i - 1) + 1):(5*i)])# varance of each group at each time point
g<-length(si)
dg<- sum(ni)
s <- sum(ni*si)/dg # total varance
C <-1 + (sum(1/ni) - 1/dg)*(1/(3*(g - 1)))
M <- dg*log(s)- sum(ni*log(si))
B <- M/C #Bartllet test
bartlet[j-1] <- pchisq(B, df = g - 1, lower.tail=FALSE)#chi square test
lin_y <- lm(con~data1[,1]) #linear model
summ <- summary(lin_y) # summary the linear model result
rse <- summ$sigma #Residual standard error
coeff <- summ$coefficients
a <- coeff[1,1]
b <- coeff[2,1]
z <- qnorm(0.99)# 99% percentile
d <- z/(1/n + (t - mean(t))^2/sum((t - mean(t))^2))^0.5
k <- qt(0.95,n-2,d,lower.tail = T)
#tolerance limit
logtole <- a + b*t + k*rse*(1/n + (t - mean(t))^2/sum((t - mean(t))^2))^0.5
tol <- exp(logtole) # withdrawal time
tolimit[,j] <- tol
# withdrawal time by FDA R codes
library(reschem)
tol_fda <- tolerance_band(lm(con ~ t), .99, .95, intervals = n ,ylog = TRUE)
names(tol_fda) <- c("v1")
tol_fda1[,j] <- t(tol_fda$v1)
}## Loading required package: stats4##
## ***********************************
## ** Welcome to package:reschem! **
## ***********************************tolerance_limit <- data.frame(time=t,scenario1 = tolimit[,2],scenario2 = tolimit[,3],
scenario3 = tolimit[,4])
t_fda <- seq(min(t), max(t), length.out = n)
tolerance_limit_fda <- data.frame(time=t_fda,scenario1 = tol_fda1[,2],scenario2 = tol_fda1[,3],
scenario3 = tol_fda1[,4])
bartlet## [1] 0.4075028 0.4075028 0.4075028tolerance_limit## time scenario1 scenario2 scenario3
## 1 24 14.045090 4.5506090 42.977974
## 2 24 14.045090 4.5506090 42.977974
## 3 24 14.045090 4.5506090 42.977974
## 4 24 14.045090 4.5506090 42.977974
## 5 24 14.045090 4.5506090 42.977974
## 6 48 5.751735 1.8635620 17.600308
## 7 48 5.751735 1.8635620 17.600308
## 8 48 5.751735 1.8635620 17.600308
## 9 48 5.751735 1.8635620 17.600308
## 10 48 5.751735 1.8635620 17.600308
## 11 72 2.567176 0.8317649 7.855557
## 12 72 2.567176 0.8317649 7.855557
## 13 72 2.567176 0.8317649 7.855557
## 14 72 2.567176 0.8317649 7.855557
## 15 72 2.567176 0.8317649 7.855557
## 16 96 1.248805 0.4046129 3.821344
## 17 96 1.248805 0.4046129 3.821344
## 18 96 1.248805 0.4046129 3.821344
## 19 96 1.248805 0.4046129 3.821344
## 20 96 1.248805 0.4046129 3.821344tolerance_limit_fda## time scenario1 scenario2 scenario3
## 1 24.00000 14.045090 4.5506090 42.977974
## 2 27.78947 12.143109 3.9343674 37.157914
## 3 31.57895 10.514423 3.4066729 32.174133
## 4 35.36842 9.119137 2.9546005 27.904561
## 5 39.15789 7.923161 2.5671042 24.244873
## 6 42.94737 6.897346 2.2347402 21.105880
## 7 46.73684 6.016766 1.9494322 18.411304
## 8 50.52632 5.260098 1.7042717 16.095899
## 9 54.31579 4.609105 1.4933499 14.103860
## 10 58.10526 4.048195 1.3116153 12.387478
## 11 61.89474 3.564053 1.1547530 10.906001
## 12 65.68421 3.145317 1.0190829 9.624672
## 13 69.47368 2.782320 0.9014716 8.513898
## 14 73.26316 2.466844 0.7992573 7.548541
## 15 77.05263 2.191926 0.7101839 6.707292
## 16 80.84211 1.951675 0.6323428 5.972126
## 17 84.63158 1.741118 0.5641221 5.327820
## 18 88.42105 1.556055 0.5041620 4.761530
## 19 92.21053 1.392948 0.4513151 4.262420
## 20 96.00000 1.248805 0.4046129 3.821344