% COMANDOS EM R REFERENTES AO EXEMPLO DELIVERY

  delivery = read.table("delivery.txt", header=TRUE)
  attach(delivery)

  require(robustbase)

  adjbox(tempo, ylab="Tempo Gasto")
  adjbox(ncaixas, ylab="Número de Caixas")
  adjbox(distancia, ylab="Distância Percorrida")

  plot(ncaixas, tempo, ylab="Tempo Gasto", xlab="Número de Caixas", pch=16,
  cex=2,cex.lab=1.5,cex.axis=1.5)
  lines(smooth.spline(ncaixas, tempo, df=3),lwd=2)

  plot(distancia, tempo, ylab="Tempo Gasto", xlab="Distância Percorrida", pch=16,
  cex=2,cex.lab=1.5,cex.axis=1.5)
  lines(smooth.spline(distancia, tempo, df=3),lwd=2)

  ajuste1.delivery = lm(tempo ~ ncaixas + distancia)
  summary(ajuste1.delivery)

  fit.model = ajuste1.delivery
  source("envel_norm")
  source("diag_resid_norm")
  source("diag_cook_norm")

  ajuste2.delivery = lm(tempo ~ ncaixas + distancia, subset=-9)
  summary(ajuste2.delivery)

% GRÁFICO DA VARIÁVEL ADICIONADA ncaixas

  ajuste3.delivery = lm(tempo ~ distancia)
  r = resid(ajuste3.delivery)
  ajuste4.delivery = lm(ncaixas ~ distancia)
  v = resid(ajuste4.delivery)
  plot(v,r,pch=16, xlab="Resíduo v", ylab="Resíduo r",
  cex=2,cex.lab=1.5,cex.axis=1.5,
  main="Variável Adicionada ncaixas",cex.main=1.5)
  lines(smooth.spline(v,r,df=3),lwd=2)

% GRÁFICO DA VARIÁVEL ADICIONADA distancia

  ajuste5.delivery = lm(tempo ~ ncaixas)
  r = resid(ajuste5.delivery)
  ajuste6.delivery = lm(distancia ~ ncaixas)
  v = resid(ajuste6.delivery)
  plot(v,r,pch=16, xlab="Resíduo v", ylab="Resíduo r",
  cex=2,cex.lab=1.5,cex.axis=1.5,
  main="Variável Adicionada distancia",cex.main=1.5)
  lines(smooth.spline(v,r,df=3),lwd=2)

% AJUSTE DE TERMO QUADRÁTICO

  ndistancia = distancia/1000

  ajuste7.delivery = lm(tempo ~ ncaixas + ndistancia + I(ndistancia^2))
  summary(ajuste7.delivery)

  ajuste8.delivery = lm(tempo ~ ncaixas + I(ndistancia^2))  
  summary(ajuste8.delivery)
  fit.model = ajuste8.delivery
  
  source("envel_norm")
  source("diag_resid_norm")
  source("diag_cook_norm")

% AJUSTE NO GAMLSS

  require(gamlss) 
  ajuste9.delivery = gamlss(tempo ~ ncaixas + I(ndistancia^2), family=NO(sigma.link="identity"))
  summary(ajuste9.delivery)
  plot(ajuste9.delivery)
  wp(ajuste9.delivery,cex=2,cex.lab=1.5,cex.axis=1.5)