#_______________________________________________________________________________________________________________________ library(tidyverse) library(gslnls) library(ggplot2) library(reshape2) library(drc) all_data <- data.frame(Conc= c( 80000.000, 10666.667, 1422.2222, 189.62963, 25.28395, 3.37119, 0.44949, 0.05993), ST_1 = c(414356, 402828, 405696, 510624 ,844864, 995120, 984952, 1076052), ST_2 = c(407280, 398916, 410140, 530372, 851828, 997708, 1001036, 1069060), ST_3 = c(371408, 424764, 397564, 498328, 808188, 910116, 995972, 974076), REF_70_1 = c( 402432, 413852, 422292, 579012, 896164, 1029552, 1048636 ,994528), REF_70_2 = c(401548, 406580, 429472, 560256, 895536, 976660, 945068, 1025772), REF_70_3 = c(390632, 386176, 418596, 536832, 865956, 888640, 985340, 1053968)) # PLA GHZ250925nM2.docx all_data <- data.frame(Conc= c( 80000.000, 10666.667, 1422.2222, 189.62963, 25.28395, 3.37119, 0.44949, 0.05993), ST_1 = c(138743, 138450, 156121, 262387, 561683, 652307, 701813, 730525), ST_2 = c(130197, 137202, 145935, 271407, 580252, 645236, 700294, 727128), ST_3 = c(127476, 130984, 150794, 252283, 559926, 648084, 701467, 683004), REF_70_1 = c(141259, 136659, 157568, 255747, 557256, 693976, 683018, 687692), REF_70_2 = c(132591, 135163, 151050, 258229, 535896, 615267, 672638, 697947), REF_70_3 = c(129988, 124818, 142684, 237172, 514545, 622616, 668689, 652239)) # PLA GHZ339_20251022_1MAK all_data <- data.frame(Conc= c( 80000.000, 10666.667, 1422.2222, 189.62963, 25.28395, 3.37119, 0.44949, 0.05993), GHZ339.01_REF_1 = c( 141179, 137438, 156331, 220605, 457358, 561383, 583467, 599980), GHZ339.01_REF_2= c( 142449, 145681, 151674, 242837, 458059, 575359, 596014, 612299), GHZ339.01_REF_3= c( 142608, 137862, 149732, 238484, 467783, 564427, 620924, 634909), GHZ339DS_1 = c( 137179, 156595, 155272, 227545, 460416, 569905, 610124, 592936), GHZ339DS_2= c( 139982, 144452, 161241, 243495, 458981, 564540, 613269, 605241), GHZ339DS_3 = c( 153642, 145111, 158542, 229435, 458629, 556951, 612850, 614633)) # for information: 4-pl-model evaluated with drc-package all_l <- melt(all_data, id.vars= "Conc", variable.name="replname", value.name="readout") all_l$logreadout <- log(all_l$readout) # add sample ID all_l$sample <- c(rep(0,nrow(all_l)/2), rep(1,nrow(all_l)/2)) all_l$isRef <- c(rep(1,nrow(all_l)/2), rep(0,nrow(all_l)/2)) #all_l$sample <- c(rep("A",24), rep("B",24)) all_l$log_dose <- log(all_l$Conc) # ' Achtung: logreadout' pot <- drm(readout ~ Conc, sample, data = all_l, fct = LL.4(names = c("b", "d", "a", "c")), pmodels = data.frame(1, 1, 1, sample)) potcomp <- drc::compParm(pot, "c", display = T) confint(pot, display = T) summary(pot) supot <- summary(pot)$coefficients var_log_s <- (supot[4,2]/supot[4,1])^2 # (se_b0 / b0)^2 # Approx variance of log(b) var_log_t <- (supot[5,2]/supot[5,1])^2 #(se_b1 / b1)^2 # Approx variance of log(b) #CoVarlog <- (CoSE/mean(c(b0,b1)))^2 # Simplified CI for log(b1/b2) (ignoring covariance for simplicity) se_log_ratio <- sqrt(var_log_s + var_log_t) #-2*CoVarlog) # log-pot supot[4,1]/supot[5,1] loguCI <- log(supot[4,1]) -log(supot[5,1])+qt(0.975,32)*se_log_ratio loglCI <- log(supot[4,1]) -log(supot[5,1])-qt(0.975,32)*se_log_ratio exp(c(loglCI,loguCI)) ################################################################################ # calculation from GHZ PLA printout # 1.00903 (0.90187 1.12894) ################################################################################ EDcomp(pot, percVec=c(50,50), interval="delta", display=T) EDcomp(pot, percVec=c(50,50), interval="fieller", display=T) #Estimate Lower Upper #0/1:50/50 1.13150 0.99098 1.27201 #SE_ <- -(0.804-0.902186)/qt(0.975,44) exp(potcomp[1]+qt(0.975, 43)*potcomp[2]) # = delta methode #[1] 1.272015 # Delta method for ratios deltaS <- 1/supot[5,1]^2*(vcov(pot)[4,4]-2*potcomp[1,1]*vcov(pot)[4,5]+potcomp[1,1]^2*vcov(pot)[5,5]) potcomp[1,1]+qt(0.975,48-5)*sqrt(deltaS) potcomp[1,1]-qt(0.975,48-5)*sqrt(deltaS) # Fieller method EDcomp(pot, percVec=c(50,50), interval="fieller", display=T) #Estimate Lower Upper #0/1:50/50 1.13150 0.99848 1.28281 potcomp[1]+qnorm(0.975)*potcomp[2] modelFit(pot) # Lack-of-fit test # # ModelDf RSS Df F value p value # ANOVA 32 0.050984 # DRC model 43 0.197391 11 8.3538 0.0000 startlist <- list(a = min(all_l$readout), b = -1.1, cs=mean(all_l$log_dose), d = max(all_l$readout), r = 0) mr <- gsl_nls(fn = readout ~ a + (d - a)/(1 + exp(b*(log_dose - (cs-r*sample)))), data = all_l, start = startlist, algorithm = "lmaccel", control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)) s_mr <- summary(mr)$coefficients (s_mr) #sum(pot$predres[,2]^2) # [1] 0.197391 RelPot_mr <- exp(s_mr[5,1]) # r exponentiated print(paste("Potency restricted",RelPot_mr)) (POTr_CI <- exp(confint(mr, "r", method = "asymptotic"))) # mit profile Likelihood: (POTr_CI <- exp(confint(mr, "r", method = "profile"))) potcomp potcomp[1]+qt(0.975,43)*potcomp[2] potcomp[1]-qt(0.975,43)*potcomp[2] print(POTr_CI) # 2.5 % 97.5 % # r 0.9994816 1.280957 # exact wie Stegmann s_mr[5,1]+qt(0.975,43)*s_mr[5,2] #[1] 0.2476073 exp(s_mr[5,1]+qt(0.975,43)*s_mr[5,2]) #[1] 1.280957 # exact wie Stegmann exp(s_mr[5,1]-qt(0.975,43)*s_mr[5,2]) #[1] 0.9994816 # exact wie Stegmann plot(all_l$log_dose, all_l$readout) # weighting w <- aggregate(all_l$readout, by=list(all_l$log_dose), FUN=mean) splt <- split(all_l, f=c(all_l$log_dose)) Res <- sapply(1:8, function(x) splt[[x]]['readout']-w[,2][x]) SD <- aggregate(all_l$readout, by=list(all_l$log_dose), FUN=sd) ResW <- sapply(1:8, function(x) Res[[x]] / SD[,2][x]) # PLA does not do weighting # mr <- gsl_nls(fn = readout ~ a + (d - a)/(1 + exp(b*(log_dose - (cs-r*sample)))), # data = all_l, # start = startlist, weights=c(rep(SD$x/100,6)), # control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)) # s_mr <- summary(mr)$coefficients # (s_mr) # exp(s_mr[5,1]) library(calculus) jac <- calculus::jacobian(f="a+(d-a)/(1+exp(b*(c-r*sample-x)))" , var=c("d","b","c","a","r")) resjac <- matrix(NA,nrow=48, ncol=5) for (n in 1:nrow(all_l)) { Vec <- c(a=s_mr[4,1], d=s_mr[1,1], c=s_mr[3,1], b=s_mr[2,1], r=s_mr[5,1], sample=all_l$sample[n], x=all_l$log_dose[n]) # resjac[n,] <- evaluate(jac, var=Vec) } identical(resjac, mr$m$gradient()) # eigentlich gleich, aber col3 und 4 vertauscht; resjac- mr$m$gradient() #_______________________________________________________________________________ colnames(resjac) <- c("dd", "db", "ddc", "da", "dr") V_V <- solve(t(resjac)%*%resjac) RMSE <- sqrt(sum(mr$m$resid()^2)/43) VCOV <- V_V*RMSE^2 SEs <- sqrt(diag(V_V)*RMSE^2) # stimmen mit summary(mr)$coefficients ueberein!! SEs2 <- sqrt(diag(vcov(mr))) #_______________________________________________________________________________ FitAnova <- anova(lm(readout ~ factor(Conc)*sample, all_l)) # Pure error pureSS <- FitAnova[4,2] pureSS_df <- FitAnova[4,1] meanPureErr <- FitAnova[4,3] SEsPure <- sqrt(diag(V_V)*meanPureErr) # PLA software: 1.0252 - 1.4242 exp(s_mr[5,1]+qt(0.975,43)*SEsPure[5]) # dr # 1.422001 exp(s_mr[5,1]-qt(0.975,43)*SEsPure[5]) # dr # 1.027014 # Was PLA macht: 16 df wegen 8+8 Mittelwerte exp(s_mr[5,1]+qt(0.975,48-16)*SEsPure[5]) exp(s_mr[5,1]-qt(0.975,48-16)*SEsPure[5]) # rel CIs exp(s_mr[5,1]+qt(0.975,48-16)*SEsPure[5])/exp(s_mr[5,1]) exp(s_mr[5,1]-qt(0.975,48-16)*SEsPure[5])/exp(s_mr[5,1]) CORro <- cor(all_l$log_dose[1:8], all_l$readout[1:8]) if (CORro<0) B <- -1 else B <- 1 startlstmu <- list(as = min(all_l$readout), bs = B, cs = mean(all_l$log_dose), ds = max(all_l$readout), at = min(all_l$readout), bt = B, ct = mean(all_l$log_dose), dt = max(all_l$readout)) mu <- tryCatch({ gsl_nls(readout ~ as*isRef + at*sample + (ds*isRef + dt*sample - as*isRef - at*sample)/ (1 + isRef*exp(bs*(cs -log_dose)) + sample*exp(bt*(ct -log_dose))), data = all_l, start = startlstmu, control = gsl_nls_control(xtol = 1e-10, ftol = 1e-10, gtol = 1e-10)) }, error = function(msg){ return(0) } ) if (length(mu)>1) { smu_coef <- tryCatch({ summary(mu)$coefficients }, error = function(msg){ return(rep(NA, 8)) } ) } su <- summary(mu) # Folgende CIs sind symmetrisch, also nicht im log berechnet gslnls::confintd(mu, quote(dt/ds)) gslnls::confintd(mu, quote(bt/bs), level = 0.95) # fit lwr upr # bt/bs 1.009691 0.8629895 1.156392 potU <- drm(readout ~ Conc, sample, data = all_l, fct = LL.4(names = c("b", "d", "a", "c")), pmodels = data.frame(sample, sample, sample, sample)) b0 <- coef(potU)[1] b1 <- coef(potU)[2] vcov(potU) se_b0 <- sqrt(vcov(potU)[1, 1]) se_b1 <- sqrt(vcov(potU)[2, 2]) CoSE <- sqrt(abs(vcov(potU)[1, 2])) # Calculate log(b) and its variance log_b0 <- log(b0) log_b1 <- log(b1) var_log_b0 <- (se_b0 / b0)^2 # Approx variance of log(b) var_log_b1 <- (se_b1 / b1)^2 # Approx variance of log(b) CoVarlog <- (CoSE/mean(c(b0,b1)))^2 # Simplified CI for log(b1/b2) (ignoring covariance for simplicity) se_log_ratio <- sqrt(var_log_b0 + var_log_b1 -2*CoVarlog) z_score <- 1.96 # For 95% CI lower_log_ratio <- log_b1 - log_b0 - z_score * se_log_ratio upper_log_ratio <- log_b1 - log_b0 + z_score * se_log_ratio # CI for the ratio itself ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) print(ci_ratio) # b:1 b:1 # 0.8708895 1.1705629 # gsl_nls package b0 <- coef(mu)[2] b1 <- coef(mu)[6] vcov(mu) se_b0 <- sqrt(vcov(mu)[2, 2]) se_b1 <- sqrt(vcov(mu)[6, 6]) #CoSE <- sqrt(abs(vcov(mu)[6, 2])) # Calculate log(b) and its variance log_b0 <- log(abs(b0)) log_b1 <- log(abs(b1)) var_log_b0 <- (se_b0 / b0)^2 # Approx variance of log(b) var_log_b1 <- (se_b1 / b1)^2 # Approx variance of log(b) #CoVarlog <- (CoSE/mean(c(b0,b1)))^2 # Simplified CI for log(b1/b2) (ignoring covariance for simplicity) se_log_ratio <- sqrt(var_log_b0 + var_log_b1) #-2*CoVarlog) z_score <- 1.96 # For 95% CI lower_log_ratio <- log_b1 - log_b0 - z_score * se_log_ratio upper_log_ratio <- log_b1 - log_b0 + z_score * se_log_ratio # CI for the ratio itself ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) print(ci_ratio) # bt bt # 0.876991 1.162469 # Pure error FitAnova <- anova(lm(readout ~ factor(Conc)*sample, all_l)) # Pure error pureSS <- FitAnova[4,2] pureSS_df <- FitAnova[4,1] meanPureErr <- FitAnova[4,3] V_V <- vcov(mu)/su$sigma^2 SEsPure <- sqrt(diag(V_V)*meanPureErr) VCOVpure <- V_V*meanPureErr se_b0 <- sqrt(VCOVpure[2, 2]) se_b1 <- sqrt(VCOVpure[6, 6]) #CoSE <- sqrt(abs(vcov(mu)[6, 2])) # Calculate log(b) and its variance log_b0 <- log(abs(b0)) log_b1 <- log(abs(b1)) var_log_b0 <- (se_b0 / b0)^2 # Approx variance of log(b) var_log_b1 <- (se_b1 / b1)^2 # Approx variance of log(b) #CoVarlog <- (CoSE/mean(c(b0,b1)))^2 # Simplified CI for log(b1/b2) (ignoring covariance for simplicity) se_log_ratio <- sqrt(var_log_b0 + var_log_b1) #-2*CoVarlog) z_score <- 1.96 # For 95% CI lower_log_ratio <- log_b1 - log_b0 - qt(0.975,32)* se_log_ratio upper_log_ratio <- log_b1 - log_b0 + qt(0.975,32) * se_log_ratio # ------------------------------------------------------------------------------ ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio)) print(ci_ratio) # BINGOOOO # 1.009691 # bt bt # 0.8937963 1.1406122 # im Vergleich: PLA printout # 0.89354 - 1.14077, aber point est. 1.00977