Profa. Dra. Suely Ruiz Giolo
Profa. Dra. Suely Ruiz Giolo
Códigos em Linguagem 
Códigos em Linguagem
Capítulo 2 - Delineamentos amostrais e modelos associados
2.1 Gráficos de colunas múltiplas e segmentadas - Figura 2.1 options(OutDec=",") # opção: valores numéricos com vírgulas
par(mfrow=c(1,2))
data<-cbind(E=c(75,45)/120, NE=c(21,56)/77)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","white"), ylim=range(c(0,1.19)), names.arg = c("Sim","Não"),
xlab="Exposição ao fator", ylab="Proporções amostrais",
legend.text = c("Doentes", "Sadios"),
args.legend = list(x = "topleft", bty="n", cex=1))
abline(h=0)
text(bp, c(0.64,0.38,0.28,0.73), round(data,3), cex=1, pos=3)
title("(a)")
mydata<-data.frame(row.names=c("Sim","Não"), Doentes=c(0.625,0.273), Sadios=c(0.375,0.727))
x <- barplot(t(as.matrix(mydata)), col=c("gray", "white"),
legend=TRUE, xlim=c(0,2.8), args.legend=list(x="topleft", bty="n", cex=1),
ylab="Proporções amostrais", xlab="Exposição ao fator", ylim=range(c(0,1.19)))
text(x, mydata$Doentes-0.15, labels=mydata$Doentes, col="black")
text(x, mydata$Doentes+0.20, labels=1-(mydata$Doentes))
abline(h=0)
title("(b)")
# Comandos p/ salvar a figura em extensão .eps e em escala de cinza
dev.copy2eps(file="Figura21.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray")
# Nota: há outras extensões disponíveis: png, tiff, jpeg, bmp etc.) 2.2 Gráficos de barras múltiplas e segmentadas - Figura 2.2 options(OutDec=",")
par(mfrow=c(1,2))
data<-cbind(E = c(75, 45)/120, NE = c(21, 56)/77)
bp<- barplot(height = data, beside = T, h = T,
col = c("gray","white"), xlim=range(c(0,1)),
names.arg = c("A", "B"), ylab="Tratamentos", xlab="Proporções amostrais",
legend.text = c("Curado", "Não curado"),
args.legend = list(x = "topright", bty="n", cex=1.2))
abline(v=0)
text(c(0.63,0.38,0.28,0.735), bp, round(data, 3), cex=1, pos=2)
title("(a)")
data<-cbind(C = c(75, 21)/96, CO = c(45, 56)/101)
bp<-barplot(height = data, beside = F, h = T,
col = c("gray", "white"), xlim=range(c(0,1.43)),
names.arg = c("Casos", "Controles"), ylab="Grupos",
xlab="Proporções amostrais",
legend.text = c("Expostos", "Não expostos"), width=1,
args.legend = list(x = "topright", bty="n", cex=1.2))
abline(v=0)
mydata <- data.frame(row.names=c("C","Co"), E=c(0.781,0.445), NE=c(0.219,0.555))
text(mydata$E-0.10, bp, labels=mydata$E, col="black")
text(c(mydata$E[1]+0.12, mydata$E[2]+0.45), bp, labels=1-(mydata$E))
title("(b)")
dev.copy2eps(file="Figura22.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray") 2.3 Gráficos de colunas e de setores - Figura 2.3 options(OutDec=",")
par(mfrow=c(1,2))
b<-c(75/197,45/197,21/197,56/197)
bp<-barplot(b, col= gray.colors(4), space=0, ylim=c(0,1), ylab="Proporções amostrais",
legend.text = c("Expostos e doentes", "Expostos e sadios", "Não expostos e doentes",
"Não expostos e sadios"), width = 1, args.legend = list(x = "topleft", bty="n", cex=1.2))
title("(a)")
text(bp, b, round(b,3), cex=1, pos=3)
abline(h=0)
slices<-c(0.381,0.228,0.107,0.284)
pielabels<-c("Expostos e \n doentes =","Expostos e \n sadios =",
"Não expostos e \n doentes =","Não expostos e \n sadios =")
pct<-slices*100
lbls<-paste(pielabels,pct)
lbls<-paste(lbls,"%", sep="")
pie(slices,labels = lbls, col=gray.colors(4), cex=1.2)
title("(b)")
dev.copy2eps(file="Figura23.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray") 2.4 Gráficos quadrúplo e mosaico - Figura 2.4 options(OutDec=",")
par(mfrow=c(1,2))
dados<-matrix(c(75,21,45,56), nc=2)
rownames(dados) <- c("exposto","não exposto")
colnames(dados) <- c("doente","sadio")
dados
fourfoldplot(dados, col=gray.colors(2), margin=2, conf.level=0)
title("(a)")
dados<-matrix(c(13,29,7,7,21,7), nc=3)
rownames(dados) <- c("Ativo","Placebo")
colnames(dados) <- c("Nenhuma","Alguma","Acentuada")
dados
mosaicplot(dados,main="", xlab="Tratamento", ylab="Melhora do paciente", col=gray.colors(3))
title("(b)")
dev.copy2eps(file="Figura24.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray")Capítulo 3 - Tabelas de contingência 2 × 2
3.1 Teste exato de Fisher (Tabela 3.2) dados<-matrix(c(2,18,6,14), nc=2)
dados
fisher.test(dados) 3.2 Exemplo sobre avaliação de medicamento (Tabela 3.4) # Gráfico de colunas múltiplas dos dados
options(OutDec=",")
par(mfrow=c(1,1))
data<-cbind(N = c(40, 20)/60, P = c(16, 48)/64)
bp<- barplot(height = data,
beside = TRUE,
col = c("gray","white"), ylim=range(c(0,1)),
names.arg = c("Novo", "Placebo"), xlab="Medicamentos", ylab="Proporções amostrais",
legend.text = c("Melhora: sim", "Melhora: não"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(0.67,0.33,0.25,0.75), round(data, 2), cex=1.2, pos=3)
# Teste qui-quadrado de Pearson
dados<-matrix(c(40,16,20,48), nc=2)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected
# Diferença entre proporções e IC(d)_95%
dados<-matrix(c(40,16,20,48), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
d<-p11-p21
vd<-((p11*(1-p11))/(sum(dados[1,])-1)) + ((p21*(1-p21))/(sum(dados[2,])-1))
dvd<-sqrt(vd)
z<-qnorm(0.975)
li<- d - (z*dvd)
ls<- d + (z*dvd)
cbind(d,li,ls)
# Risco relativo (RR) e IC(RR)_95%
dados<-matrix(c(40,16,20,48), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
RR<-p11/p21
vf1<-((1-p11)/(sum(dados[1,])*p11)) + ((1-p21)/(sum(dados[2,])*p21))
dpf1<-sqrt(vf1)
z<-qnorm(0.975)
li<-exp(log(RR) - z*dpf1)
ls<-exp(log(RR) + z*dpf1)
cbind(RR,li,ls)
# Razão de chances (OR) e IC(OR)_95%
dados<-matrix(c(40,16,20,48), nc=2)
dados
OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1])
vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2]))
dpf<-sqrt(vf)
z<-qnorm(0.975)
li<-exp(log(OR) - z*dpf)
ls<-exp(log(OR) + z*dpf)
cbind(OR,li,ls) 3.3 Exemplo sobre armadilhas de insetos (Figura 3.1) options(OutDec=",")
par(mfrow=c(1,1))
data<-cbind(AM = c(246, 17)/263, AL = c(458, 32)/490)
bp<- barplot(height = data,
beside = TRUE,
col = c("gray","white"), ylim=range(c(0,1.0)),
names.arg = c("Alaranjada", "Amarela"), xlab="Cor das armadilhas", ylab="Proporções amostrais",
legend.text = c("Insetos machos", "Insetos fêmeas"),
args.legend = list(x = "topright", bty="n", cex=1.2))
abline(h=0)
text(bp, c(246/263,17/263,458/490,32/490), round(data, 4), cex=1.2, pos=3)
# Teste qui-quadrado de Pearson
dados<-matrix(c(246,458,17,32), nc=2)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected
# Razão de chances (OR) e IC(OR)_95%
dados<-matrix(c(246,458,17,32), nc=2)
dados
OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1])
vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2]))
dpf<-sqrt(vf)
z<-qnorm(0.975)
li<-exp(log(OR) - z*dpf)
ls<-exp(log(OR) + z*dpf)
cbind(OR,li,ls) 3.4 Exemplo sobre tabagismo e câncer de pulmão (Tabela 3.5) # Gráfico de colunas múltiplas dos dados
options(OutDec=",")
par(mfrow=c(1,1))
data<-cbind(S = c(75, 45)/120, N = c(21, 56)/77)
bp<- barplot(height = data,
beside = TRUE,
col = c("gray","white"), ylim=range(c(0,1)),
names.arg = c("Sim", "Não"), xlab="Fumante", ylab="Proporções amostrais",
legend.text = c("Câncer de pulmão: sim", "Câncer de pulmão: não"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(75/120,45/120,21/77,56/77), round(data, 2), cex=1.2, pos=3)
# Teste qui-quadrado de Pearson
dados<-matrix(c(75,21,45,56), nc=2)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected
# Diferença entre proporções e IC(d)_95%
dados<-matrix(c(75,21,45,56), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
d<-p11-p21
vd<-((p11*(1-p11))/(sum(dados[1,])-1)) + ((p21*(1-p21))/(sum(dados[2,])-1))
dvd<-sqrt(vd)
z<-qnorm(0.975)
li<- d - (z*dvd)
ls<- d + (z*dvd)
cbind(d,li,ls)
# Risco relativo (RR) e IC(RR)_95%
dados<-matrix(c(75,21,45,56), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
RR<-p11/p21
vf1<-((1-p11)/(sum(dados[1,])*p11)) + ((1-p21)/(sum(dados[2,])*p21))
dpf1<-sqrt(vf1)
z<-qnorm(0.975)
li<-exp(log(RR) - z*dpf1)
ls<-exp(log(RR) + z*dpf1)
cbind(RR,li,ls)
# Razão de chances (OR) e IC(OR)_95%
dados<-matrix(c(75,21,45,56), nc=2)
dados
OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1])
vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2]))
dpf<-sqrt(vf)
z<-qnorm(0.975)
li<-exp(log(OR) - z*dpf)
ls<-exp(log(OR) + z*dpf)
cbind(OR,li,ls) 3.5 Exemplo sobre doenças respiratórias em crianças (Figura 3.2) # Gráficos de colunas múltiplas e de setores (pizza)
options(OutDec=",")
par(mfrow=c(1,2))
b<-c(0.329,0.379,0.116,0.176)
bp<-barplot(b, col= gray.colors(4), space=0, ylim=c(0,1), ylab="Proporções amostrais",
legend.text = c("Feminino com sintomas", "Masculino com sintomas",
"Feminino sem sintomas", "Masculino sem sintomas"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
text(bp, b, round(b, 3), cex=1, pos=3)
abline(h=0)
slices<-c(0.329,0.379,0.116,0.176)
pielabels<-c("Feminino com \n sintomas =", "Masculino com \n sintomas =",
"\n Feminino sem \n sintomas =", "Masculino sem \n sintomas =")
pct<-slices*100
lbls<-paste(pielabels, pct)
lbls<-paste(lbls,"%", sep="")
pie(slices, labels = lbls, col=gray.colors(4), cex=1.2)
# Teste qui-quadrado de Pearson
dados<-matrix(c(355,410,125,190), nc=2)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected
# Razão de chances (OR) e IC(OR)_95%
dados<-matrix(c(355,410,125,190), nc=2)
dados
OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1])
vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2]))
dpf<-sqrt(vf)
z<-qnorm(0.975)
li<-exp(log(OR) - z*dpf)
ls<-exp(log(OR) + z*dpf)
cbind(OR,li,ls) 3.6 Exemplo sobre medicamentos para infecções graves (Figura 3.3) # Gráfico de colunas múltiplas dos dados
options(OutDec=",")
par(mfrow=c(1,1))
data<-cbind(E = c(29, 16)/45, NE = c(14, 31)/45)
bp<- barplot(height = data,
beside = TRUE,
col = c("gray","white"), ylim=range(c(0,1)),
names.arg = c("Novo", "Padrão"), xlab="Medicamentos", ylab="Proporções amostrais",
legend.text = c("Resposta favorável", "Resposta não favorável"),
args.legend=list(x="topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(0.645,0.355,0.311,0.689), round(data,3), cex=1.2, pos=3)
# Teste qui-quadrado de Pearson
dados<-matrix(c(29,14,16,31), nc=2)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected
# Diferença entre proporções e IC(d)_95%
dados<-matrix(c(29,14,16,31), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
d<-p11-p21
vd<-((p11*(1-p11))/(sum(dados[1,])-1)) + ((p21*(1-p21))/(sum(dados[2,])-1))
dvd<-sqrt(vd)
z<-qnorm(0.975)
li<- d - (z*dvd)
ls<- d + (z*dvd)
cbind(d,li,ls)
# Risco relativo (RR) e IC(RR)_95%
dados<-matrix(c(29,14,16,31), nc=2)
dados
p11<-(dados[1,1]/(sum(dados[1,])))
p21<-(dados[2,1]/(sum(dados[2,])))
RR<-p11/p21
vf1<-((1-p11)/(sum(dados[1,])*p11)) + ((1-p21)/(sum(dados[2,])*p21))
dpf1<-sqrt(vf1)
z<-qnorm(0.975)
li<-exp(log(RR) - z*dpf1)
ls<-exp(log(RR) + z*dpf1)
cbind(RR,li,ls)
# Razão de chances (OR) e IC(OR)_95%
dados<-matrix(c(29,14,16,31), nc=2)
dados
OR<-(dados[1,1]*dados[2,2])/(dados[1,2]*dados[2,1])
vf<-(1/dados[1,1])+(1/dados[1,2])+(1/dados[2,1]+(1/dados[2,2]))
dpf<-sqrt(vf)
z<-qnorm(0.975)
li<-exp(log(OR) - z*dpf)
ls<-exp(log(OR) + z*dpf)
cbind(OR,li,ls)Nota: o pacote epiR também pode ser utilizado para obtenção de Qp, RR, d e OR
Instalar o package epiR (http://www.r-project.org/)
Consultar o manual epiR.pdf para mais detalhes require(epiR)
dados<-matrix(c(40,16,20,48), nc=2)
epi.2by2(dados, method="cohort.count", conf.level=0.95, units=1)
dados<-matrix(c(40,20,16,48), nc=2)
epi.2by2(dados, method="case.control", conf.level=0.95, units=1)
dados<-matrix(c(40,16,20,48), nc=2)
epi.2by2(dados, method="cross.sectional", conf.level=0.95, units=1)Capítulo 4 - Tabelas de contingência s × r
4.1 Ensaio clínico sobre artrite reumatoide - Figura 4.1 # Gráfico de colunas múltiplas dos dados
options(OutDec=",")
par(mfrow=c(1,1))
data<-cbind(A = c(13,7,21)/41, P = c(29,7,7)/43)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(3), ylim=range(c(0,1.19)),
names.arg = c("Ativo", "Placebo"), xlab="Tratamento", ylab="Proporções amostrais",
legend.text = c("Nenhuma melhora", "Alguma melhora", "Melhora acentuada"),
args.legend = list(x = "topleft", bty="n", cex=1))
abline(h=0)
text(bp, c(0.317,0.171,0.512,0.674,0.163,0.163), round(data, 3), cex=1, pos=3)
# Estatística QS em Tabela 2 x 3 (Y ordinal e ni+ fixos)
dados<-matrix(c(13,29,7,7,21,7), nc=3)
dados
escore<-c(0,1,2)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs = ((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1
p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p) 4.2 Estudo sobre uso de tabaco por adolescentes - Figura 4.2 # Gráfico de colunas múltiplas dos dados
par(mfrow=c(1,2))
b<-c(70,202,218,33,40,11)/574
bp<-barplot(b, col= gray.colors(6), space=0, ylim=c(0,1), ylab="Proporções amostrais",
legend.text = c("Consciência mínima e não usa tabaco",
"Consciência moderada e não usa tabaco",
"Consciência substancial e não usa tabaco",
"Consciência mínima e usa tabaco",
"Consciência moderada e usa tabaco",
"Consciência substancial e usa tabaco"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
title("(a)")
text(bp, b, round(b, 3), cex=1, pos=3)
abline(h=0)
# Gráfico de colunas condicional às marginais-linha
data<-cbind(Min = c(70, 33)/103, Mod = c(202,40)/242, Sub = c(218,11)/229)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(2), ylim=range(c(0,1.19)),
names.arg = c("Mínima", "Moderada", "Substancial"),
xlab="Consciência do risco", ylab="Proporções amostrais",
legend.text = c("Uso de tabaco: não", "Uso de tabaco: sim"),
args.legend = list(x = "topleft", bty="n", cex=1))
title("(b)")
abline(h=0)
text(bp, c(0.68,0.32,0.835,0.165,0.952,0.048), round(data, 3), cex=1, pos=3)
# Estatística QCs em Tabela 3 x 2 (X e Y ordinais e n fixo)
x<-c(rep(1,103), rep(2,242), rep(3,229))
y<-c(rep(0,70), rep(1,33), rep(0,202), rep(1,40), rep(0,218), rep(1,11))
rac<-cor(y,x)
n<-length(x)
QCS<-(n-1)*rac^2
p<-1-pchisq(QCS,1)
cbind(rac,QCS,p) 4.3 Estudo sobre afiliações político-partidárias (Tabela 4.4 - X e Y nominais) dados<-matrix(c(221,200,208,160,291,106,360,160,316,140,311,97), nc=4)
dados
Qp<-chisq.test(dados, correct=F)
Qp
Qp$observed
Qp$expected 4.4 Estudo sobre medicamentos para cefaleia (Tabela 4.5 - Y ordinal e ni+ fixos) dados<-matrix(c(6,1,2,9,4,5,6,6,6,3,6,8,1,8,6), nc=5)
dados
escore<-c(0,1,2,3,4)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
fb3<-(sum(dados[3,]*escore))/sum(dados[3,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3]),sum(dados[,4]),sum(dados[,5])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
sf<-sum(sum(dados[1,])*((fb1-mua)^2),sum(dados[2,])*((fb2-mua)^2),sum(dados[3,])*((fb3-mua)^2))
Qs<-((sum(dados)-1)*sf)/(sum(dados)*va)
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,fb3,Qs,gl,p)
# padrão versus novo
dados<-matrix(c(1,2,4,5,6,6,6,8,8,6), nc=5)
dados
escore<-c(0,1,2,3,4)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3]),sum(dados[,4]),sum(dados[,5])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs<-((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1
p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p)
# placebo versus novo
dados<-matrix(c(6,2,9,5,6,6,3,8,1,6), nc=5)
dados
escore<-c(0,1,2,3,4)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3]),sum(dados[,4]),sum(dados[,5])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs<-((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p)
# placebo versus padrão
dados<-matrix(c(6,1,9,4,6,6,3,6,1,8), nc=5)
dados
escore<-c(0,1,2,3,4)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3]),sum(dados[,4]),sum(dados[,5])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs<-((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p) 4.5 Estudo sobre produtos de limpeza - Tabela 4.6 (X e Y ordinais e ni+ fixos) x<-c(rep(1,46),rep(2,53),rep(3,67))
y<-c(rep(1,27),rep(2,14),rep(3,5),rep(1,10),rep(2,17),rep(3,26),rep(1,5),rep(2,12),rep(3,50))
rac<-cor(y,x)
n<-length(x)
QCS<-(n-1)*rac^2
p<-1-pchisq(QCS,1)
cbind(rac,QCS,p)
dados<-matrix(c(27,10,5,14,17,12,5,26,50), nc=3)
dados
escore<-c(1,2,3)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
fb3<-(sum(dados[3,]*escore))/sum(dados[3,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
sf<-sum(sum(dados[1,])*((fb1-mua)^2),sum(dados[2,])*((fb2-mua)^2),sum(dados[3,])*((fb3-mua)^2))
Qs<-((sum(dados)-1)*sf)/(sum(dados)*va)
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,fb3,Qs,gl,p)
# água versus água + dose única
dados<-matrix(c(27,10,14,17,5,26), nc=3)
dados
escore<-c(1,2,3)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs <- ((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p)
# água versus água + dose dupla
dados<-matrix(c(27,5,14,12,5,50), nc=3)
dados
escore<-c(1,2,3)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs <- ((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p)
# água + dose única versus água + dose dupla
dados<-matrix(c(10,5,17,12,26,50), nc=3)
dados
escore<-c(1,2,3)
fb1<-(sum(dados[1,]*escore))/sum(dados[1,])
fb2<-(sum(dados[2,]*escore))/sum(dados[2,])
esp<-(c(sum(dados[,1]),sum(dados[,2]),sum(dados[,3])))/sum(dados)
mua<-sum(escore*esp)
va<-sum((escore-mua)^2*esp)
vbf1<-((sum(dados) - sum(dados[1,]))/(sum(dados[1,])*(sum(dados)-1)))*va
Qs <- ((fb1-mua)^2)/vbf1
gl<-nrow(dados)-1; p<-1-pchisq(Qs,gl)
cbind(fb1,fb2,Qs,gl,p) 4.6 Estudo sobre anúncio publicitário de veículos - Tabela 4.7 dados<-matrix(c(4,0,5,0,3,5,0,3,2,2,4,2), nc=4)
dados
fisher.test(dados)Capítulo 5 - Análise estratificada
5.1 Ensaio clínico multicentros - Tabela 5.2 (Estatística QMH) tab<-array(c(29,14,16,31,37,24,8,21), dim=c(2,2,2))
mantelhaen.test(tab, correct=F)
Nota: QMH não é apropriado se as razões de chances não estiverem na mesma direção.
# Teste de Breslow-Day para testar H0: homogeneidade das razões de chances
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/breslowday.test.R.txt")
tab<-array(c(29,14,16,31,37,24,8,21), dim=c(2,2,2))
breslowday.test(tab) 5.2 Ensaio clínico duplo cego - Tabela 5.5 (Estatística QSMH) dados<-matrix(c(6,19,7,9,5,7,2,1,16,6,5,1), nc=3)
dados
escore<-c(0,1,2)
fb11<-(sum(dados[1,]*escore))/sum(dados[1,])
fb21<-(sum(dados[3,]*escore))/sum(dados[3,])
fb12<-(sum(dados[2,]*escore))/sum(dados[2,])
fb22<-(sum(dados[4,]*escore))/sum(dados[4,])
fm1<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(fb11,fb21))
esp1<-(c(sum(dados[1:2,1]),sum(dados[1:2,2]),sum(dados[1:2,3])))/sum(dados[1:2,])
mu1<-sum(escore*esp1)
esp2<-(c(sum(dados[3:4,1]),sum(dados[3:4,2]),sum(dados[3:4,3])))/sum(dados[3:4,])
mu2<-sum(escore*esp2)
mu<-sum(c(sum(dados[1,]),sum(dados[3,]))*c(mu1,mu2))
v1<- sum(((escore-mu1)^2)*esp1)
v2<- sum(((escore-mu2)^2)*esp2)
vfma<-(sum(dados[1,])*sum(dados[2,])*v1)/(sum(dados[1:2,])-1)
vfmb<-(sum(dados[3,])*sum(dados[4,])*v2)/(sum(dados[3:4,])-1)
vfm<- sum(c(vfma,vfmb))
QSMH<-((fm1-mu)^2)/vfm; p<-1-pchisq(QSMH,1)
round(c(QSMH,p),digits=5)
c(fb11,fb21,fb12,fb22)Capítulo 6 - Tabelas com dados relacionados
6.1 Teste de McNemar aplicado aos dados dispostos na Tabela 6.1 dados<-matrix(c(20,10,5,10), nc=2)
mcnemar.test(dados,correct=F) 6.2 Nomograma de Fagan # https://achekroud.github.io/nomogrammer_vignette.html
source("https://raw.githubusercontent.com/achekroud/nomogrammer/master/nomogrammer.r")
library(ggplot2) # instalar
library(scales) # instalar
nomogrammer(Prevalence = .371, Sens=0.846, Spec=0.909)
scen1<-nomogrammer(Prevalence = .371, Sens=0.846, Spec=0.909)
nomogrammer(Prevalence = .371, Sens=0.846, Spec=0.909) + ggtitle("Cenário 1") +
scale_color_manual(values=c("darkgray", "black"))
nomogrammer(Prevalence = .055, Sens=0.846, Spec=0.901) + ggtitle("Cenário 2") +
scale_color_manual(values=c("darkgray", "black"))
* Nota: nomogramas e os IC a seguir também podem ser obtidos em http://araw.mede.uic.edu/cgi-bin/testcalc.pl 6.2.1 Intervalos de confiança - Cenário 1 # IC(LR+)_95% = exp{log(LR+) +/- sqrt(var(LR+))}
a1 <- 0.846 # sensibilidade
a2 <- 0.909 # especificidade
LRp <- a1/(1-a2) # LRp = LR+
v1<-((1-a1)/(65*a1))+((a2)/(110*(1-a2))) # var(LR+) = [(1-a1)/n1+*a1] + [(a2)/(n2+*(1-a2)]
sv1<- sqrt(v1)
ln <-log(LRp)
Lrs<-exp(ln + 1.96*sv1)
Lri<-exp(ln - 1.96*sv1)
cbind(LRp,Lri,Lrs) # (LRp, IC(LR+)) = (LR+, Li(LR+); Ls(LR+))
# IC(p2)_95% = (c2i/(1+c2i); c2s/(1+c2s))
p1<-65/175 # p1 = prevalência
c1<-(p1/(1-p1))
c2<-c1*LRp
p2<-c2/(1+c2) # estimativa pontual de p2 = probabilidade a posteriori
c2i<-c1*Lri
Li<-c2i/(1+c2i)
c2s<-c1*Lrs
Ls<- c2s/(1+c2s)
cbind(p2,Li,Ls) # (p2, IC(p2)) = (p2, Li, Ls)
# IC(LR-)_95% = exp{log(LR-) +/- sqrt(var(LR-))}
a1 <- 0.846 # sensibilidade
a2 <- 0.909 # especificidade
LRm <- (1-a1)/a2 # LRm = LR-
v1<-((a1)/(65*(1-a1)))+((1-a2)/(110*a2)) # var(LR-) = [(a1)/n1+*(1-a1)] + [(1-a2)/(n2+*a2]
sv1<- sqrt(v1)
ln <-log(LRm)
Lrs<-exp(ln + 1.96*sv1)
Lri<-exp(ln - 1.96*sv1)
cbind(LRm,Lri,Lrs) # (LRm, IC(LR-)) = (LR-, Li(LR-); Ls(LR-))
# IC(p2)_95% = (c2i/(1+c2i); c2s/(1+c2s))
p1<-65/175 # p1 = prevalência
c1<-p1/(1-p1)
c2<-c1*LRm
p2<-c2/(1+c2) # estimativa pontual de p2 = probabilidade a posteriori
c2i<-c1*Lri
Li<-c2i/(1+c2i)
c2s<-c1*Lrs
Ls<- c2s/(1+c2s)
cbind(p2,Li,Ls) # (p2, IC(p2)) = (p2, Li, Ls)
# Nota: em http://araw.mede.uic.edu/cgi-bin/testcalc.pl os IC's são fornecidos diretamente nas saídas 6.3 Estatística capa aplicada aos dados dispostos na Tabela 6.4 require(vcd) # instalar o package vcd (http://www.r-project.org)
x<-c(38,5,0,1,33,11,3,0,10,14,5,6,3,7,3,10)
x<-matrix(x,4,4)
Kappa(x)
summary(Kappa(x))
confint(Kappa(x))
confint(Kappa(x), level=0.90) 6.4 Gráfico de concordância (agreement plot) require(vcd)
x<-c(38,5,0,1,33,11,3,0,10,14,5,6,3,7,3,10)
x<-matrix(x,4,4)
rownames(x) <- c("1","2","3","4")
colnames(x) <- c("1","2","3","4")
x
agreementplot(x, xlab = "Estágios: neurologista 1",
ylab = "Estágios: neurologista 2",
line_col="gray")Capítulo 7 - Regressão binomial
7.1 Estudo sobre doença coronária versus idade (dados na Tabela 7.1) resim<-c(1,2,3,5,6,5,13,8)
resnao<-c(9,13,9,10,7,3,4,2)
idade<-c(25,32,38,43,47,53,57,65)
dados<-as.data.frame(cbind(resim,resnao,idade))
dados
attach(dados)
theta<-resim/(resim+resnao)
plot(idade, theta, ylim=range(0,0.9), xlab="Idade", ylab="E(Y|x)", pch=16)
ajust<-glm(cbind(resim,resnao)~idade, family=binomial(link="logit"), data=dados)
ajust
anova(ajust, test="Chisq")
summary(ajust)
cbind(ajust$fitted.values, ajust$y)
ajust$residuals
dev<-residuals(ajust, type='deviance')
dev
QL<-sum(dev^2)
p1<-1-pchisq(QL,6)
cbind(QL,p1)
rpears<-residuals(ajust, type='pearson')
rpears
QP<-sum(rpears^2)
p2<-1-pchisq(QP,6)
cbind(QP,p2)
theta<-resim/(resim+resnao)
plot(idade, theta, ylim=range(0,0.9), xlab="Idade", ylab="E(Y|x)", pch=16)
idade<-20:70
modajust<-(exp(-5.123+0.1058*idade))/(1+exp(-5.123+0.1058*idade))
lines(idade, modajust)
par(mfrow=c(1,2))
plot(rpears, ylab="Resíduos Pearson", pch=16, ylim=c(-2,2))
abline(h=0, lty=3)
plot(dev, ylab="Resíduos deviance", pch=16, ylim=c(-2,2))
abline(h=0, lty=3)
par(mfrow=c(1,1))
ntot<-c(10,15,12,15,13,8,17,10)
fit.model<-ajust
source("http://www.ime.usp.br/~giapaula/envelr_bino")
rm(list=ls())
dados1<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/coronaria.txt",h=T)
attach(dados1) # (dados1 = indivíduos por linha)
# > install.packages("Epi", dependencies=TRUE)
require(Epi)
ROC(form=y~idade, plot="ROC") 7.2 Estudo sobre doença coronária versus sexo e ECG (Tabela 7.7) resim<-c(4,8,9,21); resnao<-c(11,10,9,6)
sexo<-c(0,0,1,1); ecg<-c(0,1,0,1)
dados<-as.data.frame(cbind(resim, resnao, sexo, ecg))
attach(dados)
ajust1<- glm(cbind(resim,resnao)~sexo+ecg+sexo*ecg, family=binomial(link="logit"), data=dados)
summary(ajust1)
anova(ajust1, test="Chisq")
ajust<-glm(cbind(resim, resnao)~sexo+ecg, family=binomial(link="logit"), data=dados)
summary(ajust)
anova(ajust, test="Chisq")
cbind(ajust$y, ajust$fitted.values, ajust$residuals)
dev<-residuals(ajust, type='deviance'); dev
QL<-sum(dev^2); p1<-1-pchisq(QL,1); cbind(QL,p1)
rpears<-residuals(ajust, type='pearson'); rpears
QP<-sum(rpears^2); p2<-1-pchisq(QP,1); cbind(QP,p2) 7.3 Estudo sobre infecções urinárias (Tabela 7.12) rm(list=ls())
resim<-c(78,101,68,40,54,34)
resnao<-c(28,11,46,5,5,6)
diag<-c(1,1,1,0,0,0)
trat<-c(1,2,3,1,2,3)
dados<-as.data.frame(cbind(resim, resnao, diag, trat))
attach(dados)
ajust1<-glm(cbind(resim,resnao)~diag+relevel(factor(trat),3)+diag*(relevel(factor(trat),3)),
family=binomial(link="logit"), data=dados)
summary(ajust1)
anova(ajust1, test="Chisq")
ajust2<-glm(cbind(resim,resnao)~diag+relevel(factor(trat),3), family=binomial(link="logit"),
data=dados)
summary(ajust2)
anova(ajust2, test="Chisq")
cbind(ajust2$fitted.values, ajust2$y)
dev<-residuals(ajust2, type='deviance'); dev
QL<-sum(dev^2)
p1<-1-pchisq(QL,2)
cbind(QL,p1)
rpears<-residuals(ajust2, type='pearson'); rpears
QP<-sum(rpears^2)
p2<-1-pchisq(QP,2)
cbind(QP,p2)
par(mfrow=c(1,2))
plot(dev, pch=16, ylim=c(-3,3), ylab="Resíduos deviance")
abline(h=0, lty=3)
plot(rpears, pch=16, ylim=c(-3,3), ylab="Resíduos Pearson")
abline(h=0, lty=3)
par(mfrow=c(1,1))
ntot<-c(106,112,114,45,59,40)
fit.model<-ajust2
source("http://www.ime.usp.br/~giapaula/envelr_bino")
dados1<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/infec.txt",h=T)
# dados1 = arquivo indivíduo por linha (476 x 3)
attach(dados1)
require(Epi)
ROC(form=y~factor(x1)+relevel(factor(x2),3), plot="ROC") 7.4 Estudo sobre bronquite (Tabela 7.18) brc<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/bronquite.txt",h=T)
attach(brc)
brc
ajust<-glm(cbind(sim, nao)~smk+ses+idade+smk*ses+smk*idade+ses*idade+smk*ses*idade,
family=binomial, data=brc)
anova(ajust, test="Chisq")
ajust<-glm(cbind(sim, nao)~smk+ses+idade+smk*ses+smk*idade+ses*idade, family=binomial, data=brc)
anova(ajust, test="Chisq")
ajust<-glm(cbind(sim, nao)~smk+ses+idade+smk*ses+smk*idade, family=binomial, data=brc)
anova(ajust, test="Chisq")
ajust<-glm(cbind(sim, nao)~smk+ses+idade+smk*idade, family=binomial, data=brc)
anova(ajust, test="Chisq")
ajust<-glm(cbind(sim, nao)~smk+ses+idade+smk*ses+smk*idade, family=binomial, data=brc)
anova(ajust, test="Chisq")
summary(ajust)
cbind(ajust$fitted.values, ajust$y)
dev<-residuals(ajust, type='deviance')
dev
QL<-sum(dev^2)
p1<-1-pchisq(QL,6)
cbind(QL,p1)
rpears<-residuals(ajust, type='pearson')
rpears
QP<-sum(rpears^2)
p2<-1-pchisq(QP,6)
cbind(QP,p2)
par(mfrow=c(1,2))
plot(rpears, pch=16, ylim=c(-0.5,0.5), ylab="Residuos de Pearson")
abline(h=0, lty=3)
plot(dev, pch=16, ylim=c(-0.5,0.5), ylab="Residuos deviance")
abline(h=0, lty=3)
par(mfrow=c(1,1))
ntot<-c(111,134,95,124,173,148,143,112)
fit.model<-ajust
source("http://www.ime.usp.br/~giapaula/envelr_bino")
brc1<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/bronquite1.txt",h=T)
# brc1 = arquivo com 1 indivíduo por linha (1040 x 4)
attach(brc1)
require(Epi)
ROC(form=y~x1+x2+x3+x1*x2+x1*x3, plot="ROC")
# Procedimentos de seleção de variáveis: forward, backward e stepwise
modp<-glm(as.matrix(brc[,c(1,2)])~smk+ses+idade+smk*ses+smk*idade+ses*idade+smk*ses*idade,
family=binomial, data=brc)
mod0<-glm(as.matrix(brc[,c(1,2)])~1, family=binomial, data=brc)
step(mod0, ~smk+ses+idade+smk*ses+smk*idade+ses*idade+smk*ses*idade, direction=c("forward"))
step(modp, ~smk+ses+idade+smk*ses+smk*idade+ses*idade+smk*ses*idade, direction=c("backward"))
step(mod0, ~smk+ses+idade+smk*ses+smk*idade+ses*idade+smk*ses*idade, direction=c("both")) 7.5 Outro estudo sobre doença coronária (Tabela 7.22) dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/chd4.txt",h=T)
dados<-as.data.frame(dados)
attach(dados); head(dados)
ajust1<-glm(cbind(dcnao,dcsim)~sexo+ecg+idade+sexo*ecg+sexo*idade+ecg*idade+sexo*ecg*idade,
family=binomial(link="logit"), data=dados)
summary(ajust1)
anova(ajust1, test="Chisq")
ajust2<-glm(cbind(dcnao,dcsim)~sexo+ecg+idade, family=binomial, data=dados)
summary(ajust2)
anova(ajust2, test="Chisq")
cbind(sexo, ecg, idade, ajust2$fitted.values)
dev<-residuals(ajust2, type='deviance'); dev
rpears<-residuals(ajust2, type='pearson'); rpears
par(mfrow=c(1,2))
plot(rpears, pch=16, ylab="Resíduos de Pearson")
plot(dev, pch=16, ylab="Resíduos deviance")
# Estatística QHL
dados1<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/chd4a.txt", h=T)
attach(dados1)
dados1
ajust2a<-glm(dc~sexo+ecg+idade, family=binomial(link="logit"), data=dados1)
summary(ajust2a)
anova(ajust2a, test="Chisq")
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/gof_bino.txt")
gof.bino(ajust2a, grupos=10)
# ou source("http://www.poleto.com/funcoes/gof.bino.txt")
# Autor da função gof.bino: Poleto, F.Z. (http://www.poleto.com)
# QQ-plot com envelope simulado
fit.model<-ajust2a
attach(dados1)
source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envel_bino.txt")
# ou > source("http://www.ime.usp.br/~giapaula/envel_bino")
# Nota: use o código envelr.bino.txt para dados com repetições
# Curva ROC
# Opção 1: Instalar o package ROCR (http://www.r-project.org/)
require(ROCR)
pred <- prediction(ajust2a$fitted.values,dc); pred
perf <- performance(pred,"tpr","fpr")
perf # tpr = sensibilidade e fpr = 1 - especificidade
plot(perf); abline(c(0,0),c(1,1),lty=3)
points(0.324,0.756, pch=20, col=2)
area <- performance(pred,"auc"); area
# Opção 2: Instalar o package Epi (http://www.r-project.org/)
ajust2a<-glm(dc~sexo+ecg+idade, family=binomial(link="logit"), data=dados1)
summary(ajust2a)
require(Epi)
ROC(form=dc~sexo+ecg+idade, plot="ROC") 7.6 Modelos para dados binários: links alternativos (Tabela 7.28) # Ajuste de 4 modelos aos dados de um bioensaio
resim<-c(0,1,4,15,28,29)
resnao<-c(30,29,26,15,2,1)
lnd<-c(0,9.21,16.12,18.42,20.72,23.02)
dados<-as.data.frame(cbind(resim,resnao,lnd))
dados
attach(dados)
ajuste1<-glm(cbind(resim,resnao)~lnd, family=binomial(link="logit"), data=dados)
ajuste1
anova(ajuste1, test="Chisq")
summary(ajuste1)
ajuste2<-glm(cbind(resim,resnao)~lnd, family=binomial(link="probit"), data=dados)
ajuste2
anova(ajuste2, test="Chisq")
summary(ajuste2)
ajuste3<-glm(cbind(resim,resnao)~lnd, family=binomial(link="cloglog"), data=dados)
ajuste3
anova(ajuste3, test="Chisq")
summary(ajuste3)
ajuste4<-glm(cbind(resim,resnao)~lnd, family=binomial(link="cauchit"), data=dados)
ajuste4
anova(ajuste4, test="Chisq")
summary(ajuste4)
# QQ-plot com envelope simulado
par(mfrow=c(1,4))
ntot<-c(30,30,30,30,30,30)
fit.model<-ajuste1
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste2
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste3
source("http://www.ime.usp.br/~giapaula/envelr_bino")
fit.model<-ajuste4
source("http://www.ime.usp.br/~giapaula/envelr_bino")
# ou source("https://docs.ufpr.br/~giolo/LivroADC/Funcoes/envelr_bino.txt")
# Gráfico com as curvas associadas aos modelos ajustados
x<-seq(0,25,0.1)
m1<-exp(-12.863+0.708*x)/(1+exp(-12.863+0.708*x))
m2<-pnorm(-6.244+0.347*x)
m3<-1-exp(-exp(-8.143+0.422*x))
m4<-0.5+(atan(-26.678+1.451*x)/pi) # ou pcauchy(-26.678+1.451*x)
par(mfrow=c(1,1))
plot(lnd,resim/(resim+resnao), pch=16, ylab="Proporção de mortes",
xlab="Ln(diluições)", xlim=c(0,28), ylim=c(0,1.05))
lines(x, m1, lty=2, lwd=2, col=2)
lines(x, m2, lty=3, lwd=2, col=4)
lines(x, m3, lty=6, lwd=2, col=3)
lines(x, m4, lty=1, lwd=2, col=1)
legend(1,0.8, lty=c(2,3,6,1), col=c(2,4,3,1), lwd=2, c("logístico","probito","clog-log","Cauchy"),
bty="n")
lines(c(18.386,18.386),c(0,0.50), lty=3)
lines(c(0,18.386),c(0.50,0.50), lty=3)
legend(17.7,0.55,c("(18.386, 0.5)"), bty="n", cex=0.8)Capítulo 8 - Regressão multinomial
8.1 Estudo sobre aprendizado (Tabela 8.1) - Modelo logitos categoria de referência (MLCR) # 8.1.1 Representação gráfica dos dados
par(mfrow=c(1,3))
data<-cbind(P = c(10,17,26)/53, I = c(5,12,50)/67)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","lightgray","white"), ylim=range(c(0,1)),
names.arg = c("Padrão", "Integral"), xlab="Escola 1", ylab="Proporções amostrais",
legend.text = c("Individual", "Grupo", "Sala de aula"),
args.legend = list(x = "topleft", bty="n", cex=1.4))
abline(h=0)
text(bp, c(0.189,0.320,0.491,0.075,0.179,0.746), round(data,3), cex=1.4, pos=3)
data<-cbind(P = c(21,17,26)/64, I = c(16,12,36)/64)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","lightgray","white"), ylim=range(c(0,1)),
names.arg = c("Padrão", "Integral"), xlab="Escola 2", ylab="Proporções amostrais",
legend.text = c("Individual", "Grupo", "Sala de aula"),
args.legend = list(x = "topleft", bty="n", cex=1.4))
abline(h=0)
text(bp, c(0.328,0.266,0.406,0.250,0.187,0.563), round(data,3), cex=1.4, pos=3)
data<-cbind(P = c(15,15,16)/46, I = c(12,12,20)/44)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","lightgray","white"), ylim=range(c(0,1)),
names.arg = c("Padrão", "Integral"), xlab="Escola 3", ylab="Proporções amostrais",
legend.text = c("Individual", "Grupo", "Sala de aula"),
args.legend = list(x = "topleft", bty="n", cex=1.4))
abline(h=0)
text(bp, c(0.326,0.326,0.348,0.273,0.273,0.454), round(data,3), cex=1.4, pos=3)
# 8.1.2 Ajuste dos modelos
# Nota: Instalar o pacote VGAM (http://www.r-project.org/)
require(VGAM)
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/escolas.txt", h=T)
mlcr0<-vglm(cbind(ind,grupo,sala)~1, multinomial, dados)
mlcr0
mlcr1<-vglm(cbind(ind,grupo,sala)~factor(escola), multinomial, dados)
mlcr1
mlcr2<-vglm(cbind(ind,grupo,sala)~factor(escola)+factor(periodo), multinomial, dados)
mlcr2
mlcr3<-vglm(cbind(ind,grupo,sala)~factor(escola)+factor(periodo)+factor(escola)*factor(periodo),
multinomial, dados)
mlcr3
AIC(mlcr0)
AIC(mlcr1)
AIC(mlcr2)
AIC(mlcr3)
# Nota: (a) modelo com interação = modelo saturado, logo deviance residual = 0
# (b) função anova não está disponível em vglm( ); construir a partir das saídas.
deviance(mlcr2)
df.residual(mlcr2)
summary(mlcr2)
coef(mlcr2, matrix=TRUE)
mlcr2@y
fitted(mlcr2)
mlcr2@y - fitted(mlcr2)
par(mfrow=c(1,2))
rp<-resid(mlcr2, type = "pearson")
plot(rp[,1], pch=16, ylim=c(-3,3), xlab="Índice", ylab="Resíduos de Pearson")
title("Logito 1")
abline(h=0, lty=3)
plot(rp[,2], pch=20, ylim=c(-3,3), xlab="Índice", ylab="Resíduos de Pearson")
abline(h=0, lty=3)
title("Logito 2") 8.2 Estudo sobre artrite (Tabela 8.6) - Modelo logitos cumulativos (MLC) # 8.2.1 Representação gráfica dos dados
par(mfrow=c(1,2))
data<-cbind(P = c(16,5,6)/27, I = c(6,7,19)/32)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","lightgray", "white"), ylim=range(c(0,1.1)),
names.arg = c("Tratamento A", "Placebo"), xlab="Feminino", ylab="proporções amostrais",
legend.text = c("Acentuada", "Alguma", "Nenhuma"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp,c(16/27, 5/27, 6/27, 6/32, 7/32, 19/32), round(data,3), cex=1.2, pos=3)
data<-cbind(P = c(5,2,7)/14, I = c(1,1,9)/11)
bp<- barplot(height = data, beside = TRUE,
col = c("gray","lightgray", "white"), ylim=range(c(0,1.1)),
names.arg = c("Tratamento A", "Placebo"), xlab="Masculino", ylab="proporções amostrais",
legend.text = c("Acentuada", "Alguma", "Nenhuma"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp,c(5/14, 2/14, 7/14, 1/11, 1/11, 9/11), round(data,3), cex=1.2, pos=3)
# 8.2.2 Ajuste e seleção de modelos
# (a) Suposição de chances proporcionais: observação das estimativas
require(VGAM)
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/artrite.txt", h=T)
mlc<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=FALSE,reverse=FALSE), dados)
mop<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=TRUE, reverse=FALSE), dados)
coef(mlc, matrix = TRUE)
coef(mop, matrix = TRUE)
# Nota: inclua a opção link para utilizar outra função de ligação (que não a logito). Isto é,
# cumulative(parallel = T, reverse = F, link = "probit")
# (b) Suposição de chances proporcionais: teste da razão de verossimilhanças
## Modelo com X1 e X2 (sexo e tratamento)
mop<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=T,reverse=F), dados)
mlc<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=F,reverse=F), dados)
TRV<- 2*(logLik(mlc)-logLik(mop))
gl <- length(coef(mlc))-length(coef(mop)); p<-1-pchisq(TRV,gl)
cbind(TRV, gl, p)
## ou equivalentemente
TRV<- deviance(mop)-deviance(mlc)
gl <- df.residual(mop)-df.residual(mlc)
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
# Modelo com X1 (sexo)
mop<-vglm(cbind(ac,alg,nenh)~factor(sexo), cumulative(parallel=T,reverse=F), dados)
mlc<-vglm(cbind(ac,alg,nenh)~factor(sexo), cumulative(parallel=F,reverse=F), dados)
TRV<-2*(logLik(mlc)-logLik(mop)); gl<-length(coef(mlc))-length(coef(mop))
p<-1-pchisq(TRV,gl)
cbind(TRV, gl, p)
# Modelo com X2 (tratamento)
mop<-vglm(cbind(ac,alg,nenh)~factor(trat), cumulative(parallel=T,reverse=F), dados)
mlc<-vglm(cbind(ac,alg,nenh)~factor(trat), cumulative(parallel=F,reverse=F), dados)
TRV<-2*(logLik(mlc)-logLik(mop)); gl<-length(coef(mlc))-length(coef(mop))
p<-1-pchisq(TRV,gl)
cbind(TRV, gl, p)
#### Escolha: MOP ####
mop0<-vglm(cbind(ac,alg,nenh)~1, cumulative(parallel=T,reverse=F), dados)
mop0
mop1<-vglm(cbind(ac,alg,nenh)~factor(sexo), cumulative(parallel=T,reverse=F), dados)
mop1
mop2<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=T,reverse=F), dados)
mop2
mop3<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat)+factor(sexo)*factor(trat),
cumulative(parallel=T,reverse=F), dados)
mop3
AIC(mop0)
AIC(mop1)
AIC(mop2)
AIC(mop3)
TRV<- deviance(mop2)-deviance(mop3)
gl <- df.residual(mop2)-df.residual(mop3)
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
TRV<- deviance(mop0)-deviance(mop1)
gl <- df.residual(mop0)-df.residual(mop1)
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
TRV<- deviance(mop1)-deviance(mop2)
gl <- df.residual(mop1)-df.residual(mop2)
p <- 1-pchisq(TRV,gl)
cbind(TRV, gl, p)
# 8.2.3 MOP selecionado: com X1 e X2
mop2<-vglm(cbind(ac,alg,nenh)~factor(sexo)+factor(trat), cumulative(parallel=T,reverse=F), dados)
summary(mop2)
coef(mop2,matrix = TRUE)
mop2@y
fitted(mop2)
mop2@y - fitted(mop2)
QL<-deviance(mop2); QL
1-pchisq(QL,4)
rp<-residuals(mop2, type="pearson")
Qp<-sum(rp[,1]^2) + sum(rp[,2]^2)
Qp
1-pchisq(Qp,4)
rp<-resid(mop2, type = "pearson")
plot(1:4,rp[,1], pch=20, ylim=c(-3,3), xlim=c(1,8), xlab="Índice", ylab="Resíduos de Pearson")
points(5:8, rp[,2], pch=20)
abline(h=0, lty=3) 8.3 Estudo sobre tratamentos para artrite (Tabela 8.12) - Modelo logitos categorias adjacentes (MLCA) # Representação gráfica dos dados
par(mfrow=c(1,2))
data<-cbind(M = c(25,30,46,27)/128, F = c(7,15,35,13)/70)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(4), ylim=range(c(0,1)),
names.arg = c("Masculino", "Feminino"), xlab="Tratamento A", ylab="Proporções amostrais",
legend.text = c("Remissão total", "Remissão parcial", "Nenhuma mudança", "Progressão da doença"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(25/128,30/128,46/128,27/128, 7/70,15/70,35/70,13/70), round(data, 3),cex=1.2, pos=3)
data<-cbind(M = c(20,19,44,42)/125, F = c(8,10,21,36)/75)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(4), ylim=range(c(0,1)),
names.arg = c("Masculino", "Feminino"), xlab="Tratamento B", ylab="Proporções amostrais",
legend.text = c("Remissão total", "Remissão parcial", "Nenhuma mudança", "Progressão da doença"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(20/125, 19/125, 44/125, 42/125, 8/75, 10/75,21/75,36/75), round(data, 3), cex=1.2, pos=3)
dev.copy2eps(file="Figura85.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray")
# Ajuste e seleção de modelos
require(VGAM)
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/trat_ar.txt", h=T)
attach(dados)
dados
mca1<-vglm(cbind(rt,rp,nm,pd)~ factor(trat)+factor(sexo), family=acat(reverse=T,parallel=T), dados)
mca2<-vglm(cbind(rt,rp,nm,pd)~ factor(trat)+factor(sexo), family=acat(reverse=T,parallel=F), dados)
coef(mca1, matrix = TRUE)
coef(mca2, matrix = TRUE)
TRV<- deviance(mca1)-deviance(mca2)
gl <- df.residual(mca1)-df.residual(mca2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# Modelo com X1 (tratamento)
mca1<-vglm(cbind(rt,rp,nm,pd)~factor(trat), acat(parallel=T,reverse=T), dados)
mca2<-vglm(cbind(rt,rp,nm,pd)~factor(trat), acat(parallel=F,reverse=T), dados)
TRV<- deviance(mca1)-deviance(mca2)
gl <- df.residual(mca1)-df.residual(mca2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# Modelo com X2 (sexo)
mca1<-vglm(cbind(rt,rp,nm,pd)~factor(sexo), acat(parallel=T,reverse=T), dados)
mca2<-vglm(cbind(rt,rp,nm,pd)~factor(sexo), acat(parallel=F,reverse=T), dados)
TRV<- deviance(mca1)-deviance(mca2)
gl <- df.residual(mca1)-df.residual(mca2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# MLCA com chances proporcionais parciais (X2 com chances proporcionais e X1 não)
mca0<-vglm(cbind(rt,rp,nm,pd)~1, acat(parallel=F~factor(trat),reverse=T), dados)
AIC(mca0)
mca0
mca1<-vglm(cbind(rt,rp,nm,pd)~factor(trat), acat(parallel=F~factor(trat),reverse=T), dados)
mca1
AIC(mca1)
mca2<-vglm(cbind(rt,rp,nm,pd)~factor(trat)+factor(sexo), acat(parallel=F~factor(trat),reverse=T), dados)
AIC(mca2)
mca2
mca3<-vglm(cbind(rt,rp,nm,pd)~factor(trat)+factor(sexo)+factor(trat)*factor(sexo),
acat(parallel= F~factor(trat),reverse=T), dados)
AIC(mca3)
mca3
# Selecionado: MLCA com chances proporcionais parciais (X2 com chances proporcionais e X1 não)
summary(mca2)
coef(mca2,matrix = TRUE)
mca2@y
fitted(mca2)
mca2@y - fitted(mca2)
rp<-resid(mca2, type = "pearson")
Qp<-sum(rp[,1]^2) + sum(rp[,2]^2)+ sum(rp[,3]^2)
cbind(Qp,1-pchisq(Qp,5))
QL<-deviance(mca2)
cbind(QL,1-pchisq(QL,5))
par(mfrow=c(1,2))
rb<- mca2@y - fitted(mca2)
plot(rb[1,],ylim=c(-0.4,0.4),xlim=c(1,16),xlab="Índice",ylab="Probabilidade observada - predita",pch=16)
points(5:8, rb[2,], pch=16); points(9:12,rb[3,], pch=16)
points(13:16, rb[4,], pch=16)
abline(h=0, lty=3)
title("(a)")
rp<-resid(mca2, type = "pearson")
plot(1:4,rp[,1], pch=16, ylim=c(-3,3),xlim=c(1,12), xlab="Índice", ylab="Resíduos de Pearson")
points(5:8, rp[,2], pch=16)
points(9:12, rp[,3], pch=16)
abline(h=0, lty=3)
title("(b)") 8.4 Estudo sobre métodos de analgesia (Tabela 8.17) - Modelo logitos razão contínua (MLRC) # Representação gráfica dos dados
par(mfrow=c(1,2))
data<-cbind(A = c(9,19,21,81,537)/667, B = c(17,33,37,134,445)/666)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(5), ylim=range(c(0,1)),
names.arg = c("Método A", "Método B"), xlab="Feminino", ylab="Proporções amostrais",
legend.text = c("Intolerável", "Intensa", "Moderada", "Fraca", "Ausente"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(9/667,19/667,21/667,81/667,537/667, 17/666,33/666,37/666,134/666,445/666),
round(data, 3), cex=1, pos=3)
data<-cbind(M = c(6,9,13,53,586)/667, F = c(12,16,25,86,528)/667)
bp<- barplot(height = data,
beside = TRUE,
col=gray.colors(5), ylim=range(c(0,1)),
names.arg = c("Método A", "Método B"), xlab="Masculino", ylab="Proporções amostrais",
legend.text = c("Intolerável", "Intensa", "Moderada", "Fraca", "Ausente"),
args.legend = list(x = "topleft", bty="n", cex=1.2))
abline(h=0)
text(bp, c(6/667,9/667,13/667,53/667,586/667, 12/667,16/667,25/667,86/667,528/667),
round(data, 3), cex=1, pos=3)
dev.copy2eps(file="Figura87g.eps", horizontal=FALSE, pointsize=12, family="Helvetica", colormodel="gray")
# Ajuste e seleção de modelos
require(VGAM)
dados<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/analgesia.txt", h=T)
attach(dados)
dados
mrc1<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo)+factor(met), family=cratio(reverse=F,parallel=T), dados)
mrc2<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo)+factor(met), family=cratio(reverse=F,parallel=F), dados)
coef(mrc1, matrix = TRUE)
coef(mrc2, matrix = TRUE)
TRV<- deviance(mrc1)-deviance(mrc2)
gl <- df.residual(mrc1)-df.residual(mrc2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# Modelo com X1 (sexo)
mrc1<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo), cratio(parallel=T,reverse=F), dados)
mrc2<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo), cratio(parallel=F,reverse=F), dados)
TRV<- deviance(mrc1)-deviance(mrc2)
gl <- df.residual(mrc1)-df.residual(mrc2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# Modelo com X2 (métodos de analgesia)
mrc1<-vglm(cbind(y1,y2,y3,y4,y5)~factor(met), cratio(parallel=T,reverse=F), dados)
mrc2<-vglm(cbind(y1,y2,y3,y4,y5)~factor(met), cratio(parallel=F,reverse=F), dados)
TRV<- deviance(mrc1)-deviance(mrc2)
gl <- df.residual(mrc1)-df.residual(mrc2)
p <- 1-pchisq(TRV,gl); cbind(TRV, gl, p)
# Modelo com X1 e X2 com chances proporcionais de ambas
mrc0<-vglm(cbind(y1,y2,y3,y4,y5)~1, cratio(parallel=T,reverse=F), dados)
mrc0
mrc1<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo),cratio(parallel=T,reverse=F), dados)
mrc1
mrc2<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo)+factor(met), cratio(parallel=T,reverse=F), dados)
mrc2
mrc3<-vglm(cbind(y1,y2,y3,y4,y5)~factor(sexo)+factor(met)+factor(sexo)*factor(met),
cratio(parallel=T,reverse=F), dados)
mrc3
AIC(mrc0)
AIC(mrc1)
AIC(mrc2)
AIC(mrc3)
# Selecionado: modelo com X1 e X2 com chances proporcionais de ambas
summary(mrc2)
coef(mrc2,matrix = TRUE)
mrc2@y
fitted(mrc2)
mrc2@y - fitted(mrc2)
rp<-resid(mrc2, type = "pearson")
Qp<-sum(rp^2)
cbind(Qp,1-pchisq(Qp,10))
QL<-deviance(mrc2)
cbind(QL,1-pchisq(QL,10))
par(mfrow=c(1,2))
rb<- mrc2@y - fitted(mrc2)
plot(1:5,rb[1,], pch=20, ylim=c(-0.05,0.05), xlim=c(1,20), xlab="Índice",
ylab="Probabilidade observada - predita")
points(6:10, rb[2,], pch=20)
points(11:15, rb[3,], pch=20)
points(16:20, rb[4,], pch=20)
abline(h=0, lty=3)
title("(a)")
rp<-resid(mrc2, type = "pearson")
plot(1:4,rp[1,], pch=20, ylim=c(-3,3),xlim=c(1,16), xlab="Índice", ylab="Resíduos de Pearson")
points(5:8, rp[2,], pch=20)
points(9:12, rp[3,], pch=20)
points(13:16,rp[4,], pch=20)
abline(h=0, lty=3)
title("(b)")Capítulo 9 - Regressão logística condicional
9.1 Ensaio clínico multicentros (arquivo de dados: skin.txt) rm(list=ls())
skin<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/skin.txt", h=T)
attach(skin)
require(survival)
model1<-clogit(melhora~trat+sexo+idade+grauini+strata(clinica))
model1
summary(model1)
par(mfrow=c(1,1))
plot(model1$residuals, pch=16)
model2<-clogit(melhora~trat+grauini+strata(clinica))
model2
summary(model2)
plot(model2$residuals, pch=16, ylab="residuos", xlab="i")
# ou equivalentemente
time<-rep(1,158)
model2<-coxph(Surv(time,melhora)~trat+grauini+strata(clinica), method="exact")
summary(model2)
anova(model2) 9.2 Estudo cruzado de dois períodos (arquivo de dados: cross4.txt) cross<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/cross4.txt", h=T)
attach(cross)
require(survival)
model1<-clogit(resp~periodo+drogaA+drogaB+gpidade+resA+resB+strata(obs), data=cross)
model1
summary(model1)
plot(model1$residuals, pch=16)
model2<-clogit(resp~periodo+drogaA+drogaB+gpidade+strata(obs), data=cross)
model2
summary(model2)
plot(model2$residuals, pch=16)
model3<-clogit(resp~periodo+drogaA+drogaB+strata(obs), data=cross)
model3
summary(model3)
# Testando H0: gama_1 = gama_2
model3$var
vardif<-model3$var[2,2] + model3$var[3,3] - 2*(model3$var[2,3])
teste<-((1.408-0.296)/sqrt(vardif))^2
p<-1-pchisq(teste,1)
cbind(teste, p) 9.3 Estudo caso-controle (arquivo de dados: match.txt) match<-read.table("https://docs.ufpr.br/~giolo/LivroADC/Dados/match.txt", h=T)
attach(match)
require(survival)
model1<-clogit(cc~hvb+est+hip+strata(par), data=match)
model1
model2<-clogit(cc~hvb+est+strata(par), data=match)
model2
summary(model2)
plot(model2$residuals, pch=16)