CE-709: Estudos Dirigidos em Estatística I - Modelagem estatística

1. Período: terceiro trimestre de 2010, no programa PGMNE (Pós Graduação em Métodos Numéricos em Engenharia)
2. Matrículas e informações: com Maristela, na secretaria do PGMNE
3. Professor Responsável: Paulo Justiniano Ribeiro Jr, (LEG: Laboratório de Estatística e Geoinformação)
4. Horários e Locais:
   • Sala Multiusuário do LEG
   • Horário: Terça, 9:00 às 13:00
   • Atenção: A primeira aula do curso na terça, 14/03/2009.
5. Avaliação: a ser definida

Na oferta do 3o trimestre de 2010 será abordado o conteúdo:
Modelagem estatística

O curso é fortemente baseado no seguinte texto:

Referências adicionais:

ATENÇÃO: arquivos/páginas poderão atualizados durante o curso.

• Programa básico do curso
  1. The R project for Statistical Computing: página do programa R
  2. Um material sobre o uso do R

Veja aqui no histórico das aulas do curso.

Atividades e materiais dos participantes

####################################################### # Extending the Linear Model with R: Generalized Linear, # Mixed Effects and #Nonparametric Regression Models # Julian J.Faraway

####################################################

# CAPÍTULO 1 library(faraway)

# Eleição presidencial dos EUA em 2000, com os dados do # estado da Florida

data(gavote) help(gavote)

# Observando as variáveis e os seis primeiros casos head(gavote)

# Resumo dos dados summary(gavote)

# Considerando os votos nulos de forma relativa # NULOS = (TOTAL - VÁLIDOS)/TOTAL

gavote$undercount ← (gavote$ballots-gavote$votes)/gavote$ballots summary(gavote$undercount)

# Isso é diferente de: # (TOTAL DE NULOS)/ TOTAL sum(gavote$ballots-gavote$votes)/sum(gavote$ballots)

## OBSERVANDO A DISTRIBUIÇÃO DAS VARIÁVEIS

par(mfrow=c(1,2)) # HISTOGRAMA hist(gavote$undercount,main=“NULOS”,xlab=“Porcentagem de Nulos”, ylab=“Frequência”)

# DENSIDADE plot(density(gavote$undercount), main=“NULOS”) rug (gavote$undercount)

## AVALIANDO AS VARIÁVEIS CATEGÓRICAS

par(mfrow=c(1,2))

# PIZZA pie(table(gavote$equip))

# PARETO barplot(sort(table(gavote$equip),decreasing=TRUE),las=2)

## GRÁFICOS EM DUAS DIMENSÕES

# PERCENTUAL DE VOTOS DE GORE gavote$pergore ← gavote$gore/gavote$votes

par(mfrow=c(1,2))

# PERCENTUAL DE VOTOS DE GORE vs PERC DE AA plot(pergore ~ perAA, gavote, xlab=“Proporção de Africanos Americanos”, ylab=“Proporção para Gore”)

#BOX PLOT DE NULOS POR EQUIPAMENTO DE VOTAÇÃO plot(undercount ~ equip, gavote, xlab=“ ”, ylab=“NULOS”, las=3)

## TABELA CRUZADA REGIÃO METROPOL vs RURAL/URBANO xtabs (~ atlanta + rural, gavote)

## CÁLCULO DA CORRELAÇÃO

head(gavote)

# selecionando as colunas com as variáveis nix ← c(3,10,11,12)

#cálculo da correlação entre as variáveis selecionadas cor(gavote[,nix])

## AJUSTANDO MODELOS LINEARES

# NULOS = f(GORE(%), AA(%)) # NULOS = b0 + b1*(GORE(%)) + b2*(AA(%))

lmod ← lm(undercount ~ pergore+perAA, gavote)

# para ML, EMV e EMQ

# APENAS OS COEFICIENTES coef(lmod)

# CALCULANDO OS VALORES PREDITOS predict ← predict(lmod) head(predict) tail(predict)

# CALCULANDO OS RESÍDUOS residuals ← residuals(lmod) head(residuals) tail(residuals)

## QUALIDADE DO AJUSTE

# Em ML, Deviance = SQRes deviance(lmod)

# GL Res df.residual(lmod)

# GL Res = Núm de Obs - Núm de Parms nrow(gavote) - length(coef(lmod))

#VARIÂNCIA DO ERRO sqrt(deviance(lmod)/df.residual(lmod))

# OUTRA FORMA DE PEDIR AS ESTATÍSTICAS names(lmod) summary(lmod) names(summary(lmod))

lmodsummary ← summary(lmod) lmodsummary$sigma lmodsummary$coefficients

## QUALIDADE DO AJUSTE DE FORMA RELATIVA

#R2 e R2ajustado lmodsummary$r.squared

#R2=quadrado da correl entre valores preditos e a resposta cor(predict(lmod), gavote$undercount)^2

# R2 ajustado para o número de parâmetros do modelo lmodsummary$adj.r.squared

##PREDITORES QUALITATIVOS

contr.treatment(5)

#CENTRANDO AS VARIÁVEIS EM RELAÇÃO À MÉDIA gavote$cpergore ← gavote$pergore - mean(gavote$pergore) gavote$cperAA ← gavote$perAA - mean (gavote$perAA)

