Migra as vignettes e códigos fonte para o diretório int/doc

inst/doc/Ovelhas.R

+## ----setup, include=FALSE------------------------------------------------
+source("_setup.R")
+
+## ---- echo = FALSE, include=FALSE----------------------------------------
+##### Carregamento e tratamento inicial dos dados
+
+dados <- read.csv('../data-raw/Dadoscomp.csv',sep=',')
+dados$tratamento <- factor(dados$tratamento)
+levels(dados$tratamento) <- c('Observ', 'Observ + Interv')
+dados$linhagem <- factor(dados$linhagem)
+levels(dados$linhagem) <- c('Pouco reativo', 'Muito reativo')
+dados2 <- dados[1:38,c(1:5,30:53)]
+dados3 <- dados[1:38,30:53] ### Somente as mudanças de postura durante a intervenção.
+r1 <- rowSums(dados3[,1:3])-1
+table(r1)
+r2 <- rowSums(dados3[,4:6])-1
+table(r2)
+r4 <- rowSums(dados3[,12:16])-1
+table(r4)
+d2 <- rowSums(dados[1:38,57:59])-1
+table(d2)
+
+datapost1 <- data.frame(r1, dados2[ ,c('tratamento', 'linhagem')])
+names(datapost1)[1] <- 'Nposturas'
+datapost2 <- data.frame(r2, dados2[ ,c('tratamento', 'linhagem')])
+names(datapost2)[1] <- 'Nposturas'
+datapost4 <- data.frame(r4, dados2[ ,c('tratamento', 'linhagem')])
+names(datapost4)[1] <- 'Nposturas'
+
+##### Pacotes requeridos
+require('lmtest')
+require('boot')
+require('car')
+require('RColorBrewer')
+require('sandwich')
+require('hnp')
+require('knitr')
+
+
+## ------------------------------------------------------------------------
+str(datapost4)
+summary(datapost4)
+
+## ------------------------------------------------------------------------
+tab <- data.frame(with(datapost4, table(tratamento, Nposturas)))
+
+myColours <- brewer.pal(2,"Blues")
+
+my.settings <- list(
+  superpose.polygon=list(col=myColours[2:5], border="transparent"),
+  strip.background=list(col=myColours[6]),
+  strip.border=list(col="black")
+)
+
+bwplot(Nposturas ~ linhagem | tratamento,
+         data=datapost4, 
+         main="Mudanças de postura vs tratamento e linhagem", 
+         xlab="Linhagem", ylab="Frequência")
+
+### Variância e média amostrais por tratamento e linhagem.
+
+mdp <- aggregate(Nposturas ~ tratamento:linhagem, datapost4, function(x) c(mean = mean(x), var = var(x)))
+mdp
+
+
+## ------------------------------------------------------------------------
+ajusteps <- glm(Nposturas ~ tratamento + linhagem, data=datapost4, family=poisson)
+summary(ajusteps)
+
+### Avaliando possível efeito de interação.
+ajustepsint <- glm(Nposturas ~ tratamento * linhagem, data=datapost4, family=poisson)
+anova(ajusteps,ajustepsint, test = 'Chisq')
+
+exp(coef(ajusteps)[2])
+### Estima-se que a taxa de variação de postura no grupo com intervenção seja aproximadamente
+### a metade em relação ao grupo sem intervenção.
+
+##### Estimação do parâmetro de dispersão. 
+X2 <- sum(resid(ajusteps,type='pearson')**2)
+phichap <- X2/ajusteps$df.residual
+phichap
+
+## ------------------------------------------------------------------------
+##### Diagnóstico do modelo - gráficos padrão do R.
+par(mfrow=c(2,2))
+plot(ajusteps, pch = 20, cex = 1.25)
+
+## ---- echo = FALSE-------------------------------------------------------
+envelope=function(modelo){
+  dados=na.omit(modelo$data)
+  nsim=100
+  n=modelo$df.null+1
+  r1=sort(rstandard(modelo,type='deviance'))
+  m1=matrix(0,nrow=n,ncol=nsim)
+  a2=simulate(modelo,nsim=nsim)
+  
+  for (i in 1:nsim){
+    dados$y=a2[,i]
+    aj=update(modelo,y~.,data=dados)
+    m1[,i]=sort(rstandard(aj,type='deviance'))}
+  
+  li=apply(m1,1,quantile,0.025)
+  m=apply(m1,1,quantile,0.5)
+  ls=apply(m1,1,quantile,0.975)
+  
+  quantis=qnorm((1:n-0.5)/n)
+  
+  plot(rep(quantis,2),c(li,ls),type='n',xlab='Percentil da N(0,1)',ylab='Resíduos')
+  title('Gráfico Normal de Probabilidades')
+  lines(quantis,li,type='l')
+  lines(quantis,m,type='l',lty=2)
+  lines(quantis,ls,type='l')
+  points(quantis,r1,pch=16,cex=0.75)
+}
+
+
+## ------------------------------------------------------------------------
+##### Gráfico quantil-quantil com envelopes simulados.
+par(mfrow=c(1,1))
+envelope(ajusteps)
+
+## ------------------------------------------------------------------------
+
+### Modelo quasi poisson (V(mu) = mu).
+ajuste12 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = 'quasipoisson')
+# summary(ajuste12)
+
+### Forma alternativa de declarar o Modelo quase-poisson (V(mu) = mu).
+ajuste13 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu', link='log'))
+# summary(ajuste13)
+
+### Modelo de quase-verossimilhança (V(mu) = mu^2).
+ajuste14 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu^2', link='log'))
+summary(ajuste14)
+
+### Gráficos de diagnóstico para o modelo de quase-verossimilhança.
+par(mfrow = c(2,2))
+plot(ajuste14, pch = 20, cex = 1.25)
+
+### Gráficos quantil-quantil para os resíduos dos modelos Poisson e de Quase-Verossimilhança.
+par(mfrow=c(1,2))
+qqnorm(resid(ajusteps,type='deviance'), pch = 20, main = 'Poisson', las = 1)
+qqline(resid(ajusteps,type='deviance'))
+qqnorm(resid(ajuste14,type='deviance'), pch = 20, main = 'QL', las = 1)
+qqline(resid(ajuste14,type='deviance'))
+
+## ------------------------------------------------------------------------
+
+estrb <- coeftest(ajusteps, vcov=sandwich) ### Estimador sanduíche
+estrb ### Estimação robusta.
+
+### Usando bootstrap (R=999 simulações)
+ajusteboot <- Boot(ajusteps)
+summary(ajusteboot) ### Resultados obtidos via bootstrap.
+
+
+
+## ------------------------------------------------------------------------
+
+erroz <- rbind(summary(ajusteps)$coefficients[2,2:3], summary(ajuste13)$coefficients[2,2:3], 
+               summary(ajuste14)$coefficients[2,2:3], estrb[2,2:3], c(summary(ajusteboot)[2,4],
+               mean(ajusteboot$t[,2]/summary(ajusteboot)[2,4])))
+
+ics <- rbind(confint.default(ajusteps)[2,],confint.default(ajuste13)[2,], confint.default(ajuste14)[2,], 
+estrb[2,1] + c(-1.96,1.96) * estrb[2,2], mean(ajusteboot$t[,2])+c(-1.96,1.96)*summary(ajusteboot)[2,4]) 
+
+quadres <- cbind(erroz, ics)
+
+rownames(quadres) <- c('Poisson', 'Quasi(mu)', 'Quasi(mu^2)', 'Robusto (sanduiche)', 'Bootstrap')
+
+### Quadro resumo para as estimativas produzidas pelos quatro modelos para o efeito de intervenção.
+kable(quadres, format = "markdown", caption = "Comparativo dos modelos ajustados")
+
+
+## ------------------------------------------------------------------------
+dadosexclud <- datapost4[-c(8,18,28),]
+
+### Ajustando o modelo Poisson sem as três observações.
+ajusteexcludpoi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = poisson)
+
+### Estimativa do parâmetro de dispersão.
+sum(resid(ajusteexcludpoi, type = 'pearson')**2)/ajusteexcludpoi$df.residual 
+
+### Gráficos de diagnóstico.
+par(mfrow=c(2,2))
+plot(ajusteexcludpoi)
+
+### Agora, o modelo quase-poisson
+ajusteexcludquasi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = quasi(variance = 'mu', link='log'))
+
+### Estimativas produzidas pelo modelo quasipoisson com e sem as três observações.
+
+c1 <- compareCoefs(ajusteps, ajuste13, ajusteexcludpoi, ajusteexcludquasi, print = FALSE)
+colnames(c1) <- c('Est. Poisson c/out','E.P. Poisson c/out',
+                  'Est. Quasi c/out','E.P. Quasi c/out',
+                   'Est. Poisson s/out','E.P. Poisson s/out',
+                  'Est. Quasi s/out','E.P. Quasi s/out') 
+kable(round(c1,3), align = 'c')
+
+### Efeito da intervenção desconsiderando as três observações.
+exp(coef(ajusteexcludpoi)[2])
+
+### O efeito de intervenção aumenta, e torna-se mais significativo, mediante exclusão dos outliers.
+
+

