Sessão Inicial -- Fundamentos da Linguagem R

O R é uma linguagem interpretada

4 + 7
sin(pi)
log(10) 
exp(1)

Dados são armazenados em objetos.
O objeto mais simples contém um único elemento que pode ser numérico ou caracter.
O objeto é impresso (exibido) quando digitamos o nome dele.

obj1 <- 3 + log(20)
obj1
obj2 <- "MCT"
obj2
obj3 <- rpois(1, lam=5)
obj3

O R é uma linguagem matricial.
Podemos definir definir vetores, matrizes etc e manipular estes objetos diretamente

v1 <- c(12, 14, 18, 21, 15, 34)
v1
v2 <- 1:6
v2
v3 <- seq(0, 100, length=6)
v3
v4 <- seq(0, 50, by=10)
v4
v5 <- (1:10)*10
v5
v6 <- (-5:5)*10 + 1
v6

Matrizes são vetores com atributos de dimensão.

m1 <- matrix(v1, nrow=3)
m1
m2 <- matrix(1:20, ncol=2)
m2

Há diferentes tipos de armazenagem de objetos

class(obj1)
class(obj2)
class(obj3)
obj3 <- as.integer(obj3)
class(obj3)
object.size(rpois(1000, lam=10))
object.size(as.integer(rpois(1000, lam=10)))

Objetos são armazenados em memória na área de trabalho (workspace) e podem ser gravados/persistidos em um arquivo de formato próprio que por padrão se chama .RData localizado no diretório de trabalho. Pode-se executar comandos do sistema operacional a partir do R.
Criando um diretório (pasta) de trabalho e mudando o workspace do R para este diretório.

getwd()
dir.create("c:\\cursoR")
setwd("c:\\cursoR")
getwd()

Vamos usar o conjunto de dados hills do pacotes MASS

require(MASS)
head(hills)
dim(hills)
colnames(hills)  # nomes das colunas
rownames(hills)
class(hills)
 
str(hills)

Duplicando e modificando o conjunto de dados em outro objeto

mh <- edit(hills)

Manipulando o objeto

dh <- dim(hills)
dh
dh[1]
dh[2]
nrow(hills)
ncol(hills)
 
vars <- colnames(hills)
vars
str(vars)
 
names(mh)
names(mh)[2]
names(mh)[2] <- "altitude"
names(mh)
 
names(mh) <- paste("variavel", 1:3, sep="-")
head(mh)

Criando e manipulando um objeto de dados.

x <- rnorm(10)
class(x)
str(x)
x[2]
x[2:5]
x[x>0]
x[x<0]
 
x <- round(rnorm(10, mean=70, sd=3), dig=2)
x
y <- sample(c("M", "F"), 10, rep=T)
y
#################################
mean(x[y == "M"])
mean(x[y == "F"])
sum(y == "M")
sum(y == "F")
sum(y == "M")/length(y)
sum(y == "F")/length(y)
mean(y == "M")
mean(y == "F")
mean(y != "M")
 
## lei da reciclagem (recycling rule)
x
x + 10
x + c(10, 20)
x + c(10, 20, 30) ## ou .. x + 10*(1:3)
 
x
y
pessoas <- data.frame(peso=x, sexo=y)
pessoas
rm(x,y)
 
head(pessoas)
names(pessoas)
pessoas$peso
pessoas$sexo
 
dim(pessoas)  # dimensão do objeto
pessoas[1:6,]
pessoas[c(3, 5, 8),]
pessoas[pessoas$sexo=="M",]
pessoas[pessoas$sexo=="F",]
 
subset(pessoas, sexo=="M")
subset(pessoas, sexo=="M" | sexo == "F")
subset(pessoas, sexo=="M" & sexo == "F")
 
pessoas[,1]
pessoas$peso
pessoas[,"peso"]
 
with(pessoas, tapply(peso, sexo, mean))
with(pessoas, tapply(peso, sexo, sd))
with(pessoas, tapply(peso, sexo, function(x) 100*sd(x)/mean(x)))
 
## tapply() ~ by() ~ aggregate()

Objetos mais genéricos: listas

## lista
ml <- list(a=1:5, b=matrix(rpois(20, lambda=3), nc=4), c=pessoas)
ml

str(ml)
length(ml)

ml$b

is.list(ml)
is.list(ml$b)
is.matrix(ml$b)
is.data.frame(ml$b)

Um tipo de objeto especial: objetos que não armazenam dados: funções

