Criando um diretório (pasta) de trabalho, e mudando o workspace do R para este diretório.
getwd() dir.create("~/cursoR") setwd("~/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)
Duplicando e modificando o conjunto de dados em outro objeto
mh <- edit(hills)
str(hills) 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") 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) ## 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() ## 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) 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) head(hills) ## 3 formas de fazxer o gráfico plot(hills$dist, hills$time) with(hills, plot(dist, time)) with(hills, plot(time~dist)) ## 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") 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) ## comentários sobre classes, funções genéricas e metodos methods(plot) ## dividindo a tela gráfica par(mfrow=c(2,2)) plot(lm.h) par(mfrow=c(1,1)) ## 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() 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))) 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) 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)) ## mecanismo 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) ## 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))