inst/doc/Sinistros.R

+## ----setup, include=FALSE------------------------------------------------
+source("_setup.R")
+
+## ---- echo = FALSE, include=FALSE----------------------------------------
+
+require(lattice)
+require(MASS)
+require(effects)
+require(knitr)
+
+dados <- read.csv('../data-raw/Baseauto2.csv', header=TRUE, sep=',', dec = ',')
+
+options(scipen = 3) ### Era zero.
+
+##### Preparando os dados.
+names(dados) <- c('Bonus','Tipo','Estado','Idade','Sexo','Civil','Valor','Exposicao','Sinistros')
+dados <- dados[,c('Idade','Sexo','Valor','Exposicao','Sinistros')]
+dados <- dados[which(dados$Valor != 0),]
+dados$lexpo <- log(dados$Exposicao)
+dados$Valor <- dados$Valor/1000
+
+
+## ------------------------------------------------------------------------
+head(dados, 10) ### Dez primeiras linhas da base.
+str(dados)
+
+## ------------------------------------------------------------------------
+table(dados$Sinistros) ### Distribuição do números de sinistros.
+taxageral <- sum(dados$Sinistros)/sum(dados$Exposicao); taxageral ### Taxa de sinistros na amostra.
+
+tab <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo, data = dados, sum)
+taxa <- with(tab, Sinistros/Exposicao)
+tab <- cbind(tab, taxa)
+### Distribuição do número de sinistros por sexo.
+kable(tab, align = 'c', caption = '**Taxa de sinistros segundo Sexo**')
+
+dados$idadecat <- cut(dados$Idade, breaks=c(18,30,60, 92), include.lowest = T)
+tab <- aggregate(cbind(Exposicao, Sinistros) ~ idadecat, data = dados, sum)
+taxa <- with(tab, Sinistros/Exposicao)
+tab <- cbind(tab, taxa)
+### Distribuição do número de sinistros por sexo.
+kable(tab, align = 'c', caption = '**Taxa de sinistros segundo Idade**')
+
+tabidsex <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo + idadecat, data = dados, sum)
+taxa <- with(tabidsex, Sinistros/Exposicao)
+tabidsex <- cbind(tabidsex, taxa)
+### Distribuição do número de sinistros por idade e sexo.
+kable(tabidsex, align = 'c', caption = '**Taxa de sinistros segundo Sexo e Idade**')
+
+
+## ------------------------------------------------------------------------
+dados <- na.omit(dados)
+ajusteps <- glm(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+offset(log(Exposicao)), data = dados, family=poisson)
+summary(ajusteps)
+
+### Estimação do parâmetro de dispersão. 
+
+X2 <- sum(resid(ajusteps,type='pearson')**2)
+phichap <- X2/ajusteps$df.residual
+phichap ### Indicador de superdispersão.
+
+## ---- echo = FALSE-------------------------------------------------------
+envelope=function(modelo){
+  dados=na.omit(modelo$data)
+  nsim=100
+  n=modelo$df.null+1
+  r1=sort(rstandard(modelo,type='deviance'))
+  m1=matrix(0,nrow=n,ncol=nsim)
+  a2=simulate(modelo,nsim=nsim)
+  
+  for (i in 1:nsim){
+    dados$y=a2[,i]
+    aj=update(modelo,y~.,data=dados)
+    m1[,i]=sort(rstandard(aj,type='deviance'))}
+  
+  li=apply(m1,1,quantile,0.025)
+  m=apply(m1,1,quantile,0.5)
+  ls=apply(m1,1,quantile,0.975)
+  
+  quantis=qnorm((1:n-0.5)/n)
+  
+  plot(rep(quantis,2),c(li,ls),type='n',xlab='Percentil da N(0,1)',ylab='Resíduos')
+  title('Gráfico Normal de Probabilidades')
+  lines(quantis,li,type='l')
+  lines(quantis,m,type='l',lty=2)
+  lines(quantis,ls,type='l')
+  points(quantis,r1,pch=16,cex=0.75)
+}
+
+## ------------------------------------------------------------------------
+##### Diagnóstico do modelo - gráficos.
+par(mfrow=c(2,2))
+plot(ajusteps)
+
+par(mfrow=c(1,1))
+envelope(ajusteps)
+
+## ------------------------------------------------------------------------
+ajusteps2 <- glm(Sinistros ~ Sexo + Idade +I(Idade**2) + Valor + log(Exposicao), data = dados, family=poisson)
+summary(ajusteps2)
+anova(ajusteps, ajusteps2, test = 'Chisq')
+
+## ------------------------------------------------------------------------
+ajustenb2 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+log(Exposicao),data= dados)
+summary(ajustenb2) 
+
+
+## ------------------------------------------------------------------------
+##### Diagnóstico do modelo - gráficos.
+par(mfrow=c(2,2))
+plot(ajustenb2)
+
+## ---- echo = FALSE-------------------------------------------------------
+dadosnb3 <- dados[,c('Sexo','Idade','Valor','Exposicao','Sinistros')]
+dadosnb3$lexpo <- log(dados$Exposicao)
+ajustenb2 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+lexpo,data= dadosnb3)
+envelope=function(modelo){
+  dados=na.omit(dadosnb3)
+  nsim=100
+  n=modelo$df.null+1
+  r1=sort(rstandard(modelo,type='deviance'))
+  m1=matrix(0,nrow=n,ncol=nsim)
+  a2=simulate(modelo,nsim=nsim)
+  
+  for (i in 1:nsim){
+    dados$y=a2[,i]
+    aj=update(modelo,y~.,data=dados)
+    m1[,i]=sort(rstandard(aj,type='deviance'))}
+  
+  li=apply(m1,1,quantile,0.025)
+  m=apply(m1,1,quantile,0.5)
+  ls=apply(m1,1,quantile,0.975)
+  
+  quantis=qnorm((1:n-0.5)/n)
+  
+  plot(rep(quantis,2),c(li,ls),type='n',xlab='Percentil da N(0,1)',ylab='Resíduos')
+  title('Gráfico Normal de Probabilidades')
+  lines(quantis,li,type='l')
+  lines(quantis,m,type='l',lty=2)
+  lines(quantis,ls,type='l')
+  points(quantis,r1,pch=16,cex=0.75)
+}
+
+## ------------------------------------------------------------------------
+par(mfrow=c(1,1))
+envelope(ajustenb2)
+
+## ------------------------------------------------------------------------
+intervalos <- confint(ajustenb2)
+estimat <- cbind(ajustenb2$coefficients, intervalos)
+colnames(estimat)[1] <- 'Estimativa pontual'
+
+### Quadro de estimativas
+kable(round(estimat, 5), align = 'c', caption = '**Estimativas pontuais e intervalos de confiança - Modelo Binomial Negativo**')
+
+## ---- echo = FALSE-------------------------------------------------------
+
+
+## ------------------------------------------------------------------------
+
+efeitos <- allEffects(ajustenb2, given.values=c(lexpo=0))
+trellis.par.set(list(axis.text = list(cex = 1.2))) 
+
+plot(efeitos[[2]], type='response',main=list(
+    label="Taxa de sinistros vs. Idade",
+    cex=1.5),
+    xlab=list(
+        label="Idade (anos)",
+        cex=1.5),
+    ylab=list(
+        label="Taxa de sinistros",
+        cex=1.5))
+
+plot(efeitos[[1]], type='response',main=list(
+    label="Taxa de sinistros vs. Sexo",
+    cex=1.5),
+    xlab=list(
+        label="Sexo",
+        cex=1.5),
+    ylab=list(
+        label="Taxa de  sinistros",
+        cex=1.5))
+
+plot(efeitos[[4]], type='response',main=list(
+    label="Taxa de sinistros vs. Valor do automóvel",
+    cex=1.5),
+    xlab=list(
+        label="Valor (x1000 reais)",
+        cex=1.5),
+    ylab=list(
+        label="Taxa de sinistros",
+        cex=1.5))
+
+
+## ---- echo = FALSE, include=F--------------------------------------------
+
+## Frequências ajustadas pelas duas distribuições, com e sem covariaveis.
+
+### Poisson sem ajuste de covariáveis.
+n <- nrow(dados)
+mediasin <- mean(dados$Sinistros)
+freqps <- round(n*dpois(0:10,mediasin))
+
+
+### Poisson com covariaveis
+pred1 <- predict(ajusteps2,type='response')
+intervalo <- 0:10
+matprob <- matrix(0, nrow=nrow(dados), ncol=length(intervalo))
+probpois <- function(interv, taxa)
+    dpois(intervalo,taxa)
+for(i in 1:nrow(dados))
+    matprob[i,] <- probpois(interv = intervalo, taxa = pred1[i])
+pbarra <- colMeans(matprob)
+freqpsaj <- round(n*pbarra)
+
+
+
+### Binomial negativa sem covariaveis.
+ajustenb <- glm.nb(Sinistros ~ 1,data=dados)
+
+media <- exp(coefficients(ajustenb))
+shape <- ajustenb$theta
+freqbn <- round(n*dnbinom(0:10, mu = media, size = shape))
+
+
+### Binomial negativa com covariaveis
+pred2 <- predict(ajustenb2,type='response')
+
+intervalo <- 0:10
+matprob <- matrix(0,nrow=nrow(dados),ncol=length(intervalo))
+probnb <- function(interv, media, shape)
+    dnbinom(intervalo, mu = media, size = shape)
+for(i in 1:nrow(dados))
+    matprob[i,] <- probnb(interv = intervalo, media = pred2[i], shape = ajustenb2$theta)
+pbarra <- colMeans(matprob)
+frebnaj <- round(n*pbarra)
+ams <- c(table(dados$Sinistros), rep(0,5))
+matfreq <- rbind(ams, freqps, freqpsaj, freqbn, frebnaj)
+colnames(matfreq) <- 0:10
+rownames(matfreq) <- c('Amostra', 'Poisson não ajustada por covariáveis', 'Poisson ajustada por covariáveis', 
+                       'BN não ajustada por covariáveis', 'BN ajustada por covariáveis')
+
+## ---- results = 'markup'-------------------------------------------------
+kable(matfreq, format = "markdown", caption = "Frequências amostrais e frequências ajustadas para o número de sinistros")
+