fc <- function(x){
  100* sd(x)/mean(x)
}
 
fc <- function(x){
  if(!is.numeric(x))
    stop("objeto não numérico")
  100* sd(x)/mean(x)
}
fc(pessoas$peso)
fc(pessoas$sexo)
 
fc1 <- function(x){
  if(is.numeric(x)) res <- 100* sd(x)/mean(x)
  else res <- table(x)
  return(res)
}
 
fc1(pessoas$peso)
fc1(pessoas$sexo)
 
fc2 <- function(x) {
  if(is.numeric(x)) {
    m <- mean(x)
    s <- sd(x)
    cv <- 100*s/m
    ai <- diff(range(x))
    res <- c(média=m, "desvio padrão"=s, CV=cv, amplitude=ai)
  }
  else{
    tb <- table(x)
    moda <- tb[which.max(tb)]
    res <- list(frequencias=tb, moda=moda)
  }
  return(res)
}
 
fc2(pessoas$peso)
fc2(pessoas$sexo)
 
lapply(pessoas, fc2)
 
df.res <- lapply(pessoas, fc2)
df.res
is.list(df.res)
str(df.res)
 
summary(pessoas)
sum.res <- summary(pessoas)
sum.res
class(sum.res)
str(sum.res)
 
## redirecionando saida para um arquivo texto
sink("saidas.txt")
summary(pessoas)
lapply(pessoas, fc2)
pessoas
sink() ## voltando para o dispositivo de saida padrão (tela)
 
summary(pessoas)

Visualizando e analisando dados

head(hills)
 
## 3 formas de fazer o gráfico
plot(hills$dist, hills$time)
with(hills, plot(dist, time))
with(hills, plot(time~dist))
 
## alguma cosmética
with(hills, plot(time~dist, pch="*", cex=2.5))
with(hills, plot(time~dist, pch=5))
 
## descobrindo os símbolos
plot(1:25, 1:25, pch=1:25)
## descobrindo as cores
plot(1:8, 1:8, col=1:8, pch=19, cex=2)
plot(1:20, 1:20, col=terrain.colors(20), pch=19, cex=2)
plot(1:10, 1:10, col=heat.colors(10), pch=19, cex=2)
plot(1:13, 1:13, col=gray(seq(0,1,len=13)), pch=19, cex=2)
par(bg="yellow")
plot(1:13, 1:13, col=gray(seq(1,0,len=13)), pch=19, cex=2)
colors()
par(bg="yellow4")
plot(1:13, 1:13, col=gray(seq(1,0,len=13)), pch=19, cex=2)
par(bg="white")

Um breve exemplo de operações matriciais - coeficientes de regressao(modelo) linear

require(MASS)
data(hills)
names(hills)
X <- cbind(1, hills$dist)
y <- hills$time
## Há formas e formas de fazer operações...
## uma forma não muito adequada:
solve(t(X) %*% X) %*% t(X) %*% y
## outra melhor:
solve(crossprod(X), crossprod(X, y)) 

Há funções que já implementam modelos estatísticos usuais.

with(hills, plot(time~dist))
 
lm.h <- lm(time ~ dist, data=hills)
lm.h
names(lm.h)
is.list(lm.h)
 
lm.h$coeffi
coef(lm.h)
lm.h$res
resid(lm.h)
lm.h$fitted
fitted(lm.h)
 
plot(resid(lm.h) ~ fitted(lm.h))
 
anova(lm.h)
summary(lm.h)
 
plot(lm.h)

Uma pausa para comentários sobre classes, funções genéricas e métodos

methods(plot)

Uma breve introdução sobre recursos gráficos

## dividindo a tela gráfica
par(mfrow=c(2,2))
plot(lm.h)
par(mfrow=c(1,1))
 
## parâmetros gráficos
par()

Exportando gráficos: alguns dispositivos gráficos

jpeg("diag.jpg")
par(mfrow=c(2,2))
plot(lm.h)
dev.off()
 
pdf("diag.pdf")
par(mfrow=c(2,2))
plot(lm.h)
dev.off()
 
postscript("diag.eps")
par(mfrow=c(2,2))
plot(lm.h)
dev.off()

Continuando…

with(hills, plot(time ~ dist))
 
lm.h <- lm(time~dist,data=hills)
 
## linhas num gráfico
abline(h=100)
abline(h=c(50,100,150), lty=2)
abline(v=10, lty=3, col=2)
 
with(hills, plot(time ~ dist, ylim=c(0, 250)))
abline(coef(lm.h))
 
