## Data from Greco et al (1995)
d1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 10, 20, 50, 2, 2, 2, 2, 2, 5, 5, 5, 5, 5,10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 50, 50, 50, 50, 50)
d2 <- c(0, 0, 0, 0.2, 0.5, 1, 2, 5, 0, 0, 0, 0, 0, 0.2, 0.5, 1, 2, 5, 0.2, 0.5, 1, 2, 5, 0.2, 0.5, 1, 2, 5, 0.2, 0.5, 1, 2, 5, 0.2, 0.5, 1, 2, 5)
effect <- c(106.00, 99.20, 115.00, 79.20, 70.10, 49.00, 21.00, 3.83, 74.20, 71.50, 48.10, 30.90, 16.30, 76.30, 48.80, 44.50, 15.50, 3.21, 56.70, 47.50,
26.80, 16.90, 3.25, 46.70, 35.60, 21.50, 11.10, 2.94, 24.80, 21.60, 17.30, 7.78, 1.84, 13.60, 11.10, 6.43, 3.34, 0.89)
greco <- data.frame(d1, d2, effect)
## Fitting a generalized Loewe additivity model with common maximum
greco.m0a <- drm(effect~d1, data=greco, fct=genLoewe(fixed=c(NA,NA,NA,NA,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1,~1))
## Checking the model fit
plot(fitted(greco.m0a), residuals(greco.m0a))
## Fitting reduced model with lower limit c equal to 0
greco.m0b <- drm(effect~d1, data=greco, fct=genLoewe(fixed=c(NA,NA,0,NA,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1,~1))
## Comparing models by means of an F-test
anova(greco.m0b, greco.m0a)
## Looking at the summary output
summary(greco.m0b)
## Below commented out to make faster check on CRAN
## Fitting a generalized Loewe additivity model with different maxima
#greco.m0c <- drm(effect~d1, data=greco, fct=genLoewe2(fixed=c(NA,NA,NA,NA,NA,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1,~1))
## Checking the model fit
#plot(fitted(greco.m0c), residuals(greco.m0c))
## Looking at the summary output
#summary(greco.m0c) # suboptimal fit!
## Fitting the URSA model
#greco.m1 <- drm(effect~d1, data=greco, fct=ursa(fixed=c(NA,NA,0,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1))
#plot(fitted(greco.m1), residuals(greco.m1))
#summary(greco.m1)
## Fitting the URSA model using values from Greco et al (1995) (p. 364)
#greco.m2 <- drm(effect~d1, data=greco, fct=ursa(fixed=c(NA,NA,0,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1), start=c(-1.05, -2.04,95.1,11.1,1.07,0.52))
# same fit as above
## Adjusting for variance heterogeneity
#greco.m3 <- boxcox(greco.m2, method = "anova")
#plot(fitted(greco.m3), residuals(greco.m3)) # improved, not great
#summary(greco.m3)
## Getting closer to Greco's estimates
#greco.m4<-drm(effect~d1, data=greco, fct=ursa(fixed=c(NA,NA,0,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1), bcVal=0) # not perfect
#greco.m4 <- drm(effect~d1, data=greco, fct=ursa(fixed=c(NA,NA,0,NA,NA,NA,NA)), pmodels=list(~1,~1,~1,~I(1/d1)-1,~I(1/d2)-1,~1), bcVal=0,
#start=coef(greco.m3))
#plot(fitted(greco.m4), residuals(greco.m4)) # better!
#summary(greco.m4)
## Fitting the Bliss independence model (with common maximal response)
#greco.bliss.m1 <- drm(effect~d1+d2, data=greco, fct=genBliss())
#summary(greco.bliss.m1)
# huge estimated f parameters!
## Fitting the generalized Bliss independence model (with different maximal responses)
#greco.bliss.m2 <- drm(effect~d1+d2, data=greco, fct=genBliss())
#summary(greco.bliss.m2)
# stable fit?
#greco.bliss.m3 <- drm(effect~d1+d2, data=greco, fct=genBliss2(fixed=c(rep(NA,7),1,1)))
#summary(greco.bliss.m2)
# looks better
## Fitting the actimL model (needs some tweaking with the starting values)
#greco.actimL.m1 <- drm(effect~d1+d2, data=greco, fct=actimL(fixed=c(rep(NA,7),1,1,NA,1,rep(NA,3))))
#greco.actimL.m2 <- drm(effect~d1+d2, data=greco, fct=actimL(fixed=c(rep(NA,7),1,1,NA,1,rep(NA,3))), start=coef(greco.actimL.m1)+rnorm(11,0,0.1))
## perhaps some more tweaking with the starting values is needed to get the standard errors
#summary(greco.actimL.m2)
Run the code above in your browser using DataLab