#### experiment simulation
####  n=10, two conditions A and B
### parameters
nrep <- 10
lambdaA <- 50
b <- 1.5
lambdaB <- lambdaA*b

### simulations
xA <- rpois(n=nrep,lambda=lambdaA)
xB <- rpois(n=nrep, lambda=lambdaB)
### data table  
  ## vector describing the conditions 
  conds <- rep(c("A","B"),each=nrep)
  ### vector of observations
  Y <- c(xA,xB)
  tab <- cbind.data.frame(conds=conds,Y=Y)
  head(tab)
  ## visualization
  boxplot(Y~conds,data=tab)

### analysis with glm
myglm <- glm(Y~conds,family=poisson(),data=tab)

  ## anova
  anova(myglm,test="Chisq")
  ### residuals
  par(mar=c(4,4,4,4))
  plot(fitted(myglm),residuals(myglm,type="pearson"),
       xlab="Fitted values",ylab="residuals",
       col=as.factor(tab$conds),pch=19)
  
# fitted values are exrpessed in the original unit of the variable
  
par(mfrow=c(2,2))
plot(myglm)

# if a linear mdoel is used - wrong choice of modelisation
mylm=lm(Y~conds, data=tab)
plot(mylm)