Skip to content
Snippets Groups Projects
Commit 803b3ff3 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Adiciona dados e uma vinheta.

parent 05180976
Branches
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
examples/Ultimate_Gamma_Count.Rmd 0 → 100644
+ 1276
− 0
View file @ 803b3ff3
---
title: "Gamma-Count na análise de dados de contagem"
author: "Walmes Zeviani"
output:
rmarkdown::html_vignette:
fig_width: 6
fig_height: 6
toc: true
toc_dep: 3
vignette: >
%\VignetteIndexEntry{Análise dos dados em Pimentel-Gomes (2009)}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
<style type="text/css">
body, td, caption {
font-family: "Palatino Linotype", "Book Antiqua", Palatino, serif;
background-color: white;
font-size: 16px;
}
tt, code, pre {
font-family: "Inconsolata", "Andale Mono", monospace;
}
code {
font-size: 16px;
}
pre code {
font-size: 14px;
}
pre:not([class]) code {
background-color: #92BFB1;
}
pre, code {
background-color: #62BFB1;
border-radius: 3px;
color: #333;
}
/* R output */
pre:not([class]) code {
background-color: #D4D4D4;
}
pre:not([class]), code {
background-color: #D4D4D4;
}
/* R input */
pre, code {
border-radius: 3px;
background-color: #EDEDED;
color: #333;
}
img {
max-width: 100% !important;
display: block;
margin: auto;
}
.MathJax {
font-size: 80% !important;
}
</style>
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(
dev.args=list(family="Palatino"))
options(width=68)
```
****
## Gamma-Count
```{r, message=FALSE, error=FALSE, results="hide"}
##----------------------------------------------------------------------
## Carrega pacotes.
library(lattice)
library(latticeExtra)
library(plyr)
library(bbmle)
library(corrplot)
library(car)
library(multcomp)
library(doBy)
```
```{r, eval=FALSE}
##----------------------------------------------------------------------
## Função de probabilidade da Gamma-Count.
dgcount <- function(y, alpha, beta, offset=1){
pgamma(offset, alpha*y, alpha*beta)-
pgamma(offset, alpha*y+alpha, alpha*beta)
}
y <- 0:30
pP <- dpois(y, lambda=5)
pGC <- dgcount(y, alpha=1, beta=5)
plot(pP~y, type="h")
points(pGC~c(y+0.2), type="h", col=2)
```
```{r, eval=FALSE}
##----------------------------------------------------------------------
## Explorando a distribuição de probabilidades.
library(rpanel)
draw <- function(input){
with(input,
{
ebeta <- exp(beta)
ealpha <- exp(alpha)
pP <- dpois(y, lambda=ebeta*off)
pGC <- dgcount(y, alpha=ealpha, beta=ebeta, offset=off)
plot(pP~y, type="h", ylim=c(0, max(c(pGC, pP))))
points(pGC~c(y+0.2), type="h", col=2)
m <- sum(y*pGC)
abline(v=ebeta*off, lty=2)
abline(v=m, col=2, lty=2)
})
return(input)
}
panel <- rp.control(y=0:100)
rp.slider(panel, variable=alpha, from=-5, to=5, initval=0,
showvalue=TRUE, resolution=0.1, action=draw)
rp.slider(panel, variable=beta, from=-5, to=8, initval=3,
showvalue=TRUE, resolution=0.1, action=draw)
rp.doublebutton(panel, variable=off, range=c(0, 3), initval=1,
showvalue=TRUE, step=0.1, action=draw)
```
****
## Reprodução de nematóides em cultivo de soja
Nesse conjunto de dados, a contagem tem um *offset* que precisa ser
considerado.
### Análise exploratória
```{r}
##----------------------------------------------------------------------
## Carrega os dados.
nematodeSoybean1 <- read.table(
"nematodeSoybean1.csv",
header=TRUE, sep=";")
nematodeSoybean1$count <- with(nematodeSoybean1, count1+count2)
str(nematodeSoybean1)
##----------------------------------------------------------------------
xyplot(count1+count2~log2(inipop)|factor(das),
data=nematodeSoybean1,
type=c("p", "smooth"), jitter.x=TRUE,
ylab=expression(Nematodes~by~unit~volume),
xlab=expression(Initial~nematode~density))
xyplot(count~log2(inipop)|factor(das), data=nematodeSoybean1,
type=c("p", "smooth"), jitter.x=TRUE,
ylab=expression(Nematodes~by~unit~volume),
xlab=expression(Initial~nematode~density))
xyplot(count/(2*volume)~log2(inipop)|factor(das), data=nematodeSoybean1,
type=c("p", "smooth"), jitter.x=TRUE,
ylab=expression(Nematodes~by~unit~volume),
xlab=expression(Initial~nematode~density))
##----------------------------------------------------------------------
## Relação média-variância.
da <- subset(nematodeSoybean1, inipop>0)
mv <- ddply(.data=da,
.variables=.(inipop, das),
.fun=summarise,
m=mean(count, na.rm=TRUE),
v=var(count, na.rm=TRUE))
xyplot(v~m, data=mv, scales="free", type=c("p", "r"))+
layer(panel.abline(a=0, b=1, lty=2))
```
### Inclusão de offset no modelo Gamma-Count
```{r}
##----------------------------------------------------------------------
## Simulando um conjunto de dados para ver como acomodar o offset na
## Gamma-Count.
y <- rpois(100, lambda=5) ## Contagem.
o <- rep(2, length(y)) ## Offset.
n0 <- glm(y~1, family=poisson)
## predict(n0, type="response"); exp(coef(n0))
## logLik(n0)
n1 <- glm(y~offset(log(o))+1, family=poisson)
## predict(n1, type="response"); exp(coef(n1))*o
## logLik(n1)
ll <- function(theta, y, X, offset){
## theta: parameter vector;
## theta[1]: dispersion parameter;
## theta[-1]: location parameters;
## y: response vector (counts);
## X: model matrix;
## offset: upper limit of the intergral;
eXb <- exp(X%*%theta[-1])
alpha <- exp(theta[1])
alpha.eXb <- alpha*eXb
alpha.y <- alpha*y
## returns the log-likelihood.
-sum(log(
pgamma(offset, alpha.y, alpha.eXb)-
pgamma(offset, alpha.y+alpha, alpha.eXb)
))
}
start <- c(alpha=1, lambda=coef(n1))
parnames(ll) <- names(start)
L <- list(y=y, offset=o, X=cbind(rep(1, length(y))))
n2 <- mle2(ll, start=start, fixed=list(alpha=0),
vecpar=TRUE, data=L)
n3 <- mle2(ll, start=start,
vecpar=TRUE, data=L)
## Log-verossimilhaça.
cbind(logLik(n0), logLik(n1), logLik(n2), logLik(n3))
## Estimativas dos parâmetros.
cbind(coef(n0), coef(n1), coef(n2)[-1], coef(n3)[-1])
```
### Ajuste dos modelos
```{r}
##----------------------------------------------------------------------
## Ajuste do Poisson e Quasi-Poisson.
da <- transform(da, IP=factor(inipop), DAS=factor(das))
## Modelo Poisson para ter referência.
m0 <- glm(count~offset(log(volume))+IP*DAS, data=da, family=poisson)
summary(m0)
anova(m0, test="Chisq")
## Modelo quasi Poisson para usar os valores iniciais.
## m1 <- glm(count/volume~IP*DAS, data=da, family=quasipoisson)
m1 <- glm(count~offset(log(volume))+IP*DAS,
data=da, family=quasipoisson)
summary(m1)
anova(m1, test="F")
par(mfrow=c(2,2)); plot(m1); layout(1)
cbind(coef(m0), coef(m1))
##----------------------------------------------------------------------
## Ajuste do Gamma-Count.
start <- c("alpha"=log(1), coef(m0))
parnames(ll) <- names(start)
## y <- with(da, count/vol)
y <- with(da, count/volume)
X <- model.matrix(~IP*DAS, data=da)
L <- with(da,
list(y=count, offset=volume,
X=model.matrix(~IP*DAS)))
m2 <- mle2(ll, start=start, fixed=list(alpha=0),
data=L, vecpar=TRUE)
m3 <- mle2(ll, start=start,
data=L, vecpar=TRUE)
```
### Comparação dos ajustes
```{r}
##----------------------------------------------------------------------
## Log-verossimilhaça.
cbind(logLik(m0), logLik(m1), logLik(m2), logLik(m3))
## Estimativas dos parâmetros.
cbind(coef(m0), coef(m1), coef(m2)[-1], coef(m3)[-1])
## Estimativa do parâmetro de dispersão.
exp(coef(m3)[1])
## Perfil de log-verossmilhança.
plot(profile(m3, which=1))
## Teste para o parâmetro de dispersão.
anova(m3, m2)
## Teste para interação (Wald).
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a==3
L <- t(replicate(sum(ai), rbind(coef(m3)*0), simplify="matrix"))
L[,ai] <- diag(sum(ai))
## Teste de Wald explicito.
## t(L%*%coef(m3))%*%
## solve(L%*%vcov(m3)%*%t(L))%*%
## (L%*%coef(m3))
crossprod(L%*%coef(m3), solve(L%*%vcov(m3)%*%t(L), L%*%coef(m3)))
## Teste de Wald para interação (poderia ser LRT, claro).
linearHypothesis(model=m0, ## É necessário um objeto glm.
hypothesis.matrix=L,
vcov.=vcov(m3),
coef.=coef(m3))
## Teste pelo modelo quasi Poisson.
anova(m1, test="Chisq")
anova(m1, test="F")
##----------------------------------------------------------------------
## Ortogonalidade entre dispersão e locação.
V <- cov2cor(vcov(m3))
corrplot.mixed(V, upper="ellipse", col="gray50")
dev.off()
V[1,-1]
sum(abs(V[1,-1]))
```
### Predição
```{r}
##----------------------------------------------------------------------
## Predição das médias com IC.
pred <- expand.grid(IP=levels(da$IP),
DAS=levels(da$DAS),
volume=100,
KEEP.OUT.ATTRS=FALSE)
pred <- list(pois=pred, quasi=pred, cg=pred)
## Quantil normal.
qn <- qnorm(0.975)*c(lwr=-1, fit=0, upr=1)
## Preditos pela Poisson.
aux <- predict(m0, newdata=pred$pois, se.fit=TRUE)
aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
pred$pois <- cbind(pred$pois, aux)
str(pred$pois)
aux <- predict(m1, newdata=pred$quasi, se.fit=TRUE)
aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
pred$quasi <- cbind(pred$quasi, aux)
alpha <- coef(m3)[1] ## Locação.
theta <- coef(m3)[-1] ## Dispersão.
## Preditor linear.
X <- model.matrix(~IP*DAS, data=pred$cg)
pred$cg$eta <- c(X%*%theta)
## Matrix de covariância completa.
V <- vcov(m3)
## Marginal em theta.
Vm <- V[-1,-1]
## Condicional.
Vc <- V[-1,-1]-V[-1,1]%*%solve(V[1,1])%*%V[1,-1]
max(abs(Vm-Vc))
## Preditos pela Gamma-count.
U <- chol(Vm)
aux <- sqrt(apply(X%*%t(U), MARGIN=1,
FUN=function(x){ sum(x^2) }))
## Para calcular a média de uma Gamma-Count pelo preditor linear.
eta2mean <- function(eta, offset, alpha, tol=1e-5){
stopifnot(length(eta)==1L)
stopifnot(length(offset)==1L)
stopifnot(length(alpha)==1L)
change <- 1
S <- 0; x <- 1; p <- 1
while (p>tol){
p <- pgamma(offset,
shape=exp(alpha)*x,
rate=exp(alpha)*exp(eta))
S <- S+p
x <- x+1
}
return(S)
}
aux <- pred$cg$eta+outer(aux, qn, FUN="*")
aux <- apply(aux, MARGIN=2,
FUN=function(col){
sapply(col, FUN=eta2mean,
offset=pred$cg$volume[1], alpha=alpha)
})
pred$cg <- cbind(pred$cg, aux)
##----------------------------------------------------------------------
pred <- ldply(pred, .id="modelo")
pred <- arrange(pred, DAS, IP, modelo)
str(pred)
source("http://git.leg.ufpr.br/leg/legTools/raw/devel/R/panel.segplot.by.R")
segplot(IP~lwr+upr|DAS, centers=fit, data=pred,
draw=FALSE, scales=list(x="free"),
panel=panel.segplot.by, groups=pred$modelo, f=0.1,
pch=pred$modelo, layout=c(1, NA), as.table=TRUE,
key=list(type="o", divide=1,
lines=list(pch=1:nlevels(pred$modelo), lty=1, col=1),
text=list(c("Poisson", "Quasi-Poisson",
"Gamma-Count"))))
```
****
## Número de vagens viáveis em soja
Nesse conjunto de dados, **o efeito de bloco pode ser aleatório**. Teria
que implementar uma versão da Gamma-Count para acomodar termos de efeito
aleatório.
### Análise exploratória
```{r}
##----------------------------------------------------------------------
## Carrega os dados.
db <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
header=TRUE, sep="\t", dec=",")
names(db)[1:2] <- c("k", "a")
db <- db[-74,c(1:3, 8, 10)]
str(db)
xyplot(nv~k|a, data=db)
xyplot(ts~k|a, data=db)
##----------------------------------------------------------------------
## Gráfico de média-variância não é possível pois não existem repetições
## com o mesmo preditor.
```
### Ajuste dos modelos
```{r}
##----------------------------------------------------------------------
## Ajuste do Poisson e Quasi-Poisson.
db <- transform(db, A=factor(a), K=factor(k))
## Modelo Poisson para ter referência.
m0 <- glm(nv~bloco+A*K, data=db, family=poisson)
summary(m0)
anova(m0, test="Chisq")
## Modelo quasi Poisson para usar os valores iniciais.
## m1 <- glm(count/volume~IP*DAS, data=da, family=quasipoisson)
m1 <- update(m0, family=quasipoisson)
summary(m1)
anova(m1, test="F")
par(mfrow=c(2,2)); plot(m1); layout(1)
cbind(coef(m0), coef(m1))
##----------------------------------------------------------------------
## Ajuste do Gamma-Count.
start <- c("alpha"=log(1), coef(m0))
parnames(ll) <- names(start)
y <- with(db, nv)
X <- model.matrix(~bloco+A*K, data=db)
L <- with(db,
list(y=nv, offset=1,
X=model.matrix(~bloco+A*K)))
m2 <- mle2(ll, start=start, fixed=list(alpha=0),
data=L, vecpar=TRUE)
m3 <- mle2(ll, start=start,
data=L, vecpar=TRUE)
```
### Comparação dos ajustes
```{r}
##----------------------------------------------------------------------
## Log-verossimilhaça.
cbind(logLik(m0), logLik(m1), logLik(m2), logLik(m3))
## Estimativas dos parâmetros.
cbind(coef(m0), coef(m1), coef(m2)[-1], coef(m3)[-1])
## Estimativa do parâmetro de dispersão.
exp(coef(m3)[1])
## Perfil de log-verossmilhança.
plot(profile(m3, which=1))
## Teste para o parâmetro de dispersão.
anova(m3, m2)
## Teste para interação (Wald).
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a==4
L <- t(replicate(sum(ai), rbind(coef(m3)*0), simplify="matrix"))
L[,ai] <- diag(sum(ai))
## Teste de Wald explicito.
## t(L%*%coef(m3))%*%
## solve(L%*%vcov(m3)%*%t(L))%*%
## (L%*%coef(m3))
crossprod(L%*%coef(m3), solve(L%*%vcov(m3)%*%t(L), L%*%coef(m3)))
## Teste de Wald para interação (poderia ser LRT, claro).
linearHypothesis(model=m0, ## É necessário um objeto glm.
hypothesis.matrix=L,
vcov.=vcov(m3),
coef.=coef(m3))
## Teste pelo modelo quasi Poisson.
anova(m1, test="Chisq")
anova(m1, test="F")
##----------------------------------------------------------------------
## Ortogonalidade entre dispersão e locação.
V <- vcov(m3)
corrplot.mixed(V, upper="ellipse", col="gray50")
dev.off()
V[1,-1]
sum(abs(V[1,-1]))
```
### Predição
```{r}
##----------------------------------------------------------------------
## Predição das médias com IC.
pred <- expand.grid(A=levels(db$A),
K=levels(db$K),
bloco=factor("I", levels=levels(db$bloco)),
KEEP.OUT.ATTRS=FALSE)
pred <- list(pois=pred, quasi=pred, cg=pred)
## Quantil normal.
qn <- qnorm(0.975)*c(lwr=-1, fit=0, upr=1)
## Preditor linear que considera o efeito médio dos blocos.
X <- LSmatrix(m0, effect=c("A", "K"))
## Preditos pela Poisson.
## aux <- predict(m0, newdata=pred$pois, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
## aux <- predict(m1, newdata=pred$quasi, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$quasi <- cbind(pred$quasi, aux)
aux <- confint(glht(m1, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
alpha <- coef(m3)[1] ## Locação.
theta <- coef(m3)[-1] ## Dispersão.
## Preditor linear.
pred$cg$eta <- c(X%*%theta)
## Matrix de covariância completa.
V <- vcov(m3)
## Marginal em theta.
Vm <- V[-1,-1]
## Condicional.
Vc <- V[-1,-1]-V[-1,1]%*%solve(V[1,1])%*%V[1,-1]
max(abs(Vm-Vc))
## Preditos pela Gamma-count.
U <- chol(Vm)
aux <- sqrt(apply(X%*%t(U), MARGIN=1,
FUN=function(x){ sum(x^2) }))
aux <- pred$cg$eta+outer(aux, qn, FUN="*")
aux <- apply(aux, MARGIN=2,
FUN=function(col){
sapply(col, FUN=eta2mean,
offset=1, alpha=alpha)
})
pred$cg <- cbind(pred$cg, aux)
##----------------------------------------------------------------------
pred <- ldply(pred, .id="modelo")
pred <- arrange(pred, A, K, modelo)
str(pred)
xyplot(nv~K|A, data=db, layout=c(NA, 1), as.table=TRUE,
col=1, cex=0.7, pch=19)+
as.layer(
segplot(K~lwr+upr|A, centers=fit, data=pred,
draw=FALSE, horizontal=FALSE, ## scales=list(x="free"),
panel=panel.segplot.by, groups=pred$modelo, f=0.15,
pch=pred$modelo, layout=c(NA, 1), as.table=TRUE,
key=list(type="o", divide=1,
lines=list(pch=1:nlevels(pred$modelo),
lty=1, col=1),
text=list(c("Poisson", "Quasi-Poisson",
"Gamma-Count"))))
)
```
## Ocorrência de mosca branca
<http://www.leg.ufpr.br/~walmes/data/ninfas.txt>
### Análise exploratória
```{r}
##----------------------------------------------------------------------
## Análise exploratória.
dc <- read.table("moscaBranca.txt", header=TRUE, sep="\t")
dc <- transform(dc,
data=as.Date(dc$data),
bloc=factor(bloc),
tot=inf+med+sup,
off=100)
dc$aval <- with(dc, c(data-min(data)))
dc$Aval <- factor(dc$aval)
names(dc)[2:3] <- c("Vari","Bloc")
dc <- droplevels(subset(dc, grepl("BRS", x=Vari)))
str(dc)
## xyplot(sup+med+inf~data, groups=Vari, outer=TRUE, data=dc)
xyplot(tot~data|Vari, data=dc)
```
### Ajuste dos modelos
```{r}
##----------------------------------------------------------------------
## Ajuste do Poisson e Quasi-Poisson.
## IMPORTANT: Foi usado um offset de 100 por problemas de
## over/underflow. O problema não é com a Poisson mas com a
## Gamma-Count. Todas as análises consideram um offset de 100, portanto,
## permanecem comparáveis e igualmente interpretáveis.
## Modelo Poisson para ter referência.
m0 <- glm(tot~offset(log(off))+Bloc+Vari*Aval,
data=dc, family=poisson)
summary(m0)
anova(m0, test="Chisq")
## Modelo quasi Poisson para usar os valores iniciais.
## m1 <- glm(count/volume~IP*DAS, data=da, family=quasipoisson)
m1 <- update(m0, family=quasipoisson)
summary(m1)
anova(m1, test="F")
par(mfrow=c(2,2)); plot(m1); layout(1)
anova(m1, test="F")
##----------------------------------------------------------------------
## Ajuste do Gamma-Count.
start <- c("alpha"=log(0.05), coef(m0))
parnames(ll) <- names(start)
y <- with(dc, tot)
X <- model.matrix(~Bloc+Vari*Aval, data=dc)
L <- with(dc, list(y=tot, offset=100, X=X))
m3 <- mle2(ll, start=start,
data=L, vecpar=TRUE)
## summary(m3)
```
### Comparação dos ajustes
```{r}
##----------------------------------------------------------------------
## Log-verossimilhaça.
cbind(logLik(m0), logLik(m1), logLik(m3))
## Estimativas dos parâmetros.
cbind(coef(m0), coef(m1), coef(m3)[-1])
## Estimativa do parâmetro de dispersão.
exp(coef(m3)[1])
## Perfil de log-verossmilhança.
plot(profile(m3, which=1))
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a==4
L <- t(replicate(sum(ai), rbind(coef(m3)*0), simplify="matrix"))
L[,ai] <- diag(sum(ai))
## Teste de Wald explicito.
## t(L%*%coef(m3))%*%
## solve(L%*%vcov(m3)%*%t(L))%*%
## (L%*%coef(m3))
crossprod(L%*%coef(m3), solve(L%*%vcov(m3)%*%t(L), L%*%coef(m3)))
## Teste de Wald para interação (poderia ser LRT, claro).
linearHypothesis(model=m0, ## É necessário um objeto glm.
hypothesis.matrix=L,
vcov.=vcov(m3),
coef.=coef(m3))
## Teste pelo modelo quasi Poisson.
anova(m1, test="Chisq")
anova(m1, test="F")
##----------------------------------------------------------------------
## Ortogonalidade entre dispersão e locação.
V <- cov2cor(vcov(m3))
corrplot.mixed(V, upper="ellipse", col="gray50")
dev.off()
```
### Predição
```{r}
##----------------------------------------------------------------------
## Predição das médias com IC.
pred <- expand.grid(Vari=levels(dc$Vari),
Aval=unique(dc$Aval),
Bloc=factor("1", levels=levels(dc$Bloc)),
off=100,
KEEP.OUT.ATTRS=FALSE)
pred <- list(pois=pred, quasi=pred, cg=pred)
## Quantil normal.
qn <- qnorm(0.975)*c(lwr=-1, fit=0, upr=1)
X <- LSmatrix(lm(tot~Bloc+Vari*Aval, dc), effect=c("Vari", "Aval"))
## Preditos pela Poisson.
## aux <- predict(m0, newdata=pred$pois, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, 100*exp(aux))
## aux <- predict(m1, newdata=pred$quasi, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$quasi <- cbind(pred$quasi, aux)
aux <- confint(glht(m1, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, 100*exp(aux))
alpha <- coef(m3)[1] ## Locação.
theta <- coef(m3)[-1] ## Dispersão.
## Preditor linear.
pred$cg$eta <- c(X%*%theta)
## Matrix de covariância completa.
V <- vcov(m3)
## Marginal em theta.
Vm <- V[-1,-1]
## Condicional.
Vc <- V[-1,-1]-V[-1,1]%*%solve(V[1,1])%*%V[1,-1]
max(abs(Vm-Vc))
## Preditos pela Gamma-count.
U <- chol(Vm)
aux <- sqrt(apply(X%*%t(U), MARGIN=1,
FUN=function(x){ sum(x^2) }))
aux <- pred$cg$eta+outer(aux, qn, FUN="*")
aux <- apply(aux, MARGIN=2,
FUN=function(col){
sapply(col, FUN=eta2mean,
offset=100, alpha=alpha)
})
pred$cg <- cbind(pred$cg, aux)
##----------------------------------------------------------------------
pred <- ldply(pred, .id="modelo")
pred <- arrange(pred, Aval, Vari, modelo)
str(pred)
## source("http://git.leg.ufpr.br/leg/legTools/raw/devel/R/panel.segplot.by.R")
xyplot(tot~Aval|Vari, data=dc, layout=c(NA, 1), as.table=TRUE,
col=1, cex=0.7, pch=19)+
as.layer(
segplot(Aval~lwr+upr|Vari, centers=fit, data=pred,
horizontal=FALSE,
draw=FALSE,
panel=panel.segplot.by, groups=pred$modelo, f=0.15,
pch=pred$modelo, layout=c(NA, 1), as.table=TRUE,
key=list(type="o", divide=1,
lines=list(pch=1:nlevels(pred$modelo),
lty=1, col=1),
text=list(c("Poisson", "Quasi-Poisson",
"Gamma-Count"))))
)
```
****
## Produção de ovos
### Análise exploratória
```{r}
##----------------------------------------------------------------------
dd <- read.table("ovos.txt", header=TRUE, sep="\t")
str(dd)
ftable(xtabs(~periodo+luz+box, data=dd))
## Cada box tem 6 gaiolas. Cada gaiola tem 10 aves.
## Modelar a produção por ave dia usando o total de por gaiola em duas
## semanas. O offset é 10*14=140.
useOuterStrips(
xyplot(ovos~dia|luz+periodo, data=dd, jitter.x=TRUE,
groups=gaiola,
type=c("p", "smooth")))
useOuterStrips(
xyplot(massa~dia|luz+periodo, data=dd, jitter.x=TRUE,
groups=gaiola,
type=c("p", "smooth")))
## Agregar para o total quinzenal.
dd <- aggregate(ovos~periodo+box+luz+gaiola, data=dd, FUN=sum)
str(dd)
xyplot(ovos~periodo, groups=luz, data=dd, jitter.x=TRUE,
## groups=gaiola,
type=c("p", "smooth"))
```
### Ajuste dos modelos
```{r}
##----------------------------------------------------------------------
## Ajuste do Poisson e Quasi-Poisson.
## IMPORTANT: Foi usado um offset de 10*14=140.
dd <- transform(dd,
off=10*14, ## 10 aves por gaiola, soma de 14 dias.
Per=factor(periodo),
Box=factor(box))
names(dd)[3] <- "Luz"
## WARNING: Esta sendo negligenciada a variância de gaiolas. Isso requer
## um modelo de efeitos aleatórios que ainda não dispomos.
## Modelo Poisson para ter referência.
m0 <- glm(ovos~offset(log(off))+Per+Box+Luz,
data=dd, family=poisson)
summary(m0)
anova(m0, test="Chisq")
## Modelo quasi Poisson para usar os valores iniciais.
## m1 <- glm(count/volume~IP*DAS, data=da, family=quasipoisson)
m1 <- update(m0, family=quasipoisson)
summary(m1)
anova(m1, test="F")
par(mfrow=c(2,2)); plot(m1); layout(1)
##----------------------------------------------------------------------
## Ajuste do Gamma-Count.
start <- c("alpha"=log(1), coef(m0))
parnames(ll) <- names(start)
y <- with(dd, ovos)
X <- model.matrix(~Per+Box+Luz, data=dd)
L <- list(y=y, offset=140, X=X)
m3 <- mle2(ll, start=start,
data=L, vecpar=TRUE)
## summary(m3)
```
### Comparação dos ajustes
```{r}
##----------------------------------------------------------------------
## Log-verossimilhaça.
cbind(logLik(m0), logLik(m1), logLik(m3))
## Estimativas dos parâmetros.
cbind(coef(m0), coef(m1), coef(m3)[-1])
## Estimativa do parâmetro de dispersão.
exp(coef(m3)[1])
## Perfil de log-verossmilhança.
plot(profile(m3, which=1))
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a==3
L <- t(replicate(sum(ai), rbind(coef(m3)*0), simplify="matrix"))
L[,ai] <- diag(sum(ai))
## Teste de Wald explicito.
## t(L%*%coef(m3))%*%
## solve(L%*%vcov(m3)%*%t(L))%*%
## (L%*%coef(m3))
crossprod(L%*%coef(m3), solve(L%*%vcov(m3)%*%t(L), L%*%coef(m3)))
## Teste de Wald para interação (poderia ser LRT, claro).
linearHypothesis(model=m0, ## É necessário um objeto glm.
hypothesis.matrix=L,
vcov.=vcov(m3),
coef.=coef(m3))
## Teste pelo modelo quasi Poisson.
anova(m1, test="Chisq")
anova(m1, test="F")
##----------------------------------------------------------------------
## Ortogonalidade entre dispersão e locação.
V <- cov2cor(vcov(m3))
corrplot.mixed(V, upper="ellipse", col="gray50")
dev.off()
```
### Predição
```{r}
##----------------------------------------------------------------------
## Predição das médias com IC.
pred <- expand.grid(Per=factor(levels(dd$Per)[1],
levels=levels(dd$Per)),
Box=factor(levels(dd$Box)[1],
levels=levels(dd$Box)),
Luz=levels(dd$Luz),
off=10*14,
KEEP.OUT.ATTRS=FALSE)
pred <- list(pois=pred, quasi=pred, cg=pred)
## Quantil normal.
qn <- qnorm(0.975)*c(lwr=-1, fit=0, upr=1)
X <- LSmatrix(lm(ovos~Per+Box+Luz, dd), effect=c("Luz"))
## Preditos pela Poisson.
## aux <- predict(m0, newdata=pred$pois, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, 140*exp(aux))
## aux <- predict(m1, newdata=pred$quasi, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$quasi <- cbind(pred$quasi, aux)
aux <- confint(glht(m1, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, 140*exp(aux))
alpha <- coef(m3)[1] ## Locação.
theta <- coef(m3)[-1] ## Dispersão.
## Preditor linear.
pred$cg$eta <- c(X%*%theta)
## Matrix de covariância completa.
V <- vcov(m3)
## Marginal em theta.
Vm <- V[-1,-1]
## Condicional.
Vc <- V[-1,-1]-V[-1,1]%*%solve(V[1,1])%*%V[1,-1]
max(abs(Vm-Vc))
## Preditos pela Gamma-count.
U <- chol(Vm)
aux <- sqrt(apply(X%*%t(U), MARGIN=1,
FUN=function(x){ sum(x^2) }))
aux <- pred$cg$eta+outer(aux, qn, FUN="*")
aux <- apply(aux, MARGIN=2,
FUN=function(col){
sapply(col, FUN=eta2mean,
offset=140, alpha=alpha)
})
pred$cg <- cbind(pred$cg, aux)
##----------------------------------------------------------------------
pred <- ldply(pred, .id="modelo")
pred <- arrange(pred, Luz, modelo)
str(pred)
## source("http://git.leg.ufpr.br/leg/legTools/raw/devel/R/panel.segplot.by.R")
segplot(Luz~lwr+upr, centers=fit, data=pred,
horizontal=FALSE,
draw=FALSE,
panel=panel.segplot.by, groups=pred$modelo, f=0.15,
pch=pred$modelo, layout=c(NA, 1), as.table=TRUE,
key=list(type="o", divide=1,
lines=list(pch=1:nlevels(pred$modelo),
lty=1, col=1),
text=list(c("Poisson", "Quasi-Poisson",
"Gamma-Count"))))
```
****
## Capulhos no algodão em função da pressão da mosca branca
<http://www.leg.ufpr.br/~walmes/data/mosca_algodao_prod.txt>
### Análise exploratória
```{r}
##----------------------------------------------------------------------
da <- read.table(
"http://www.leg.ufpr.br/~walmes/data/mosca_algodao_prod.txt",
header=TRUE, sep="\t")
## da <- aggregate(ncapu~vaso+dexp, data=da, FUN=sum)
str(da)
xyplot(ncapu~dexp, data=da, jitter.x=TRUE,
type=c("p", "smooth"))
```
### Ajuste dos modelos
```{r}
##----------------------------------------------------------------------
## Ajuste do Poisson e Quasi-Poisson.
## da <- transform(da, Dexp=factor(dexp))
da <- transform(da, dexp=dexp-mean(range(dexp)))
## WARNING: Esta sendo negligenciada a variância entre vasos. Isso
## requer um modelo de efeitos aleatórios que ainda não dispomos.
## Modelo Poisson para ter referência.
m0 <- glm(ncapu~dexp+I(dexp^2), data=da, family=poisson)
summary(m0)
anova(m0, test="Chisq")
## Modelo quasi Poisson para usar os valores iniciais.
m1 <- update(m0, family=quasipoisson)
summary(m1)
anova(m1, test="F")
par(mfrow=c(2,2)); plot(m1); layout(1)
##----------------------------------------------------------------------
## Ajuste do Gamma-Count.
start <- c("alpha"=log(2), coef(m0))
parnames(ll) <- names(start)
y <- with(da, ncapu)
X <- model.matrix(m0)
L <- list(y=y, offset=1, X=X)
m3 <- mle2(ll, start=start,
data=L, vecpar=TRUE)
summary(m3)
```
### Comparação dos ajustes
```{r}
##----------------------------------------------------------------------
## Log-verossimilhaça.
cbind(logLik(m0), logLik(m1), logLik(m3))
## Estimativas dos parâmetros.
cbind(coef(m0), coef(m1), coef(m3)[-1])
## Estimativa do parâmetro de dispersão.
exp(coef(m3)[1])
## Perfil de log-verossmilhança.
plot(profile(m3, which=1))
a <- c(0, attr(model.matrix(m0), "assign"))
ai <- a!=0
L <- t(replicate(sum(ai), rbind(coef(m3)*0), simplify="matrix"))
L[,ai] <- diag(sum(ai))
## Teste de Wald explicito.
## t(L%*%coef(m3))%*%
## solve(L%*%vcov(m3)%*%t(L))%*%
## (L%*%coef(m3))
crossprod(L%*%coef(m3), solve(L%*%vcov(m3)%*%t(L), L%*%coef(m3)))
## Teste de Wald para interação (poderia ser LRT, claro).
linearHypothesis(model=m0, ## É necessário um objeto glm.
hypothesis.matrix=L,
vcov.=vcov(m3),
coef.=coef(m3))
## Teste pelo modelo quasi Poisson.
anova(m1, test="Chisq")
anova(m1, test="F")
##----------------------------------------------------------------------
## Ortogonalidade entre dispersão e locação.
V <- cov2cor(vcov(m3))
corrplot.mixed(V, upper="ellipse", col="gray50")
dev.off()
```
### Predição
```{r}
##----------------------------------------------------------------------
## Predição das médias com IC.
pred <- data.frame(dexp=with(da, seq(min(dexp), max(dexp), l=30)))
pred <- list(pois=pred, quasi=pred, cg=pred)
## Quantil normal.
qn <- qnorm(0.975)*c(lwr=-1, fit=0, upr=1)
X <- model.matrix(formula(m0)[-2], data=pred$cg)
## Preditos pela Poisson.
## aux <- predict(m0, newdata=pred$pois, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$pois <- cbind(pred$pois, exp(aux))
## aux <- predict(m1, newdata=pred$quasi, se.fit=TRUE)
## aux <- exp(aux$fit+outer(aux$se.fit, qn, FUN="*"))
## pred$quasi <- cbind(pred$quasi, aux)
aux <- confint(glht(m1, linfct=X), calpha=univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred$quasi <- cbind(pred$quasi, exp(aux))
alpha <- coef(m3)[1] ## Locação.
theta <- coef(m3)[-1] ## Dispersão.
## Preditor linear.
pred$cg$eta <- c(X%*%theta)
## Matrix de covariância completa.
V <- vcov(m3)
## Marginal em theta.
Vm <- V[-1,-1]
## Condicional.
Vc <- V[-1,-1]-V[-1,1]%*%solve(V[1,1])%*%V[1,-1]
max(abs(Vm-Vc))
## Preditos pela Gamma-count.
U <- chol(Vm)
aux <- sqrt(apply(X%*%t(U), MARGIN=1,
FUN=function(x){ sum(x^2) }))
aux <- pred$cg$eta+outer(aux, qn, FUN="*")
aux <- apply(aux, MARGIN=2,
FUN=function(col){
sapply(col, FUN=eta2mean,
offset=1, alpha=alpha)
})
pred$cg <- cbind(pred$cg, aux)
##----------------------------------------------------------------------
pred <- ldply(pred, .id="modelo")
pred <- arrange(pred, dexp)
str(pred)
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/prepanel.cbH.R")
source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/panel.cbH.R")
xyplot(ncapu~dexp, data=da, jitter.x=TRUE, pch=19, col=1)+
as.layer(
under=TRUE,
xyplot(fit~dexp, data=pred, groups=modelo, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands",
## fill="gray90",
alpha=0.5,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)
)
```
****
## Mais dados para se considerar
* http://www.leg.ufpr.br/~walmes/data/banana_culttecido.txt
* http://www.leg.ufpr.br/~walmes/data/frango_comportamento.txt
* http://www.leg.ufpr.br/~walmes/data/mosca_algodao_aval.txt
examples/Ultimate_Gamma_Count.html 0 → 100644
+ 2378
− 0
View file @ 803b3ff3
This diff is collapsed.
examples/moscaBranca.txt 0 → 100644
+ 241
− 0
View file @ 803b3ff3
data vari bloc sup med inf
2009-12-11 BRS 245 RR 1 168 71 39
2009-12-11 BRS 245 RR 2 30 46 22
2009-12-11 BRS 245 RR 3 28 36 7
2009-12-11 BRS 245 RR 4 6 1 31
2009-12-11 BRS 243 RR 1 31 22 28
2009-12-11 BRS 243 RR 2 4 6 10
2009-12-11 BRS 243 RR 3 3 12 27
2009-12-11 BRS 243 RR 4 48 76 13
2009-12-11 BRS 246 RR 1 100 50 34
2009-12-11 BRS 246 RR 2 54 15 7
2009-12-11 BRS 246 RR 3 1 50 12
2009-12-11 BRS 246 RR 4 9 2 13
2009-12-11 BRS 239 1 40 32 44
2009-12-11 BRS 239 2 76 38 16
2009-12-11 BRS 239 3 52 43 25
2009-12-11 BRS 239 4 23 41 16
2009-12-11 EMBRAPA 48 1 84 61 56
2009-12-11 EMBRAPA 48 2 129 189 133
2009-12-11 EMBRAPA 48 3 84 35 12
2009-12-11 EMBRAPA 48 4 34 17 30
2009-12-11 CD 214 RR 1 17 137 3
2009-12-11 CD 214 RR 2 107 95 31
2009-12-11 CD 214 RR 3 130 120 26
2009-12-11 CD 214 RR 4 0 34 25
2009-12-11 CD 202 1 16 21 12
2009-12-11 CD 202 2 260 34 44
2009-12-11 CD 202 3 91 67 19
2009-12-11 CD 202 4 5 2 7
2009-12-11 M 7908 RR 1 27 37 29
2009-12-11 M 7908 RR 2 103 72 32
2009-12-11 M 7908 RR 3 50 60 32
2009-12-11 M 7908 RR 4 53 26 41
2009-12-11 VMAX RR 1 24 12 11
2009-12-11 VMAX RR 2 67 6 15
2009-12-11 VMAX RR 3 0 0 0
2009-12-11 VMAX RR 4 12 26 36
2009-12-11 CD 219 RR 1 0 22 16
2009-12-11 CD 219 RR 2 391 208 88
2009-12-11 CD 219 RR 3 53 81 23
2009-12-11 CD 219 RR 4 38 39 93
2009-12-19 BRS 245 RR 1 147 54 37
2009-12-19 BRS 245 RR 2 98 65 17
2009-12-19 BRS 245 RR 3 0 53 4
2009-12-19 BRS 245 RR 4 16 11 6
2009-12-19 BRS 243 RR 1 79 96 38
2009-12-19 BRS 243 RR 2 4 20 19
2009-12-19 BRS 243 RR 3 19 2 12
2009-12-19 BRS 243 RR 4 8 3 5
2009-12-19 BRS 246 RR 1 77 42 32
2009-12-19 BRS 246 RR 2 72 110 23
2009-12-19 BRS 246 RR 3 64 17 14
2009-12-19 BRS 246 RR 4 25 7 13
2009-12-19 BRS 239 1 53 14 7
2009-12-19 BRS 239 2 144 38 12
2009-12-19 BRS 239 3 136 48 42
2009-12-19 BRS 239 4 49 23 30
2009-12-19 EMBRAPA 48 1 69 38 19
2009-12-19 EMBRAPA 48 2 172 158 83
2009-12-19 EMBRAPA 48 3 150 104 32
2009-12-19 EMBRAPA 48 4 22 12 8
2009-12-19 CD 214 RR 1 36 25 10
2009-12-19 CD 214 RR 2 59 101 27
2009-12-19 CD 214 RR 3 376 115 82
2009-12-19 CD 214 RR 4 363 170 62
2009-12-19 CD 202 1 14 21 19
2009-12-19 CD 202 2 175 160 51
2009-12-19 CD 202 3 77 51 23
2009-12-19 CD 202 4 38 34 6
2009-12-19 M 7908 RR 1 52 15 12
2009-12-19 M 7908 RR 2 138 261 139
2009-12-19 M 7908 RR 3 87 48 19
2009-12-19 M 7908 RR 4 21 10 16
2009-12-19 VMAX RR 1 53 17 9
2009-12-19 VMAX RR 2 129 65 68
2009-12-19 VMAX RR 3 11 18 7
2009-12-19 VMAX RR 4 12 19 14
2009-12-19 CD 219 RR 1 11 40 18
2009-12-19 CD 219 RR 2 85 1157 278
2009-12-19 CD 219 RR 3 95 75 24
2009-12-19 CD 219 RR 4 24 34 135
2009-12-24 BRS 245 RR 1 138 54 39
2009-12-24 BRS 245 RR 2 47 56 14
2009-12-24 BRS 245 RR 3 32 14 10
2009-12-24 BRS 245 RR 4 24 21 5
2009-12-24 BRS 243 RR 1 85 38 11
2009-12-24 BRS 243 RR 2 21 7 19
2009-12-24 BRS 243 RR 3 8 3 14
2009-12-24 BRS 243 RR 4 10 30 4
2009-12-24 BRS 246 RR 1 99 85 40
2009-12-24 BRS 246 RR 2 84 56 23
2009-12-24 BRS 246 RR 3 17 23 13
2009-12-24 BRS 246 RR 4 51 19 23
2009-12-24 BRS 239 1 156 42 22
2009-12-24 BRS 239 2 60 66 10
2009-12-24 BRS 239 3 144 57 11
2009-12-24 BRS 239 4 75 23 15
2009-12-24 EMBRAPA 48 1 135 44 43
2009-12-24 EMBRAPA 48 2 87 169 62
2009-12-24 EMBRAPA 48 3 268 107 111
2009-12-24 EMBRAPA 48 4 20 14 3
2009-12-24 CD 214 RR 1 85 47 25
2009-12-24 CD 214 RR 2 76 108 36
2009-12-24 CD 214 RR 3 265 155 146
2009-12-24 CD 214 RR 4 464 223 170
2009-12-24 CD 202 1 25 4 9
2009-12-24 CD 202 2 182 100 23
2009-12-24 CD 202 3 77 53 22
2009-12-24 CD 202 4 64 65 10
2009-12-24 M 7908 RR 1 39 13 16
2009-12-24 M 7908 RR 2 216 152 151
2009-12-24 M 7908 RR 3 62 32 22
2009-12-24 M 7908 RR 4 82 24 19
2009-12-24 VMAX RR 1 57 9 16
2009-12-24 VMAX RR 2 88 176 26
2009-12-24 VMAX RR 3 0 0 0
2009-12-24 VMAX RR 4 71 4 11
2009-12-24 CD 219 RR 1 89 99 14
2009-12-24 CD 219 RR 2 15 751 423
2009-12-24 CD 219 RR 3 131 79 13
2009-12-24 CD 219 RR 4 140 159 215
2010-01-02 BRS 245 RR 1 25 11 30
2010-01-02 BRS 245 RR 2 8 11 3
2010-01-02 BRS 245 RR 3 9 2 0
2010-01-02 BRS 245 RR 4 18 8 3
2010-01-02 BRS 243 RR 1 42 19 3
2010-01-02 BRS 243 RR 2 6 3 6
2010-01-02 BRS 243 RR 3 2 0 3
2010-01-02 BRS 243 RR 4 8 25 3
2010-01-02 BRS 246 RR 1 24 30 14
2010-01-02 BRS 246 RR 2 18 7 2
2010-01-02 BRS 246 RR 3 7 2 1
2010-01-02 BRS 246 RR 4 25 14 11
2010-01-02 BRS 239 1 38 20 55
2010-01-02 BRS 239 2 21 4 5
2010-01-02 BRS 239 3 21 23 7
2010-01-02 BRS 239 4 8 12 5
2010-01-02 EMBRAPA 48 1 10 7 13
2010-01-02 EMBRAPA 48 2 27 13 12
2010-01-02 EMBRAPA 48 3 101 23 18
2010-01-02 EMBRAPA 48 4 11 8 6
2010-01-02 CD 214 RR 1 27 16 27
2010-01-02 CD 214 RR 2 12 35 11
2010-01-02 CD 214 RR 3 62 24 54
2010-01-02 CD 214 RR 4 36 48 3
2010-01-02 CD 202 1 18 8 6
2010-01-02 CD 202 2 41 6 17
2010-01-02 CD 202 3 28 7 10
2010-01-02 CD 202 4 21 12 14
2010-01-02 M 7908 RR 1 22 11 2
2010-01-02 M 7908 RR 2 78 30 18
2010-01-02 M 7908 RR 3 41 12 1
2010-01-02 M 7908 RR 4 101 83 60
2010-01-02 VMAX RR 1 34 6 4
2010-01-02 VMAX RR 2 42 30 10
2010-01-02 VMAX RR 3 0 4 2
2010-01-02 VMAX RR 4 43 7 5
2010-01-02 CD 219 RR 1 166 39 33
2010-01-02 CD 219 RR 2 135 265 65
2010-01-02 CD 219 RR 3 26 43 1
2010-01-02 CD 219 RR 4 77 88 9
2010-01-11 BRS 245 RR 1 2 7 10
2010-01-11 BRS 245 RR 2 3 5 2
2010-01-11 BRS 245 RR 3 1 1 0
2010-01-11 BRS 245 RR 4 4 0 3
2010-01-11 BRS 243 RR 1 0 5 3
2010-01-11 BRS 243 RR 2 0 1 1
2010-01-11 BRS 243 RR 3 1 1 0
2010-01-11 BRS 243 RR 4 4 0 2
2010-01-11 BRS 246 RR 1 4 9 6
2010-01-11 BRS 246 RR 2 3 4 0
2010-01-11 BRS 246 RR 3 0 0 0
2010-01-11 BRS 246 RR 4 2 0 0
2010-01-11 BRS 239 1 1 10 2
2010-01-11 BRS 239 2 1 0 0
2010-01-11 BRS 239 3 4 3 1
2010-01-11 BRS 239 4 2 0 1
2010-01-11 EMBRAPA 48 1 0 2 1
2010-01-11 EMBRAPA 48 2 2 0 0
2010-01-11 EMBRAPA 48 3 3 1 2
2010-01-11 EMBRAPA 48 4 1 0 1
2010-01-11 CD 214 RR 1 2 0 1
2010-01-11 CD 214 RR 2 3 7 4
2010-01-11 CD 214 RR 3 4 16 8
2010-01-11 CD 214 RR 4 1 0 4
2010-01-11 CD 202 1 12 2 0
2010-01-11 CD 202 2 7 1 4
2010-01-11 CD 202 3 1 0 0
2010-01-11 CD 202 4 3 1 0
2010-01-11 M 7908 RR 1 2 9 3
2010-01-11 M 7908 RR 2 11 6 4
2010-01-11 M 7908 RR 3 4 1 1
2010-01-11 M 7908 RR 4 55 10 0
2010-01-11 VMAX RR 1 1 1 1
2010-01-11 VMAX RR 2 23 7 0
2010-01-11 VMAX RR 3 5 1 4
2010-01-11 VMAX RR 4 10 0 0
2010-01-11 CD 219 RR 1 65 12 4
2010-01-11 CD 219 RR 2 217 145 27
2010-01-11 CD 219 RR 3 10 2 1
2010-01-11 CD 219 RR 4 11 4 0
2010-01-18 BRS 245 RR 1 5 9 12
2010-01-18 BRS 245 RR 2 3 0 1
2010-01-18 BRS 245 RR 3 0 0 0
2010-01-18 BRS 245 RR 4 3 2 1
2010-01-18 BRS 243 RR 1 6 0 7
2010-01-18 BRS 243 RR 2 1 0 0
2010-01-18 BRS 243 RR 3 0 0 0
2010-01-18 BRS 243 RR 4 0 0 0
2010-01-18 BRS 246 RR 1 5 2 10
2010-01-18 BRS 246 RR 2 1 0 0
2010-01-18 BRS 246 RR 3 0 0 0
2010-01-18 BRS 246 RR 4 0 0 0
2010-01-18 BRS 239 1 3 1 0
2010-01-18 BRS 239 2 1 0 0
2010-01-18 BRS 239 3 0 0 0
2010-01-18 BRS 239 4 2 9 3
2010-01-18 EMBRAPA 48 1 1 6 4
2010-01-18 EMBRAPA 48 2 0 1 0
2010-01-18 EMBRAPA 48 3 0 0 0
2010-01-18 EMBRAPA 48 4 2 1 4
2010-01-18 CD 214 RR 1 5 2 4
2010-01-18 CD 214 RR 2 3 2 2
2010-01-18 CD 214 RR 3 0 0 2
2010-01-18 CD 214 RR 4 0 0 0
2010-01-18 CD 202 1 1 9 6
2010-01-18 CD 202 2 0 0 1
2010-01-18 CD 202 3 0 0 0
2010-01-18 CD 202 4 0 2 1
2010-01-18 M 7908 RR 1 3 1 13
2010-01-18 M 7908 RR 2 5 0 0
2010-01-18 M 7908 RR 3 1 0 0
2010-01-18 M 7908 RR 4 2 0 0
2010-01-18 VMAX RR 1 0 1 0
2010-01-18 VMAX RR 2 13 3 0
2010-01-18 VMAX RR 3 1 0 0
2010-01-18 VMAX RR 4 6 1 3
2010-01-18 CD 219 RR 1 49 12 2
2010-01-18 CD 219 RR 2 110 42 28
2010-01-18 CD 219 RR 3 7 3 0
2010-01-18 CD 219 RR 4 4 3 4
examples/ovos.txt 0 → 100644
+ 2101
− 0
View file @ 803b3ff3
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment