% COMANDOS EM R REFERENTES AO EXEMPLO ATAQUES EPILÉPTICOS

  ataques  <- read.table("ataques.txt")
  ataques1 <- read.table("ataques1.txt")


% USANDO gee
  require(gee)

  fit1.ataques <- gee(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas)), data=ataques, id=ind, 
  family=poisson, corstr="exchangeable")
  summary(fit1.ataques)

  fit2.ataques <- gee(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas)), data=ataques1, id=ind, 
  family=poisson, corstr="exchangeable")
  summary(fit2.ataques)

  source("diag_gee_pois")
  analise_diag_poisson <- diag_gee_poisson(fit1.ataques, ataques, umaJanela=F, selGraficos=c(0,0,1,0), identifica=c(0,0,3,0))

  source("envel_gee_pois")
  envelope_gee_poisson(fit1.ataques, ataques, analise_diag_poisson, opcaoEnvelope="Normal")

  analise_diag_poisson <- diag_gee_poisson(fit2.ataques, ataques1, umaJanela=F, selGraficos=c(0,0,1,0), identifica=c(0,0,3,0))
  envelope_gee_poisson(fit2.ataques, ataques1, analise_diag_poisson, opcaoEnvelope="Normal")


% CÁLCULO DO QIC

  require(MuMIn)
  model.sel(fit1.ataques, fit2.ataques, rank = QIC)

% USANDO geepack
  require(geepack)

  fit3.ataques <- geeglm(nataques~factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas)), data=ataques, id=ind, 
family=poisson, corstr="exchangeable")
  summary(fit3.ataques)

  fit4.ataques <- geeglm(nataques~factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas)), data=ataques1, id=ind, 
family=poisson, corstr="exchangeable")
  summary(fit4.ataques)

% USANDO gamlss

  require(gamlss)

  fit5.ataques <- gamlss(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas))
+ random(as.factor(ind)), data=ataques, family = PO)  
  summary(fit5.ataques)
  plot(fit5.ataques)
  wp(fit5.ataques)

  fit6.ataques <- gamlss(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas))
+ random(as.factor(ind)), data=ataques1, family = PO)  
  summary(fit6.ataques)
  plot(fit6.ataques)

  fit7.ataques <- gamlss(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas))
+ random(as.factor(ind)), data=ataques, family = NBI)  
  summary(fit7.ataques)
  plot(fit7.ataques)
  wp(fit7.ataques)

  fit8.ataques = gamlss(nataques ~ factor(tratamento) + factor(periodo) + factor(tratamento)*factor(periodo) + offset(log(nsemanas)) +
    re(random=~1|ind),data=ataques, family = NBI) 
  summary(fit8.ataques)

  ranef0 = ranef(getSmo(fit8.ataques))$"(Intercept)"
  plot(density(ranef0))