with(hills,(segments(x0=dist, y0=fitted(lm.h),
                    x1=dist, y1=time, lty=2)))

A idéia é portanto acrescentar elementos no gráfico.
Alguns exemplos de anotação em gráficos.

title("Regressão entre tempo e distâncias")
text(2, 250, substitute(hat(beta)[1] == b1,
                        list(b1=round(coef(lm.h)[2], dig=2))),
     pos=4, cex=1.5)

Um pouco mais sobre modelos (com e sem intercepto)

lm0.h <- lm(time ~ dist-1, data=hills)
lm0.h
 
with(hills, plot(time ~ dist, ylim=c(0, 250), xlab="distância", ylab="tempo"))
abline(lm.h)
abline(lm0.h, lty=2, col=4)
legend("topleft", c("2 parâmetros","1 parâmetro"), lty=1:2, col=c(1,4))

Mas como saber ???
Alguns mecanismos de ajuda

help(plot)
## help no navegador
help.start(browser="firefox") ## veja no navegador acesso a documentação!!!!
help(plot)
## procura no seu computador
help.search("cluster")
# procura no site do R
RSiteSearch("cluster")
## procura um objeto
find("mean")
find("glm")
## pedaco de palavra (nome do objeto)
apropos("mean")

De volta a regressão…
Predizendo (valores preditos) usado objetos de ajuste de modelos.

with(hills, plot(time~dist, ylim=c(0,250)))
pred.df <- data.frame(dist=seq(0,30,length=100))
pred.h <- predict(lm.h, newdata=pred.df)
pred.h
lines(pred.df$dist, pred.h)
 
## agora acrescentando intervalo de confiança...
predc.h <- predict(lm.h, newdata=pred.df, int="c")
dim(predc.h)
head(predc.h)
matlines(pred.df$dist,predc.h, lty=2, col=1)
## ... e o intervalo de predição...
predp.h <- predict(lm.h, newdata=pred.df, int="p")
matlines(pred.df$dist,predp.h, lty=3, col=1)
 
legend("topleft", c("modelo ajustado", "intervalo de confiança", "intervalo de predição"), lty=1:3)

Análises retirando o ponto mais influente

args(lm)
 
lmo.h <- lm(time ~ dist, data=hills, subset=(dist<25))
summary(lmo.h)
predo.h <- predict(lmo.h, newdata=pred.df)
 
with(hills, plot(time~dist, ylim=c(0,250)))
lines(pred.df$dist, pred.h)
lines(pred.df$dist, predo.h, col=2)
names(hills)
points(hills[hills$dist>25,c(1,3)], pch=19, col=2)

Mostrando os resultados dos modelos com e sem o ponto mais atípico em dois gráficos separados.

par(mfrow=c(2,2))
with(hills, plot(time~dist, ylim=c(0,250)))
lines(pred.df$dist, pred.h)
matlines(pred.df$dist,predc.h, lty=2, col=1)
matlines(pred.df$dist,predp.h, lty=3, col=1)
 
with(hills, plot(time~dist, ylim=c(0,250)))
points(hills[hills$dist>25,c(1,3)], pch=19, col=2)
lines(pred.df$dist, predict(lmo.h, new=pred.df), col=2)
matlines(pred.df$dist,predict(lmo.h, new=pred.df, int="c"), lty=2, col=1)
matlines(pred.df$dist,predict(lmo.h, new=pred.df, int="p"), lty=3, col=1)
 
hist(resid(lm.h))
hist(resid(lmo.h))

Regressão múltipla, fórmulas, transformação e comparação e seleção de modelos.

##
head(hills)
 
lm2.h <- lm(time ~ dist + climb, data=hills)
anova(lm2.h)
summary(lm2.h)
par(mfrow=c(2,2))
plot(lm2.h)
 
shapiro.test(resid(lm2.h))
 
par(mfrow=c(1,1))
boxcox(time ~dist+climb, data=hills)
 
 
lm2r.h <- lm(sqrt(time) ~ dist + climb, data=hills)
coef(lm2r.h)
## ou... aproveitando o modelo previamente definido...
lm2r.h <- update(lm2.h, sqrt(time) ~ ., data=hills)
coef(lm2r.h)
 
summary(lm2r.h)
 
## modelo retirando variavel climb
lm3r.h <- update(lm2r.h, . ~ . - climb)
coef(lm3r.h)
 
anova(lm3r.h, lm2r.h)
 
stepAIC(lm(time ~ dist*climb, data=hills))