% COMANDOS EM R REFERENTES AO EXEMPLO ROTIFERS
  
  rotifers = read.table("rotifers.txt", header=TRUE)
  attach(rotifers)
  summary(rotifers)

  dens0 = split(dens,espec)$"0"
  dens1 = split(dens,espec)$"1"

  susp0 = split(susp,espec)$"0"
  susp1 = split(susp,espec)$"1"

  exp0 = split(exp,espec)$"0"
  exp1 = split(exp,espec)$"1"
 
  freq0 = susp0/exp0
  freq1 = susp1/exp1

  plot(dens0,freq0, type="n", xlab="Densidade Relativa",
  ylab="Frequência Relativa", main="Espécie 0", pch=16,cex=2, cex.lab=1.5, cex.axis=1.5,cex.main=1.5)
  lines(smooth.spline(dens0,freq0,df=3), lwd=1.5)
  points(dens0,freq0,pch=16,cex=2)

  plot(dens1,freq1, type="n", xlab="Densidade Relativa",
  ylab="Frequência Relativa", main="Espécie 1", pch=16,cex=2, cex.lab=1.5, cex.axis=1.5, cex.main=1.5)
  lines(smooth.spline(dens1,freq1,df=3),lwd=1.5)
  points(dens1,freq1,pch=16,cex=2)

  y0.rotifers = cbind(susp0, exp0-susp0)
  y1.rotifers = cbind(susp1, exp1-susp1)

  fit1.esp0 = glm(y0.rotifers ~ dens0, family=binomial)
  summary(fit1.esp0)
  ntot = exp0
  fit.model = fit1.esp0
  source("envelr_bino")
  
  fit1.esp1 = glm(y1.rotifers ~ dens1, family=binomial)
  summary(fit1.esp1)
  ntot = exp1
  fit.model = fit1.esp1
  source("envelr_bino")

  require(gamlss)
   
  fit2.esp1 = gamlss(y1.rotifers ~dens1, family=BB)
  summary(fit2.esp1)
  plot(fit2.esp1)
  rqres.plot(fit2.esp1, howmany=6, type="wp",cex=2, cex.lab=1.5, cex.axis=1.5)
 
  fit2.esp0 = gamlss(y0.rotifers ~dens0, family=BB)
  summary(fit2.esp0)
  plot(fit2.esp0)
  rqres.plot(fit2.esp0, howmany=6, type="wp",cex=2, cex.lab=1.5, cex.axis=1.5)

  fit3.esp0 = gamlss(y0.rotifers ~dens0, ~dens0, family=BB)
  summary(fit3.esp0)
  plot(fit3.esp0)
  rqres.plot(fit3.esp0, howmany=6, type="wp",cex=2, cex.lab=1.5, cex.axis=1.5)

  fit4.esp0 = gamlss(y0.rotifers ~dens0 + I(dens0^2), ~dens0,
  family=BB)
  summary(fit4.esp0)
  plot(fit4.esp0)
  rqres.plot(fit4.esp0, howmany=6, type="wp",cex=2, cex.lab=1.5, cex.axis=1.5)

% CURVAS AJUSTADAS

  dens = seq(1.015,1.075,0.001)
  pred1 = -90.19 + 86.86*dens
  pred1 = exp(pred1)/(1+exp(pred1))
  plot(dens,pred1, ylab="Probabilidade de Ficar Suspenso",
  xlab="Densidade Relativa",ylim=c(0,1),xlim=c(1.014,1.076),
  pch=16,type="n",cex=2, cex.lab=1.5, cex.axis=1.5,main="Espécie 1",cex.main=1.5)
  lines(dens,pred1,lwd=2)
  points(dens1,freq1,pch=16,cex=2)

  dens = seq(1.015,1.075,0.001)
  pred0 = 2449.5 - 4790.5*dens + 2340.1*I(dens^2)
  pred0 = exp(pred0)/(1+exp(pred0))
  plot(dens,pred0, ylab="Probabilidade de Ficar Suspenso",
  xlab="Densidade Relativa",ylim=c(0,1), xlim=c(1.014,1.076),
  pch=16,type="n", main="Espécie 0",cex=2, cex.lab=1.5, cex.axis=1.5, cex.main=1.5)
  lines(dens,pred0,lwd=1.5)
  points(dens0,frq0,pch=16,cex=2)