diff --git a/inst/slides/images/graf1.pdf b/inst/slides/images/graf1.pdf new file mode 100644 index 0000000000000000000000000000000000000000..9cb41388fa25ad83e2b8643ccbda8ab6ec8235b5 Binary files /dev/null and b/inst/slides/images/graf1.pdf differ diff --git a/inst/slides/images/graf2.pdf b/inst/slides/images/graf2.pdf new file mode 100644 index 0000000000000000000000000000000000000000..998e9068ec962d67a2038d93ecd273adbd7782f8 Binary files /dev/null and b/inst/slides/images/graf2.pdf differ diff --git a/inst/slides/images/processos14.pdf b/inst/slides/images/processos14.pdf new file mode 100644 index 0000000000000000000000000000000000000000..fb02e99a3893f869fcbdd8d52ede6b4a2439ebf0 Binary files /dev/null and b/inst/slides/images/processos14.pdf differ diff --git a/vignettes/Ovelhas.Rmd b/vignettes/Ovelhas.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..d8b6cdd10dae2d696bc9d7d5f9d54f5685cf503e --- /dev/null +++ b/vignettes/Ovelhas.Rmd @@ -0,0 +1,240 @@ +--- +title: "Análise de Contagens - Quase-Verossimilhança" +author: > + Walmes M. Zeviani, + Eduardo E. Ribeiro Jr & + Cesar A. Taconeli +vignette: > + %\VignetteIndexEntry{Análise de Contagens - Quase-Verossimilhança} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include=FALSE} +source("_setup.R") +``` + +Dados referentes a um experimento delineado em blocos, com o objetivo de +investigar o efeito de uma intervenção, por parte do cuidador, no +comportamento de ovelhas. + +Para isso, foram considradas ovelhas de duas +linhagens distintas (pouco ou muito reativas), submetidas a dois tipos +diferentes de intervenção (observação ou observação + intervenção). + +A variável resposta aqui considerada é o número de mudanças na postura +corporal do animal ao longo do período de observação. + + +```{r, echo = FALSE, include=FALSE} +##### Carregamento e tratamento inicial dos dados +dados <- read.csv2('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') +``` + + +### Verificação do conteúdo e a estrutura dos dados. + +```{r} +str(datapost4) +summary(datapost4) +``` + + + +### Análise descritiva + +```{r} +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 + +``` + + +### Regressão poisson com estimação por máxima verossimilhança. + +```{r} +ajusteps <- glm(Nposturas ~ tratamento + linhagem, data=datapost4, family=poisson) +summary(ajusteps) + +exp(coef(ajusteps)[2]) +### Estima-se que a frequência média 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 ajuste (gráficos). + +```{r} +##### Diagnóstico do modelo - gráficos padrão do R. +par(mfrow=c(2,2)) +plot(ajusteps, pch = 20, cex = 1.25) +``` + +```{r, 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) +} + +``` + +```{r} +##### Gráfico quantil-quantil com envelopes simulados. +par(mfrow=c(1,1)) +envelope(ajusteps) +``` + + +### Ajustando modelos por quase-verossimilhança. + +```{r} + +### 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')) +``` + + +### Usando estimação robusta e bootstrap. + +```{r} + +estrb <- coeftest(ajusteps, vcov=sandwich) ### Estimador sanduíche +estrb + +### Usando bootstrap (R=999 simulações) +ajusteboot <- Boot(ajusteps) +plot(ajusteboot, index = 2) ### Distribuição bootstrap para o efeito de intervenção. +summary(ajusteboot) + + +``` + + +### Apanhado geral dos resultados. + +```{r} + +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. +kable(quadres, format = "markdown", caption = "Comparativo dos modelos ajustados") + +### Vamos avaliar o efeito das observações com maiores resíduos nos achados do modelo +dadosexclud <- datapost4[-c(8,18,28),] +ajusteexclud <- 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(ajuste13, ajusteexclud, print = FALSE) +colnames(c1) <- c('Est. com outliers','E.P. com outliers','Est. sem outliers','E.P. sem outliers') +kable(c1) + +### Efeito da intervenção desconsiderando as três observações. +exp(coef(ajusteexclud)[2]) + +### O efeito de intervenção aumenta, e torna-se mais significativo, mediante exclusão dos outliers. + +```` + + diff --git a/vignettes/Sinistros.Rmd b/vignettes/Sinistros.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..c6ec8708cbc379c47a1f9c87844116c8d07ec4da --- /dev/null +++ b/vignettes/Sinistros.Rmd @@ -0,0 +1,352 @@ +--- +title: "Análise de Contagens - Poisson e Binomial Negativa" +author: > + Walmes M. Zeviani, + Eduardo E. Ribeiro Jr & + Cesar A. Taconeli +vignette: > + %\VignetteIndexEntry{Análise de Contagens - Poisson e Binomial Negativa} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include=FALSE} +source("_setup.R") +``` + + +Dados referentes ao número de sinistros registrados por 16483 clientes de +uma seguradora de automóveis ao longo de um ano, contemplando as +seguintes variáveis: + +* **Sinistros**: Número de sinistros registrados; +* **Exposicao**: Período de cobertura do cliente durante o ano sob análise; +* **Tipo**: Tipo de veículo (hatch, monobloco, picape, sedan, wagon e suv); +* **Idade**: Idade do cliente (em anos); +* **Sexo**: M para masculino e F para feminino; +* **Civil**: Estado civil (Solteiro, Casado, Divorciado e Viuvo); +* **Valor**: Valor do veículo segurado (em reais). + + + +```{r, echo = FALSE, include=FALSE} + +require(lattice) +require(hnp) +require(MASS) +require(effects) + +# setwd("C:/Users/CCE/Dropbox/Backup Parana/Parana/Curso - dados discretos") +dados <- read.csv2('Baseauto2.csv', header=TRUE, sep=',') + +options(scipen = 3) ### Era zero. + +##### Preparando os dados. +names(dados) <- c('Bonus','Tipo','Estado','Idade','Sexo','Civil','Valor','Exposicao','Sinistros') +dados <- dados[,c('Tipo','Idade','Sexo','Civil','Valor','Exposicao','Sinistros')] +levels(dados$Tipo) <- c('outros','outros','outros','outros','hatch','outros','mono','picape','picape','picape','sedan','wagon','suv') +dados <- dados[-which(dados$Tipo=='outros'),] +dados$Tipo <- factor(dados$Tipo) +dados <- dados[which(dados$Valor != 0),] +dados$Civil <- relevel(dados$Civil, ref = 'Solteiro') +dados$lexpo <- log(dados$Exposicao) +dados$Valor <- dados$Valor/1000 + +``` + +### 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. + +```{r} +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); tab ### Distribuição do número de sinistros por 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); tab ### Distribuição do número de sinistros por sexo. + +tabidsex <- aggregate(cbind(Exposicao, Sinistros) ~ Sexo + idadecat, data = dados, sum) +Taxaidsex <- with(tabidsex, Sinistros/Exposicao) +tabidsex <- cbind(tabidsex, Taxaidsex); tabidsex ### Distribuição do número de sinistros por idade e sexo. + +tab <- aggregate(cbind(Exposicao, Sinistros) ~ Tipo, data = dados, sum) +taxa <- with(tab, Sinistros/Exposicao) +tab <- cbind(tab, taxa); tab ### Distribuição do número de sinistros por tipo de veículo. + +tab <- aggregate(cbind(Exposicao, Sinistros) ~ Civil, data = dados, sum) +taxa <- with(tab, Sinistros/Exposicao) +tab <- cbind(tab, taxa); tab ### Distribuição do número de sinistros por estado civil. + +dados$valorcat <- cut(dados$Valor, breaks=quantile(dados$Valor), include.lowest = T) +tab <- aggregate(cbind(Exposicao, Sinistros) ~ valorcat, data = dados, sum) +taxa <- with(tab, Sinistros/Exposicao) +tab <- cbind(tab, taxa); tab ### Distribuição do número de sinistros por valor do veículo. +``` + +## Regressão usando o modelo log-linear poisson. + +```{r} +dados <- na.omit(dados) +ajusteps <- glm(Sinistros ~ Tipo+Sexo+Idade+I(Idade**2)+Valor+Civil+offset(log(Exposicao)), data = dados, family=poisson) +summary(ajusteps) + +##### Estimação do parâmetro de dispersão. + +exp(coef(ajusteps)) + +X2 <- sum(resid(ajusteps,type='pearson')**2) +phichap <- X2/ajusteps$df.residual +phichap ### Indicador de superdispersão. +``` + +### Diagnóstico do ajuste (gráficos). + +```{r} +##### Diagnóstico do modelo - gráficos. +par(mfrow=c(2,2)) +plot(ajusteps) + + +par(mfrow=c(1,1)) +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) +} + +envelope(ajusteps) + +``` + + +### Re-ajuste do modelo associando um parâmetro ao termo offset (log-exposicao) + +```{r} +ajusteps2 <- glm(Sinistros ~ Tipo + Sexo + Idade +I(Idade**2) + Valor + Civil+log(Exposicao), data = dados, family=poisson) +summary(ajusteps2) +anova(ajusteps, ajusteps2, test = 'Chisq') +``` + + +### Regressão usando a distribuição binomial negativa. + +```{r} +ajustenb2 <- glm.nb(Sinistros ~ Tipo+Sexo+Idade+I(Idade**2)+Valor+Civil+log(Exposicao),data= dados) +summary(ajustenb2) + +### Verificando possibilidade de excluir estado civil e tipo de veículo do modelo. +ajustenb3 <- update(ajustenb2, ~.-(Civil+Tipo)) +anova(ajustenb2, ajustenb3) + +## Estimação do parâmetro de dispersão +X2 <- sum(resid(ajustenb3,type='pearson')**2) +phichap <- X2/ajustenb3$df.residual +phichap +``` + + +### Diagnóstico do ajuste (gráficos). + +```{r} +##### Diagnóstico do modelo - gráficos. +par(mfrow=c(2,2)) +plot(ajustenb3) +``` + + +```{r, echo = FALSE} +dadosnb3 <- dados[,c('Sexo','Idade','Valor','Exposicao','Sinistros')] + +par(mfrow=c(1,1)) +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) +} +``` + + +```{r} +par(mfrow=c(1,1)) +envelope(ajustenb3) +``` + + + +### Explorando os efeitos das covariáveis. Estimativas pontuais e ICs (95%) + +```{r} +intervalos <- confint(ajustenb3) +estimat <- cbind(ajustenb3$coefficients, intervalos) +colnames(estimat)[1] <- 'Estimativa pontual' + +### Quadro de estimativas +round(estimat, 5) +``` + + +### Gráficos de efeitos + +```{r} + +dados$lexpo <- log(dados$Exposicao) +ajustenb3 <- glm.nb(Sinistros ~ Sexo+Idade+I(Idade**2)+Valor+lexpo,data= dados) + +efeitos <- allEffects(ajustenb3, 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)) + +``` + + + + +```{r, echo = FALSE} + +## 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(ajusteps,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)); freqbn + + +### Binomial negativa com covariaveis +pred2 <- predict(ajustenb3,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 = ajustenb3$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') +``` + + +### Frequências amostrais e frequências ajustadas pelas distribuições Poisson e Binomial Negativa (sem e com inclusão de covariáveis) + +```{r, results = 'markup'} +kable(matfreq, format = "markdown", caption = "Frequências amostrais e freuências ajustadas para o número de sinistros") +``` + + + + + + + + + + +