vignettes/Ovelhas.Rmd

@@ -29,7 +29,7 @@ corporal do animal ao longo do período de observação (3 minutos).
 ```{r, echo = FALSE, include=FALSE}
 ##### Carregamento e tratamento inicial dos dados
 
-dados <- read.csv('Dadoscomp.csv',sep=',')
+dados <- read.csv('../data-raw/Dadoscomp.csv',sep=',')
 dados$tratamento <- factor(dados$tratamento)
 levels(dados$tratamento) <- c('Observ', 'Observ + Interv')
 dados$linhagem <- factor(dados$linhagem)
@@ -60,6 +60,7 @@ require('RColorBrewer')
 require('sandwich')
 require('hnp')
 require('knitr')
+
 ```
 
 

vignettes/Sinistros.Rmd

@@ -34,7 +34,7 @@ require(MASS)
 require(effects)
 require(knitr)
 
-dados <- read.csv('Baseauto2.csv', header=TRUE, sep=',', dec = ',')
+dados <- read.csv('../data-raw/Baseauto2.csv', header=TRUE, sep=',', dec = ',')
 
 options(scipen = 3) ### Era zero.
 
@@ -47,14 +47,14 @@ dados$Valor <- dados$Valor/1000
 
 ```
 
-### Verificação do conteúdo e a estrutura dos dados.
+### Verificação do conteúdo e a estrutura dos dados
 
 ```{r}
 head(dados, 10) ### Dez primeiras linhas da base.
 str(dados)
 ```
 
-### Análise descritiva da distribuição do número de sinistros. 
+### Análise descritiva da distribuição do número de sinistros
 
 ```{r}
 table(dados$Sinistros) ### Distribuição do números de sinistros.
