Ajuste de modelos lineares e mistos no ambiente R |
|
|
09 e 10 de Outubro de 2014 - Piracicaba - SP |
Prof. Dr. Walmes M. Zeviani |
|
Escola Superior de Agricultura “Luiz de Queiroz” - USP |
Lab. de Estatística e Geoinformação - LEG |
| Pós Graduação em Genética e Melhoramento de Plantas | Departamento de Estatística - UFPR |
##-----------------------------------------------------------------------------
## Sequências regulares.
1:5
2*(1:5)-1
3*(1:5)-2
x <- seq(from=1, to=12, by=2); x
x <- seq(from=1, to=12, length.out=5); x
x <- seq(from=1, by=3, length.out=5); x
x <- seq(to=20, by=3, length.out=6); x
##-----------------------------------------------------------------------------
## Repetições regulares.
x <- rep(1, times=5); x
x <- rep(1:3, times=3); x
x <- rep(1:3, each=3); x
##-----------------------------------------------------------------------------
## Grides.
da <- expand.grid(x1=1:3, x2=1:3)
da
da <- expand.grid(trat=c("tratado","controle"),
bloco=c("I","II","III"))
da##-----------------------------------------------------------------------------
## Números uniformes.
x <- runif(10); x
##-----------------------------------------------------------------------------
## Números da distribuição normal padrão.
x <- rnorm(10, mean=0, sd=1)
##-----------------------------------------------------------------------------
## Amostra aleatória de população discreta finita com e sem reposição.
u <- c("Superman","Batman","Aquaman","Flash","Green Lantern")
sample(u, size=3)
u <- c("Iron Man","Thor","Captain America","Hulk","Black Widow","Hawkeye")
sample(u, size=10, replace=TRUE)##-----------------------------------------------------------------------------
## Um vetor com número uniformes.
x <- runif(6)
## As posições do menor para o maior valor.
order(x)
## O vetor ordenado do menor para o maior.
sort(x, decreasing=FALSE)
x[order(x)]
## Os valores em ordem inversa.
rev(x)
## Ordenando um vetor por valores em outro.
u <- c("Iron Man","Thor","Captain America","Hulk","Black Widow","Hawkeye")
u[order(x)]
## O vetor com os elementos aleatorizados.
sample(u)##-----------------------------------------------------------------------------
## Um vetor com dados de precipitação.
precip
x <- precip
sum(x)
mean(x)
median(x)
min(x)
max(x)
range(x)
fivenum(x)
IQR(x)
quantile(x, probs=c(0.05, 0.95))
var(x)
sd(x)
mad(x)
diff(range(x))
100*sd(x)/mean(x)
require(fBasics)## Loading required package: fBasics
## Loading required package: MASS
## Loading required package: methods
## Loading required package: timeDate
## Loading required package: timeSeries
##
## Attaching package: 'fBasics'
##
## The following object is masked from 'package:base':
##
## normskewness(x)
kurtosis(x)
require(EnvStats)## Loading required package: EnvStats
##
## Attaching package: 'EnvStats'
##
## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness
##
## The following object is masked from 'package:MASS':
##
## boxcox
##
## The following object is masked from 'package:stats':
##
## predict.lm## Algumas funções são substituídas.
skewness(x)
kurtosis(x)
cv(x)
summaryStats(x)
summaryFull(x)A família *apply e agregados representa um conjunto de funções destinadas à tarefas por estrato e/ou margem (índice). Tarefas assim foram recentemente batizadas de SAC (split-apply-combine). São essas funções do pacote básico do R que permiter obter médias por tratamento, coeficientes de variação por resposta (por coluna). Isoladamente cada uma delas têm seu papel e juntas resolvem uma série de problemas.
##-----------------------------------------------------------------------------
## Os membros da família *apply.
## De uso simples e frequente.
## apply : on Arrays margins (não é aaply para evitar cacofonia).
## lapply : on Lists itens or vector elements or data.frame columns.
## sapply : lapply that Simplifies when possible.
## tapply : on a ragged array, return as Tabular format.
## De uso menos frequente.
## mapply : lapply over Multiple list or vector arguments.
## eapply : on Environments.
## rapply : Recursive.
## vapply : Vectorized.
## dendrapply : related to dendrogramns.
## kernapply : related kernel.
## Além destas, tem-se outras funções que são os agregados da família.
## by(): divide de acordo com estrato e aplica algo.
## aggragate(): divide por estrato e aplica por coluna.
## sweep(): operações entre matrizes e vetores.
##-----------------------------------------------------------------------------
## tapply.
str(npk)
## Média de yield para cada nível de N.
r <- tapply(npk$yield, npk$N, mean)
r
class(r)
## Média de yield para combinando os níveis de N e P.
r <- tapply(npk$yield, list(npk$N, npk$P), mean)
r
## Melhor usa com with() para simplificar a declaração.
r <- with(npk, tapply(yield, list(N, P), mean))
r
## O mesmo considerando N, P e K.
r <- with(npk, tapply(yield, list(N, P, K), mean))
r
## Nomes na lista geram nomes para as dimensões do array.
r <- with(npk, tapply(yield, list(Nitro=N, Phos=P, Pot=K), mean))
r
##-----------------------------------------------------------------------------
## aggregate.
## A aggragate funciona com uso de formula, além de poder ser usada como a
## tapply(). O resultado é em data.frame.
s <- with(npk,
aggregate(cbind(Y=yield), list(Nitro=N, Phos=P, Pot=K), mean))
s
## Será obtido o mesmo, mas usando uma fórmula para representar o que se
## deseja.
s <- aggregate(yield~N+P+K, data=npk, mean)
s
## Diferente da tapply, a aggregate pode ter mais de uma variáveis
## resposta. Por falta de outra variável resposta, será usando o log de
## yield.
s <- aggregate(cbind(y=yield, log.y=log(yield))~N+P+K, data=npk, mean)
s
##-----------------------------------------------------------------------------
## by.
by(data=npk, INDICES=with(npk, N), FUN=nrow)
r <- with(npk, by(yield, INDICES=N, FUN=mean)); r
str(r)
c(is.list(r), is.array(r))
r <- with(npk, by(yield, INDICES=list(N=N, P=P, K=K), FUN=mean)); r
str(r)
c(is.list(r), is.array(r))
class(r) ## Como é de classe by ele é mostrado de forma diferente.
unclass(r) ## Se a classe é removida, então é mostrado como array comum.Em resumo, as funções tapply, aggregate e by fazem tarefas por estrato. Em outras palavras, separam os valores em um (ou mais) vetores respeitando valores em outro (ou mais) e em seguida aplicam uma função. A diferença é como declarar e o que é retornado.
##-----------------------------------------------------------------------------
## apply.
Titanic
str(Titanic)
is.array(Titanic)
dimnames(Titanic)
sum(Titanic[ ,1, , ]) ## Total de homens.
sum(Titanic[ ,2, , ]) ## Total de mulheres.
apply(Titanic, MARGIN=2, sum) ## Totais das margens para Sex.
apply(Titanic, MARGIN=1, sum) ## Para Class.
apply(Titanic, MARGIN=c(1,2), sum) ## Class e Sex.
apply(Titanic, MARGIN=c(2,4), sum) ## Sex e Survived.
apply(Titanic, MARGIN=c(3,4), sum) ## Age e Survived.
str(HairEyeColor)
dimnames(HairEyeColor)
apply(HairEyeColor, MARGIN=1, sum) ## Por cor de cabelo.
apply(HairEyeColor, MARGIN=2, sum) ## Por cor de olhos.
apply(HairEyeColor, MARGIN=3, sum) ## Por cor de sexo.
##-----------------------------------------------------------------------------
## lapply e sapply.
is.list(rock)
str(rock) ## Todas as colunas tem conteúdo numérico.
lapply(rock, mean) ## Média
lapply(rock, range) ## Extremos.
## Porque trata-se de um data.frame, dá pra usar apply também.
apply(rock, MARGIN=2, mean)
apply(rock, MARGIN=2, range)
sapply(rock, mean) ## Foi possível simplificar para um vetor.
sapply(rock, range) ## Foi possível simplificar para uma matriz.
str(iris)
lapply(iris, is.numeric) ## Quais colunas tem conteúdo numérico?
lapply(iris, class) ## Qual a classe?
## Também se pode usar apply pois iris é um data.frame.
apply(iris, MARGIN=2, is.numeric)
sapply(iris, class) ## Foi possível simplificar.
## Não foi possível simplificar.
sapply(iris, summary)
## Separar as colunas que são númericas e então pedir o summary.
i <- sapply(iris, is.numeric); i
sapply(iris[,i], summary) ## Simplificou.Em resumo, as funções apply, lapply e sapply aplicam funções sobre índices, ou seja, em arranjos as operações são varrendo linhas, colunas, etc e em listas são varrendo os ítens. A diferença de lapply e sapply é que a última tenta simplificar o resultado acomodanto em um tipo de objeto mais simples que uma lista quando possível.
##-----------------------------------------------------------------------------
## Função para calcular a hipotenusa de um triângulo retângulo.
hipo <- function(a, b){ ## Argumentos.
h <- sqrt(a^2+b^2) ## Processo/procedimento.
return(h) ## Retorno.
}
## Explorando o objeto.
class(hipo)
mode(hipo)
args(hipo)
body(hipo)
str(hipo)
## Usando a função.
hipo(3, 4)
hipo(3, 4.6)
hipo(3, 3)
## Note que ela opera vetorialmente.
hipo(a=3, b=c(4, 4.6, 3))
##-----------------------------------------------------------------------------
## Abrindo algumas funções.
body(plot) ## Indica que é uma função genérica.
methods(generic.function="plot") ## As várias 'faces' da plot.
## Abrindo a função.
plot.default
plot.ecdf
## plot.ts
## As com asteriscos podem ser abertas assim.
getAnywhere(plot.histogram)
## getS3method(f="plot", class="histogram")
## Abrindo a função sequência.
methods(generic.function="seq")
getAnywhere(seq.default)
## Abrindo a matrix.
matrix
runif ## Chama uma função escrita em C.
##-----------------------------------------------------------------------------
## Incluindo o argumento ângulo.
hipo2 <- function(a, b, angle){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo2(3, 4, angle=pi/2)
hipo(3, 4)
hipo2(3, 4, angle=pi/2+0.3)
hipo2(3, 4, angle=pi/2-0.3)
##-----------------------------------------------------------------------------
## Ângulo com valor default de 90 graus (pi/2 radianos).
hipo3 <- function(a, b, angle=pi/2){
h2 <- a^2+b^2-2*cos(angle)
return(sqrt(h2))
}
hipo3(3, 4)
hipo3(3, 4, angle=pi/2+0.5)
##-----------------------------------------------------------------------------
## Três implementações.
baskara1 <- function(a, b, c){
x1 <- (-b-sqrt(b^2-4*a*c))/(2*a)
x2 <- (-b+sqrt(b^2-4*a*c))/(2*a)
return(c(x1, x2))
}
baskara2 <- function(a, b, c){
s <- sqrt(b^2-4*a*c)
d <- 2*a
x1 <- (-b-s)/d
x2 <- (-b+s)/d
return(c(x1, x2))
}
baskara3 <- function(a, b, c){
x <- (-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a)
return(x)
}
a <- 2; b <- 1; c <- -3
curve(a*x^2+b*x+c, -3, 3)
abline(h=0, lty=2)
baskara1(a, b, c)
baskara2(a, b, c)
baskara3(a, b, c)
##-----------------------------------------------------------------------------
## Tempo para 50 mil excussões de cada uma.
system.time(replicate(50000, baskara1(a, b, c)))
system.time(replicate(50000, baskara2(a, b, c)))
system.time(replicate(50000, baskara3(a, b, c)))
## As implementação diferentes implicam em custos diferentes.
##-----------------------------------------------------------------------------
## Parábola sem raízes.
a <- 2; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
abline(h=0, lty=2)
baskara1(a, b, c)
##-----------------------------------------------------------------------------
## a==0, então não tem curvatura.
a <- 0; b <- 1; c <- 3
curve(a*x^2+b*x+c, -5, 5)
baskara1(a, b, c)
## Como notificar da origem desses erros?
##-----------------------------------------------------------------------------
## Colocando mensagens de erro para testes feitos nos argumentos.
baskara4 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
delta <- b^2-4*a*c
if (delta<=0) {
stop("Função sem raíz pois `delta` é não positivo.")
}
x <- (-b+c(-1,1)*sqrt(delta))/(2*a)
return(x)
}
baskara4(2, 1, -3)
baskara4(0, 1, -3)
baskara4(2, 1, 3)
##-----------------------------------------------------------------------------
## Função que calcula as raízes da equação e ainda os valores no ponto
## crítico (x, y). O ponto crítico sempre existe se a!=0. As raízes
## podem não existir.
baskara5 <- function(a, b, c){
if (a==0){
stop("O argumento `a` tem que ser diferente de zero.")
}
d <- 2*a
delta <- b^2-4*a*c
if (delta<=0){
warning("Função sem raíz pois `delta` é não positivo.")
r <- c(NA, NA)
} else {
r <- (-b+c(-1,1)*sqrt(delta))/d
}
x <- -b/d
y <- a*x^2+b*x+c
L <- list(x=x, y=y, r=r)
return(L)
}
a <- baskara5(2, 1, -3)
str(a)
curve(2*x^2+1*x-3, -2, 2)
abline(v=a$x, h=a$y, col=2, lty=2)
abline(v=a$r, h=0, col=3, lty=2)
##-----------------------------------------------------------------------------
## Argumentos como vetores, usa posição.
baskara6 <- function(abc){
if (length(abc)==3L){
x <- (-abc[2]+c(-1,1)*sqrt(abc[2]^2-4*abc[1]*abc[3]))/(2*abc[1])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3.")
}
}
baskara6(c(2,1,-3))
baskara6(c(2,1,-3,500))
baskara6(c(2,1))
##-----------------------------------------------------------------------------
## Argumentos como vetores nomeados.
baskara7 <- function(abc){
if (length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- (-abc["b"]+c(-1,1)*
sqrt(abc["b"]^2-4*abc["a"]*abc["c"]))/(2*abc["a"])
return(x)
} else {
stop("O vetor `abc` deve ter comprimento 3 e ser nomeado.")
}
}
baskara7(c(a=2, b=1, c=-3))
baskara7(c(b=1, a=2, c=-3))
baskara7(c(b=1, a=2, m=-3))
##-----------------------------------------------------------------------------
## Argumentos como vetores ou lista nomeados.
baskara8 <- function(abc){
if (is.vector(abc)) abc <- as.list(abc)
if (is.list(abc) & length(abc)==3L & all(names(abc)%in%c("a","b","c"))){
x <- with(abc,
(-b+c(-1,1)*sqrt(b^2-4*a*c))/(2*a))
return(x)
} else {
stop("O objeto `abc` deve ser vetor/lista nomeado de 3 elementos.")
}
}
baskara8(c(a=2, b=1, c=-3))
baskara8(list(a=2, b=1, c=-3))
baskara8(data.frame(a=2, b=1, c=-3))print(sessionInfo(), locale=FALSE)## R version 3.1.1 (2014-07-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
##
## attached base packages:
## [1] methods stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] EnvStats_1.0.2 fBasics_3010.86 timeSeries_3010.97 timeDate_3010.98
## [5] MASS_7.3-34 rmarkdown_0.2.68 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 stabledist_0.6-6
## [6] stringr_0.6.2 tools_3.1.1Sys.time()## [1] "2014-10-07 14:52:37 BRT"