head(gavote)

#ADICIONANDO VARIÁVEIS QUALITATIVAS AO MODELO lmodi ← lm(undercount ~ cperAA+cpergore*rural+equip,gavote) summary(lmodi)

## TESTANDO HIPÓTESES

#Testa a hipótese de que não há ganho com o modelo maior anova(lmod,lmodi)

#Esta função calcula F para comparar o modelo sem cada # um dos preditores em relação ao modelo completo drop1(lmodi, test=“F”)

## INTERVALOS DE CONFIANÇA

#Intervalo de valores plausíveis para o parâmetro # pode-se julgar o tamanho do efeito do preditor confint(lmodi)

## DIAGNÓSTICOS DO MODELO

plot (lmodi)

#OBSERVANDO OS DOIS CASOS INFLUENTES gavote[cooks.distance(lmodi) > 0.1,]

#OBSERVANDO PONTOS EXTREMOS halfnorm(influence(lmodi)$hat)

#OBSERVANDO OS CASOS INFLUENTES MAIORES QUE 0.3 gavote[influence(lmodi)$hat>0.3,]

#GRÁFICO DOS RESÍDUOS PARCIAIS

termplot(lmodi,partial=T,terms=1)

REGRESSAO ROBUSTA

(MQ funciona bem para erros normais,desempenho insatisfatorio para erros de calda longa. Outliers (varios)….'e #melhor usar uma alternativa robusta funcao rlm) #Metodos de inferencia sao mais dificies de aplicar na estimacao robusta. #O ajuste robusto tem reduzido o efeito dos dois municipios distantes

library(MASS) rlmodi← rlm(undercount ~ cperAA+cpergore*rural+equip,gavote) summary(rlmodi)

#MINIMOS QUADRADOS PONDERADOS### #Os tamanhos dos municipios nesse conjunto de dados variam muito com o numeros de votos em cada municipio. #Uso dos minimos quadrados ponderado- criacao dos pesos proporcionais ao numeros de votos.

wlmodi← lm(undercount ~ cperAA+cpergore*rural+equip,gavote,weights=ballots) #usando o desvio padrao usando a binomial sqrt(0.035*(1-0.035)/881)

#Menor que o erro padrao residual de 0,0233. Logo os efeitos serao menores para os outros municipios.

##TRANSFORMACAO DAS VARIAVEIS

#### #Melhor ajuste Box-Cox- transformacao na variavel resposta ##Considerar o preditor (Africano)no atual modelo## Uso do polinomio ortogonal 'e recomentavel (estavel numericamente e ##facil de selecionar o correto grau)####

plmodi ← lm(undercount ~ poly(cperAA,4)+cpergore*rural+equip, gavote) summary(plmodi)

##Vantagem dos polinomios ortogonais 'e que o coeficientes de baixa ordem nao muda conforme o grau maximo do modelo## #RESIDUOS PARCIAIS#

termplot(plmodi,partial=TRUE,terms=1) #Ajuste estavel usa-se splines

B-SPLINES CUBICOS

#Obs: Complexidade do ajuste B_spline pode ser controlada especificando os graus de liberdade. Usou-se 4. library(splines) blmodi ← lm(undercount ~ cperAA+bs (cpergore,4)+rural+equip, gavote) termplot(blmodi,partial=TRUE,terms=2) #Obs:O ajuste da curva nao 'e muito diferente de uma constante?????

#SELECAO DE VARIAVEIS

#Modelo lineares-previsao -prever medidas de incerteza-remocao de variaveis -pouco para modelo -modelo otimizado-AIC- #AIC =-2 maxima log vero+ 2p onde p= numeros de paramettros

biglm ← lm(undercount ~ (equip+econ+rural+atlanta)^2+ (equip+econ+rural+atlanta)* (perAA+pergore), gavote) ##eliminar termos-minimizar AIC smallm ← step (biglm, trace=F)

##Comparar os modelos AIC### dropl(smallm,test=“F”)??????????

#Interacao do modelo final

## finalm← lm(undercount~equip + econ + perAA + equip:econ + equip:perAA, gavote) summary(finalm) #Obs: Vimos que a rural perAA pode ser descartado. Dando um modelo final#####

##CONCLUSOES-EXPLORAR MODELO FINAL##### pdf ← data.frame(econ=rep(levels(gavote$econ),5),equip=rep(levels(gavote$equip), rep(3,5)),perAA=0.233)

#tabela predict para todas as 15 combinacoes# pp ← predict (finalm, new=pdf) xtabs(round(pp,3) ~ econ + equip,pdf) #obs:menor nos municipios mais ricos e maior nos municipios mais pobres #economia variavel a media# pdf ← data.frame(econ=rep(“middle”,15),equip=rep(levels(gavote$equip), rep(3,5)),perAA=rep(c(.11,0.23,0.35),5)) pp ← predict(finalm, new=pdf)

#fator de tres niveis-baixo,medio e alto propAA ← gl(3,1,15,labels=c(“low”,“medium”,“hiht”)) xtabs(round(pp,3) ~ propAA + equip,pdf) ##Em resumo descobrimos que a situa'cao economica 'e o fator mais claro para determinar a proporcao de votos#

Página aberta para edição pelos participantes do curso.