Não foi possível enviar o arquivo. Será algum problema com as permissões?
Essa é uma revisão anterior do documento!
Ridículas - dicas curtas sobre R
Ridículas é a página do LEG dedicada à fornecer dicas curtas sobre R, e.g. condução de análises, operação com dados e confecção de gráficos. As dicas estão organizadas pelo título, seguido de descrição, palavras-chave e CMR (código mínimo reproduzível). Se você deseja contribuir com a nossa página de Ridículas, envie e-mail para walmes@ufpr.br.
Regressão na análise de variância
Temporariamente sem descrição.
palavras-chave: #regressão, #polinômio, #anova.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados
sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE)
sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar))
str(sorgo)
#
#------------------------------------------------------------------------------------------
# gráficos exploratórios
require(lattice)
xyplot(indice~dose|cultivar, groups=bloco, data=sorgo,
jitter.x=TRUE, type=c("p","l"), layout=c(3,1))
xyplot(indice~dose, groups=cultivar, data=sorgo, jitter.x=TRUE, type=c("p","a"))
#
#------------------------------------------------------------------------------------------
# análise de variância do modelo de fatores
m0 <- aov(indice~bloco+cultivar*ordered(dose), data=sorgo)
summary(m0)
#
#------------------------------------------------------------------------------------------
# checagem
par(mfrow=c(2,2))
plot(m0)
layout(1)
#
#------------------------------------------------------------------------------------------
# desdobrando as somas de quadrados de doses dentro de cultivar
# dicas: forneça para ’by’ o número de níveis de cultivar (3)
# forneça para ’length.out’ os graus de liberdade de dose (6-1)
m1 <- aov(indice~bloco+cultivar/ordered(dose), data=sorgo)
summary(m1)
coef(m1)
summary(m1, split=list("cultivar:ordered(dose)"=list(
"Ag-1002"=seq(1, by=3, length.out=5),
"BR-300"=seq(2, by=3, length.out=5),
"Pioneer-B815"=seq(3, by=3, length.out=5)
)))
#
#------------------------------------------------------------------------------------------
# desdobrando somas de quadrados de cultivar dentro das doses
# dicas: forneça para ’by’ o número de níveis de dose (6)
# forneça para ’length.out’ os graus de liberdade de cultivar (3-1)
m2 <- aov(indice~bloco+ordered(dose)/cultivar, data=sorgo)
coef(m2)
summary(m2, split=list("ordered(dose):cultivar"=list(
"N.0"=seq(1, by=6, length.out=2),
"N.60"=seq(2, by=6, length.out=2),
"N.120"=seq(3, by=6, length.out=2),
"N.180"=seq(4, by=6, length.out=2),
"N.240"=seq(5, by=6, length.out=2),
"N.300"=seq(6, by=6, length.out=2)
)))
#
#------------------------------------------------------------------------------------------
# desdobrando efeitos dos graus polinômio dentro de dose dentro de cultivar
# lof é falta de ajuste (lack of fit)
summary(m1, split=list("cultivar:ordered(dose)"=list(
"Ag-1002.L"=1,
"Ag-1002.Q"=4,
"Ag-1002.C"=7,
"Ag-1002.lof"=c(10,13),
"BR-300.L"=2,
"BR-300.Q"=5,
"BR-300.C"=8,
"BR-300.lof"=c(11,14),
"Pioneer-B815.L"=3,
"Pioneer-B815.Q"=6,
"Pioneer-B815.C"=9,
"Pioneer-B815.lof"=c(12,15)
)))
#
#------------------------------------------------------------------------------------------
# obter as equações de regressão e R^2 para os modelos linear, quadrático e cúbico
# dica: usar contraste tipo soma zero para blocos para se anularem na fórmula
# e remover o intercepto especificando o ’-1’, trocar a ordem dos termos no modelo
# linear (estimativas corretas mas erros padrões e p-valores precisam de correção)
m3 <- aov(indice~-1+cultivar/dose+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m3)
#
#------------------------------------------------------------------------------------------
# quadrático (estimativas corretas mas erros padrões e p-valores precisam de correção)
m4 <- aov(indice~-1+cultivar/(dose+I(dose^2))+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m4)
#
#------------------------------------------------------------------------------------------
# cúbico (estimativas corretas mas erros padrões e p-valores precisam de correção)
m5 <- aov(indice~-1+cultivar/(dose+I(dose^2)+I(dose^3))+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m5)
#
#------------------------------------------------------------------------------------------
# calcular os R^2
sapply(c(linear=1, quadrático=2, cúbico=3),
function(degree){
sapply(levels(sorgo$cultivar),
function(i){
da <- with(subset(sorgo, cultivar==i),
aggregate(indice, list(dose=dose), mean))
summary(lm(x~poly(dose, degree, raw=TRUE), da))$r.squared
})})
#
#------------------------------------------------------------------------------------------
Experimento com dois fatores de efeito aditivo e perda de muitas parcelas
Temporariamente sem descrição.
palavras-chave: #parcela_perdida, #desbalanceamento, #médias_ajustadas, #aditivo.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados
da <- expand.grid(rept=1:5, ep=factor(1:5), tr=factor(1:4))
da$y <- c(58.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
68.4, NA, NA, NA, NA, 258.8, 265.6, NA, NA, NA, NA, NA, 250, NA, 278.8,
268.8, NA, NA, NA, 309.6, NA, NA, NA, NA, NA, NA, NA, NA, NA, 254, 598.8,
NA, NA, NA, NA, 250, 399.6, 260, NA, NA, NA, 288.4, NA, NA, NA, 397.2, NA,
NA, 337.6, NA, 415.2, NA, 450.8, NA, NA, NA, NA, 393.2, NA, NA, NA, NA,
NA, NA, NA, 380.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 634, 417.2,
NA, NA, NA, NA, NA)
#------------------------------------------------------------------------------------------
# ajuste do modelo aditivo com teste F marginal
m0 <- lm(y~ep+tr, data=da)
drop1(m0, test="F")
#------------------------------------------------------------------------------------------
# análise gráfica dos resíduos
par(mfrow=c(2,2))
plot(m0)
layout(1)
#------------------------------------------------------------------------------------------
# estimativas dos efeitos (solução) sob a restrição do R
summary(m0)
#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas dos níveis de tratamento
require(contrast)
lapply(levels(da$tr),
function(i){
contrast(m0, type="average", list(tr=i, ep=levels(da$ep)))
}
)
#------------------------------------------------------------------------------------------
# comparação múltipla de efeitos
require(multcomp)
summary(glht(m0, linfct=mcp(tr="Tukey")))
#------------------------------------------------------------------------------------------
Experimento em parcelas subdivididas com resposta do tipo binomial
Temporariamente sem descrição.
palavras-chave: #binomial, #subdividida, #verossimilhança, #experimento, #glm.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# pacote que contem a função glmmPQL. pode-se usar a lme4::glmer()
require(MASS)
#------------------------------------------------------------------------------------------
# simulandos com a estrutura de um modelo de parcelas subdivididas e resposta binomial
# melhor forma de entender como são gerados os valores observados num experimento
nA <- 4; # número de níveis de um fator A de efeito fixo
nB <- 8; # número de níveis de um fator B de efeito aleatório
nC <- 2; # número de níveis de um fator C de efeito fixo
da <- expand.grid(A=factor(1:nA), B=factor(1:nB), C=factor(1:nC))
da$parcela <- with(da, factor(paste(A,B,sep="-"))) # níveis de parcela aleatório
XA <- model.matrix(~A-1, data=da) # matriz dos efeitos fixos para A
XB <- model.matrix(~B-1, data=da) # matriz dos efeitos aleatórios de B
XC <- model.matrix(~C-1, data=da) # matriz dos efeitos fixos para C
Xp <- model.matrix(~parcela-1, data=da) # matriz dos efeitos aleatórios de parcela
bA <- c(-1,0,1,1) # vetor de efeitos fixos A
bB <- rnorm(nlevels(da$B),0,1) # vetor de efeitos aleatórios B
bC <- c(0,-1) # vetor de efeitos fixos C
bp <- rnorm(nlevels(da$parcela),0,1) # vetor de efeitos aleatórios parcela
da$eta <- XA%*%bA+XB%*%bB+XC%*%bC+Xp%*%bp # vetor de valores preditos
summary(exp(da$eta)/(1+exp(da$eta))) # vetor das probabilidades preditas
da$yobs <- rbinom(nrow(da),
prob=exp(da$eta)/(1+exp(da$eta)),
size=100) # vetor de dados observados
summary(da$yobs)
str(da)
#------------------------------------------------------------------------------------------
# estimação de parâmetros via máxima verossimilhança (usando Penalized Quasi Likelihood)
m0 <- glmmPQL(cbind(yobs, 100-yobs)~A*C, random=~1|B/parcela, data=da, family=binomial)
summary(m0)
#------------------------------------------------------------------------------------------
Duas ou mais funções em um mesmo gráfico
Temporariamente sem descrição.
palavras-chave: #gráfico, #função, #curva.
#------------------------------------------------------------------------------------------ # por Walmes #------------------------------------------------------------------------------------------ # cria os vetores usando as funções abaixo e plota 2 num mesmo gráfico x <- seq(0, 10, by=0.1) y1 <- 4-0.3*x+0.02*x^2 y2 <- 3+exp(-0.9*x) matplot(x, cbind(y1, y2), lty=1, col=1:2, type="l") #------------------------------------------------------------------------------------------