@@ -81,7 +81,7 @@ kable(tabidsex, align = 'c', caption = '**Taxa de sinistros segundo Sexo e Idade
 
 ```
 
-## Regressão usando o modelo log-linear poisson.
+## Regressão usando o modelo log-linear poisson
 
 ```{r}
 dados <- na.omit(dados)
@@ -125,7 +125,7 @@ envelope=function(modelo){
 ```
 
 
-### Diagnóstico do ajuste (gráficos).
+### Diagnóstico do ajuste (gráficos)
 
 ```{r}
 ##### Diagnóstico do modelo - gráficos.
@@ -146,7 +146,7 @@ anova(ajusteps, ajusteps2, test = 'Chisq')
 ```
 
 
-## Regressão usando a distribuição binomial negativa.
+## Regressão usando a distribuição binomial negativa
 
 ```{r}
 ajustenb2 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+log(Exposicao),data= dados)
@@ -155,7 +155,7 @@ summary(ajustenb2)
 ```
 
 
-### Diagnóstico do ajuste (gráficos).
+### Diagnóstico do ajuste (gráficos)
 
 ```{r}
 ##### Diagnóstico do modelo - gráficos.
@@ -316,14 +316,3 @@ rownames(matfreq) <- c('Amostra', 'Poisson não ajustada por covariáveis', 'Poi
 ```{r, results = 'markup'}
 kable(matfreq, format = "markdown", caption = "Frequências amostrais e frequências ajustadas para o número de sinistros")
 ```
-
-
-
-
-
-
-
-
-
-
-