# Bicycle example from Box/Hunter/Hunter
# y: time to cycle up a hill
# A: seat up/down
# B: dynamo off/on
# C: handlebars up/down
# D: gear low/medium
# E: raincoat on/off
# F: breakfast yes/no
# G: tyres hard/soft

library(FrF2)

planf <- FrF2(nruns=2^7,nfactors=7) # too many?

plan <- FrF2(nruns=8,nfactors=7,generators=c("AB","AC","BC","ABC"),randomize=FALSE)
plan
y <- c(69,52,60,83,71,50,59,88)
plan <- add.response(plan,y)

fit <- lm(y~.,data=plan)
summary(fit)

DanielPlot(fit)

lenth <- function(lmobj,alpha=0.05)
{
  abseffects <- abs(2*lmobj$coeff)[-1]
  medabseffects <- median(abseffects)
  medabseffects
  s0 <- 1.5*medabseffects
  keepeffects <- abseffects[abseffects< 2.5*s0]
  PSE <- 1.5*median(keepeffects)
  signlimit <-qt(1-alpha/2,length(abseffects)/3)*PSE
  return(list("PSE"=PSE,"Signlimit"=signlimit,"Number of sign"=sum(abseffects>signlimit)))
}

lenthres <- lenth(fit)
lenthres
# see that factors with effects larger than lenthres$Signlimit should be assumed
# to have effect different from zero
effects <- 2*fit$coeff
effects>lenthres$Signlimit

# How to make a relatively ok Paretoplot for the effects
barplot(sort(abs(effects[-1]),decreasing=FALSE),las=1,horiz=TRUE,cex.names=1.0)
abline(v=lenthres$Signlimit,col=2,lwd=2)