R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require("DoseFinding") Loading required package: DoseFinding > data(IBScovars) > lmfit <- lm(resp~factor(dose)+gender, data=IBScovars) > cf <- coef(lmfit)[-c(6)] > vcv <- vcov(lmfit)[-c(6), -c(6)] > lmfit2 <- lm(resp~as.factor(dose)-1+gender, data=IBScovars) > cf2 <- coef(lmfit2)[-c(6)] > vcv2 <- vcov(lmfit2)[-c(6), -c(6)] > dose <- c(0:4) > > ## test fitting all available models > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="linear", placAdj=TRUE,type="general") Dose Response Model Model: linear Fit-type: general Coefficients dose-response model delta 0.07484 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 3.472 > fitMod(dose, cf2, S=vcv2, model="linear", placAdj=FALSE,type="general") Dose Response Model Model: linear Fit-type: general Coefficients dose-response model e0 delta 0.33400 0.07484 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 3.472 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="quadratic", placAdj=TRUE,type="general") Dose Response Model Model: quadratic Fit-type: general Coefficients dose-response model b1 b2 0.22800 -0.03811 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.9281 > fitMod(dose, cf2, S=vcv2, model="quadratic", placAdj=FALSE,type="general") Dose Response Model Model: quadratic Fit-type: general Coefficients dose-response model e0 b1 b2 0.24950 0.22800 -0.03811 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.9281 > fitMod(dose, cf2, S=vcv2, model="linlog", placAdj=FALSE, + addArgs=list(off=0.01*max(dose)),type="general") Dose Response Model Model: linlog Fit-type: general Coefficients dose-response model e0 delta 0.47790 0.07728 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.1605 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="emax", placAdj=TRUE, + bnds=defBnds(max(dose))$emax,type="general") Dose Response Model Model: emax Fit-type: general Coefficients dose-response model eMax ed50 0.3771 0.3630 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.09815 > fitMod(dose, cf2, S=vcv2, model="emax", placAdj=FALSE, + bnds=defBnds(max(dose))$emax,type="general") Dose Response Model Model: emax Fit-type: general Coefficients dose-response model e0 eMax ed50 0.2221 0.3771 0.3630 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.09815 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="sigEmax", placAdj=TRUE, + bnds=defBnds(max(dose))$sigEmax,type="general") Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model eMax ed50 h 0.4714 0.4880 0.5000 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.07825 > fitMod(dose, cf2, S=vcv2, model="sigEmax", placAdj=FALSE, + bnds=defBnds(max(dose))$sigEmax,type="general") Dose Response Model Model: sigEmax Fit-type: general Coefficients dose-response model e0 eMax ed50 h 0.2223 0.4714 0.4880 0.5000 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.07825 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="exponential", + placAdj=TRUE, bnds=defBnds(max(dose))$exponential,type="general") Dose Response Model Model: exponential Fit-type: general Coefficients dose-response model e1 delta 0.4391 8.0000 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 4.102 > fitMod(dose, cf2, S=vcv2, model="exponential", placAdj=FALSE, + bnds=defBnds(max(dose))$exponential,type="general") Dose Response Model Model: exponential Fit-type: general Coefficients dose-response model e0 e1 delta 0.3510 0.4391 8.0000 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 4.102 > fitMod(dose, cf2, S=vcv2, model="logistic", placAdj=FALSE, + bnds=defBnds(max(dose))$logistic,type="general") Dose Response Model Model: logistic Fit-type: general Coefficients dose-response model e0 eMax ed50 delta -0.1095 0.6665 0.0040 0.4217 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.1951 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="betaMod", placAdj=TRUE, + bnds=defBnds(max(dose))$betaMod, addArgs=list(scal=1.2*4),type="general") Dose Response Model Model: betaMod Fit-type: general Coefficients dose-response model eMax delta1 delta2 0.3468 0.2216 0.0500 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.0783 > fitMod(dose, cf2, S=vcv2, model="betaMod", placAdj=FALSE, + bnds=defBnds(max(dose))$betaMod, addArgs=list(scal=1.2*4),type="general") Dose Response Model Model: betaMod Fit-type: general Coefficients dose-response model e0 eMax delta1 delta2 0.2222 0.3468 0.2216 0.0500 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.0783 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="linInt", placAdj=TRUE, type="general") Dose Response Model Model: linInt Fit-type: general Coefficients dose-response model d1 d2 d3 d4 0.2846 0.2965 0.3502 0.3480 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0 > fitMod(dose, cf2, S=vcv2, model="linInt", placAdj=FALSE, type="general") Dose Response Model Model: linInt Fit-type: general Coefficients dose-response model d0 d1 d2 d3 d4 0.2220 0.5065 0.5185 0.5722 0.5700 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 1.484e-30 > ## test using starting value (instead of grid search) > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="emax", + placAdj=TRUE, bnds=defBnds(max(dose))$emax, start = 0.5,type="general") Dose Response Model Model: emax Fit-type: general Coefficients dose-response model eMax ed50 0.3771 0.3630 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.09815 > fitMod(dose, cf2, S=vcv2, model="emax", placAdj=FALSE, + bnds=defBnds(max(dose))$emax, start = 0.2,type="general") Dose Response Model Model: emax Fit-type: general Coefficients dose-response model e0 eMax ed50 0.2221 0.3771 0.3630 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.09815 > fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="betaMod", + placAdj=TRUE, bnds=defBnds(max(dose))$betaMod, + addArgs=list(scal=1.2*4),type="general") Dose Response Model Model: betaMod Fit-type: general Coefficients dose-response model eMax delta1 delta2 0.3468 0.2216 0.0500 Fitted to: 1 2 3 4 0.2846 0.2965 0.3502 0.3480 Generalized residual sum of squares: 0.0783 > fitMod(dose, cf2, S=vcv2, model="betaMod", placAdj=FALSE, + bnds=defBnds(max(dose))$betaMod, start = c(1, 1), + addArgs=list(scal=1.2*4),type="general") Dose Response Model Model: betaMod Fit-type: general Coefficients dose-response model e0 eMax delta1 delta2 0.2222 0.3468 0.2216 0.0500 Fitted to: 0 1 2 3 4 0.2220 0.5065 0.5185 0.5722 0.5700 Generalized residual sum of squares: 0.0783 > > ## test predict, vcov, coef, intervals, plot, summary > ggI <- fitMod(dose, cf2, S=vcv2, model="betaMod", + placAdj=FALSE, bnds=defBnds(max(dose))$betaMod, + addArgs=list(scal=1.2*4),type="general") > ggNI <- fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="betaMod", + placAdj=TRUE, bnds=defBnds(max(dose))$betaMod, + addArgs=list(scal=1.2*4),type="general") > predict(ggI, se.fit=TRUE, predType = "e") $fit [1] 0.0000000 0.2756421 0.3165257 0.3387105 0.3466573 $se.fit [1] 0.0000000 0.1235753 0.1106156 0.1135736 0.1263900 > predict(ggNI, se.fit=TRUE, predType = "e") $fit [1] 0.2756421 0.3165257 0.3387105 0.3466573 $se.fit [1] 0.1235753 0.1106156 0.1135736 0.1263900 > vcov(ggI) e0 eMax delta1 delta2 e0 0.0118410628 -0.008186899 0.003722514 -0.0008419768 eMax -0.0081868987 0.014334651 -0.016984400 -0.0239850502 delta1 0.0037225143 -0.016984400 0.475384247 0.3641344062 delta2 -0.0008419768 -0.023985050 0.364134406 0.3356088151 > vcov(ggNI) eMax delta1 delta2 eMax 0.01433465 -0.0169844 -0.02398505 delta1 -0.01698440 0.4753842 0.36413441 delta2 -0.02398505 0.3641344 0.33560882 > plot(ggI, CI=T, plotData = "meansCI") > plot(ggNI, CI=T, plotData = "meansCI") > > ggI <- fitMod(dose, cf2, S=vcv2, model="linInt", + placAdj=FALSE,type="general") > ggNI <- fitMod(dose[-1], cf[-1], S=vcv[-1,-1], model="linInt", + placAdj=TRUE,type="general") > predict(ggI, se.fit=TRUE, predType = "full-model") Message: Setting predType to "ls-means" for "type = general" $fit [1] 0.2219856 0.5065390 0.5185323 0.5721580 0.5699850 $se.fit [1] 0.1088272 0.1048057 0.1044027 0.1046862 0.1088405 > predict(ggI, se.fit=TRUE, predType = "effect-curve") $fit [1] 0.0000000 0.2845534 0.2965467 0.3501724 0.3479994 $se.fit [1] 0.0000000 0.1252903 0.1265494 0.1279290 0.1273279 > predict(ggNI, se.fit=TRUE, predType = "full-model") Message: Setting predType to "effect-curve" for placebo-adjusted data $fit [1] 0.2845534 0.2965467 0.3501724 0.3479994 $se.fit [1] 0.1252903 0.1265494 0.1279290 0.1273279 > vcov(ggI) d0 d1 d2 d3 d4 d0 0.011843358 0.003564965 0.003364270 0.003218371 0.003738601 d1 0.003564965 0.010984234 0.003307337 0.003163906 0.003675332 d2 0.003364270 0.003307337 0.010899929 0.002985790 0.003468425 d3 0.003218371 0.003163906 0.002985790 0.010959203 0.003318008 d4 0.003738601 0.003675332 0.003468425 0.003318008 0.011846244 > vcov(ggNI) d1 d2 d3 d4 d1 0.015697662 0.008221459 0.008223928 0.008215124 d2 0.008221459 0.016014746 0.008246506 0.008208911 d3 0.008223928 0.008246506 0.016365819 0.008204394 d4 0.008215124 0.008208911 0.008204394 0.016212400 > plot(ggI, CI=T, plotData = "meansCI") > plot(ggNI, CI=T, plotData = "meansCI") > > ## even more tests for the linInt model > data(IBScovars) > ## without covariates > fit <- fitMod(dose, resp, data=IBScovars, model="linInt") > plot(fit, CI=TRUE, plotData="meansCI") > fit <- fitMod(dose, resp, data=IBScovars, model="linInt", + addCovars=~gender) > plot(fit, CI=TRUE, plotData="meansCI") > vcov(fit) d0 d1 d2 d3 d4 d0 0.011843358 0.003564965 0.003364270 0.003218371 0.003738601 d1 0.003564965 0.010984234 0.003307337 0.003163906 0.003675332 d2 0.003364270 0.003307337 0.010899929 0.002985790 0.003468425 d3 0.003218371 0.003163906 0.002985790 0.010959203 0.003318008 d4 0.003738601 0.003675332 0.003468425 0.003318008 0.011846244 gender2 -0.005149393 -0.005062250 -0.004777264 -0.004570087 -0.005308813 gender2 d0 -0.005149393 d1 -0.005062250 d2 -0.004777264 d3 -0.004570087 d4 -0.005308813 gender2 0.007312139 > fit <- lm(resp~as.factor(dose)-1, data=IBScovars) > cf <- coef(fit) > vc <- vcov(fit) > doseVec <- 0:4 > fit <- fitMod(doseVec, cf, model="linInt", S=vc, type = "general") > plot(fit, CI=TRUE, plotData="meansCI") > vcov(fit) d0 d1 d2 d3 d4 d0 0.00819461 0.000000000 0.000000000 0.000000000 0.0000000 d1 0.00000000 0.007459197 0.000000000 0.000000000 0.0000000 d2 0.00000000 0.000000000 0.007757564 0.000000000 0.0000000 d3 0.00000000 0.000000000 0.000000000 0.008080796 0.0000000 d4 0.00000000 0.000000000 0.000000000 0.000000000 0.0079701 > fit <- lm(resp~as.factor(dose)+gender, data=IBScovars) > cf <- coef(fit)[2:5] > vc <- vcov(fit)[2:5,2:5] > doseVec <- 1:4 > fit <- fitMod(doseVec, cf, model="linInt", S=vc, type = "general", placAdj=TRUE) > vcov(fit) d1 d2 d3 d4 d1 0.015697662 0.008221459 0.008223928 0.008215124 d2 0.008221459 0.016014746 0.008246506 0.008208911 d3 0.008223928 0.008246506 0.016365819 0.008204394 d4 0.008215124 0.008208911 0.008204394 0.016212400 > plot(fit, CI=TRUE, plotData="meansCI") > predict(fit, predType = "effect-curve", se.fit=TRUE) $fit [1] 0.2845534 0.2965467 0.3501724 0.3479994 $se.fit [1] 0.1252903 0.1265494 0.1279290 0.1273279 > > proc.time() user system elapsed 1.50 0.32 1.81