SAUDAÇÕES!

Seja bem vindo à página do professor Pedro Albuquerque. Para saber mais sobre meu currículo, disciplinas ministradas e interesses de pesquisa, navegue no menu disponível no topo da página.

sexta-feira, 11 de dezembro de 2015

Cálculo do tamanho de amostras: Equações Estruturais.


Já vimos em dois posts como calcular o Tamanho Amostral necessário sobre a abordagem de populações finitas para Proporções e Médias.

Existe em estatística, dois paradigmas acerca do cálculo do tamanho amostral:

  • Abordagem de Populações Finitas:
  • Nessa proposta, o objetivo é construir uma amostra representativa de uma população finita de modo a manter "mais ou menos" as mesmas características da população alvo para o conjunto de variáveis de interesse.
  • Abordagem de Populações Infinitas (modelo de superpopulação):
  • Já na abordagem de Populações Infinitas o interesse está em obter uma amostra representativa de uma população teoricamente infinita, para que um determinado modelo (Regressão, Correlação, Equações Estruturais, etc.) possa ser aplicado a determinados níveis de Erro Tipo I e Erro Tipo II.

Nesse sentido, o objetivo desse post é apresentar como o tamanho amostral deve ser calculado quando o interesse é a realização de um modelo de Equações Estruturais, o qual engloba modelos de regressão e análise fatorial confirmatória, por exemplo.

Para isso, utilizaremos como base dois textos: MacCallum, Browne e Sugawara (1996) - Power Analysis and Determination of Sample Size for Covariance Structure Modeling e Rigdon(1994) - Calculating degrees of freedom for a structural equation model. Especificamente, MacCallum, Browne e Sugawara (1996) considera que o interesse do analista que executa um modelo de Equações Estruturais é avaliar o adequamento global do modelo segundo alguma medida de ajuste, em especial o RMSEA - The Root Mean Square Error of Approximation. Em outras palavras o interesse é testar algo como:

$$
\begin{cases}
H_{0}:RMSEA>\kappa_{1}\\
H_{a}:RMSEA\leq\kappa_{2}\\
\end{cases}
$$

onde $\kappa_{1}$ e $\kappa_{2}$ são níveis considerados para o RMSEA. MacCallum, Browne e Sugawara (1996) sugere os seguintes níveis:

  • $RMSEA \geq 0.10$ (Ajuste ruim).
  • $RMSEA \leq 0.10$ (Ajuste medíocre)
  • $RMSEA \leq 0.08$ (Ajuste aceitável).
  • $RMSEA \leq 0.05 $ (Ajuste adequado)
  • $RMSEA = 0.00 $ (Ajuste exato)

Assim, para se calcular o tamanho amostral desejado para um modelo de Equações Estruturais são necessários alguns insumos: níveis máximos admitidos para os Erros Tipo I e Tipo II (usualmente $\alpha=0.05$ e $\beta=0.2$, respectivamente), limites considerados para o teste de hipótese do RMSEA (usualmente, $\kappa_{1}=0.08$ e $\kappa_{2}=0.05$) e o número de graus de liberdade do modelo.

O cálculo do tamanho amostral necessário necessita de algum conhecimento prévio teórico sobre os Modelos de Equações Estruturais, uma boa referência é o livro Bollen (2014) - Structural equations with latent variables.

Considere o seguinte modelo proposto por Rigdon(1994):


O cálculo do número de graus de liberdade desse modelo é a diferença entre o número de informações manifestas não redundantes disponíveis (isto é, número de variâncias e covariâncias possíveis de serem calculadas diretamente com base nos dados $\frac{6\times(6+1)}{2}=21$) e o número de parâmetros livres a ser estimados ($\boldsymbol\theta=(\lambda_{1},\dots,\lambda_{6},\phi_{21},\psi_{1},\dots,\psi_{6})^{T})$ onde $\lambda_{1},\dots,\lambda_{6}$ representam as cargas fatoriais, $\phi_{21}$ a covariância entre as variáveis latentes e $\psi_{1},\dots,\psi_{6}$ as variâncias associadas aos termos de erro $\delta_{1},\dots,\delta_{6}$ tal que para o modelo apresentado temos $dim(\boldsymbol\theta)=13$ e portanto o número de graus de liberdade é igual a $21-13=8$).

Portanto, para o teste de hipótese na forma:

$$
\begin{cases}
H_{0}:RMSEA>0.08\\
H_{a}:RMSEA\leq 0.05\\
\end{cases}
$$

No R basta fazer:

#Chama a biblioteca semTools
library(semTools)
#Define o erro do tipo 1:
erro1<-0.05
#Define o erro do tipo 2:
erro2<-0.20
#Define o número de graus de liberdade
gl<-8
#Define o limite para a hipótese nula do RMSEA
k1REMSEA<-0.08
#Define o limite para a hipótese alternativa do RMSEA
k2REMSEA<-0.05
#Calcula o tamanho amostral
findRMSEAsamplesize(rmsea=k1REMSEA,
                    rmseaA=k2REMSEA, 
                    df=gl, power=(1-erro2), alpha=erro1)
Após a execução do código anterior, o pacote semTools fornece um tamanho amostral mínimo igual a 961 observações. É importante ressaltar que cada modelo pode possuir graus de liberdade diferentes, e portanto, o tamanho amostral dependerá da estrutura gráfica considerada para o modelo de Equações Estruturais.

segunda-feira, 16 de novembro de 2015

All your Bayes are belong to us!


A campanha do Facebook para fortalecer sua equipe de pesquisadores em Aprendizado de Máquina obteve uma grande vitória.

Dia 22/11/2014 o laboratório de Inteligência Artificial do Facebook postou no Facebook, seu mais novo colaborador: Vladimir Vapnik. Ele foi responsável pelo desenvolvimento do que hoje é conhecido como Máquina de Suporte Vetorial.

Já a piadinha geek com o termo "All your Bayes are belong to us!" refere-se a um jogo antigo de Mega-Drive (Zero Wing) que teve uma tradução "quebrada" do japonês para o inglês na forma "All your base are belong to us"!, o base trocado por bayes é uma provocação, sugerindo a superioridade do SVM (Suport Vector Machine) sobre as abordagens bayesianas tradicionais.


Note ainda que no topo do quadro há a seguinte desigualdade:


Esta fórmula representa o limite probabilístico (com probabilidade $1-\eta$) do risco esperado de um classificador sobre a abordagem de Empirical Risk Minimization, base da construção das Máquinas de Suporte Vetorial. Em termos simples, essa desigualdade relaciona o erro populacional do classificador com o erro amostral (empírico) em função de um conjunto de treinamento (amostra) de tamanho $l$ e em que a cardinalidade do conjunto de funções de perda é $N$.

Bem, a foto é basicamente uma declaração sobre a superioridade da teoria da aprendizagem estatística de Vapnik sobre a alternativa bayesiana. Em poucas palavras, o paradigma bayesiano começa com a proposição de uma distribuição a priori ao longo de um conjunto de hipóteses (nossas crenças) e então é atualizada por meio dos dados (distribuição a posteriori). Nesse paradigma, a regra de decisão ótima para um classificador é baseado na distribuição a posteriori.

Já no paradigma de Empirical Risk Minimization não há necessidade de pressupostos a cerca da distribuição de probabilidade dos dados. Isto é motivado pelo fato de que a estimativa de densidade é um problema mal-posto (ill-posed), e, portanto, desejamos evitar esse passo intermediário. O objetivo é minimizar diretamente a probabilidade de se fazer uma má decisão no futuro.

quinta-feira, 15 de outubro de 2015

Automação de tarefas usando rSelenium.


A automação de tarefas é uma maneira simples de evitar trabalhos repetitivos, em especial, desejamos coletar dados de sites e por algum motivo, cada observação (ou grupo de observações) precisa ser "manualmente" solicitada.

Essa solicitação pode ser exigida por meio de captchas ou outro mecanismo, impedindo assim muitas consultas ao servidor alvo e garantindo estabilidade ao serviço.

Entretanto, algumas ferramentas auxiliam na construção de bots para a coleta dessas informações, processo esse conhecido como Web scraping. Dentre as ferramentas mais utilizadas, destaco o iMacro e rSelenium (baseado na solução Selenium).

Nesse post trabalharemos com o driver do navegador Chrome, para isso, o primeiro passo é obter esse driver:

#Limpa o workspace
rm(list=ls())

#Define o  Working Directory
setwd("C:/Users/RSelenium")

#Instalação do RSelenium
library(RSelenium)

#Confere se o RSelenium está instalado
checkForServer()

#Baixa o driver: http://chromedriver.storage.googleapis.com/index.html
bool<-FALSE
bool<-file.exists("chromedriver.exe")
if(!bool)
{
  #Baixa o arquivo zip
  down<-"http://chromedriver.storage.googleapis.com/2.13/chromedriver_win32.zip"
  download.file(down,paste(getwd(),"/chromedriver.zip",sep=""))
  
  #Extrai o arquivo
  unzip(paste(getwd(),"/chromedriver.zip",sep=""))
  
  #Deleta o arquivo
  file.remove((paste(getwd(),"/chromedriver.zip",sep="")))
}
O código anterior obtêm o driver e armazena na pasta escolhida, em seguida, iniciamos o navegador:
#Inicia o servidor
startServer(args = c(paste("-Dwebdriver.chrome.driver=",getwd(),"/chromedriver.exe -Dwebdriver.chrome.args='--disable-logging'",sep="")), log = FALSE, invisible = FALSE)
remDr <- remoteDriver(browserName = "chrome")

#Abre o navegador
remDr$open()

#Maximiza a janela
remDr$maxWindowSize()
Nas etapas seguintes é possível navegar por páginas específicas e coletar as informações de interesse, por exemplo, um site muito útil para análise fundamentalista é o sítio http://fundamentus.com.br/. Se desejarmos obter os dados da PETR4, basta acessar:
#Vai para a pagina de interesse
site<-"http://fundamentus.com.br/balancos.php?papel=PETR4"
remDr$navigate(site)
Note que para obter os dados é necessário preencher o valor do captcha, para isso, podemos usar algum sistema de Optical character recognition (OCR) para reconhecer as letras indicadas no captcha. É necessário obter a imagem do captcha, uma solução é fazer um printscreen do site e editar a imagem:
#Faz um printscreen do site
library(base64enc)
img<-remDr$screenshot(display = FALSE, useViewer = TRUE, file = NULL)
writeBin(base64Decode(img, "raw"), 'site.png')
Uma vez obtido o valor do captcha, podemos passar o resultado para o site da seguinte forma:
#Exemplo de captcha solucionado
txt<-"TOWIT"
  
#Encontra o objeto da caixa de texto
webElem <- remDr$findElement(using = "name", "codigo_captcha")

#Manda o resultado do captcha
webElem$sendKeysToElement(list(txt))

#Encontra o objeto do botão de submit
webElem <- remDr$findElement(using = "name", "submit")

#Clica no botão
webElem$clickElement()

#Fecha as conexões
remDr$close()
remDr$closeServer()
Note que os comandos das linhas 5 e 11 precisam ser preenchidos com base nas características do objeto de interesse na página alvo, isso pode ser feito averiguando o código fonte da página por meio de qualquer navegador. Quanto a solução do captcha existem ferramentas de Optical character recognition (OCR) e também de serviços como o http://www.deathbycaptcha.com/ que realizam esse processo.

terça-feira, 15 de setembro de 2015

Credit Score usando o R.


Credit Scoring é definido como sendo um modelo estatístico/econométrico o qual atribui uma medida de risco aos clientes de uma instituição financeira ou aos futuros clientes.

Usualmente, a avaliação dos clientes é realizada por meio de um Credit Scorecard o qual é um modelo estatístico para avaliação do risco estruturado de maneira a facilitar a tomada de decisão quanto a liberação do crédito ou não.

O uso de credit scorecard é muito popular principalmente para as organizações que lidam empréstimos como bancos.

Dentre as vantagens da utilização de um credit scorecard podemos listar:

  1. Credit Scorecard é implementado facilmente e pode ser monitorado ao longo do tempo.
  2. Pessoas sem o conhecimento técnico em estatística ou econometria podem utilizar facilmente o credit scorecard para tomar decisões.

As principais questões em Credit Score são:

  1. Quem receberá o crédito ?
  2. Quanto deverá ser esse crédito ?
  3. Quais as estratégias para a distribuição e cobrança do crédito ?

Em outras palavras, Credit Score é um conjunto de modelos de decisão e técnicas estatísticas subjacentes que auxiliam os credores na tomada de decisão quanto a concessão de crédito ao consumidor.

Estas técnicas informam quem receberá o crédito, quanto crédito deverá ser fornecido e quais estratégias operacionais melhorarão a rentabilidade dos devedores para os credores (Thomas, Eldelman e Crook, 2002).

Credit Score usando o R.

Para demonstrar como o Credit Score pode ser formulado usando o R, iremos trabalhar com os dados German.csv o qual representa um conjunto de dados de crédito para uma instituição financeira Alemã.

Os dados são compostos por 300 empréstimos "ruins" (por exemplo, ausência de pagamento ou atraso) e 700 empréstimos "bons" (por exemplo, pagamentos sem atraso). O objetivo é fornecer insumos para a tomada de decisão quanto aos futuros empréstimos com base nos padrões anteriormente observados.

É comum em Credit Score classificar as contas "ruins" como aquelas contas que em algum período de tempo apresentaram inadimplência por 60 dias ou mais (em empréstimos hipotecários, 90 dias ou mais é por vezes utilizado).

Vamos importar os dados para o R:

#Limpa o Workspace
rm(list=ls())

#Importa os dados German.csv 
dados.df<-read.csv("https://dl.dropboxusercontent.com/u/36068691/Blog/german.csv")

#Apresenta as variáveis do DataFrame
names(dados.df)

#Apresenta a estrutura do DataFrame
str(dados.df)
No R as variáveis categóricas ou binárias são usualmente tratadas como fatores enquanto as variáveis contínuas ou discretas são tratadas como valores numéricos. Nesse caso, algumas conversões são necessárias:
#Transforma em fatores as variáveis categóricas e "dummies"
dados.df[,"CHK_ACCT"]     <-as.factor(dados.df[,"CHK_ACCT"])
dados.df[,"HISTORY"]      <-as.factor(dados.df[,"HISTORY"])
dados.df[,"NEW_CAR"]      <-as.factor(dados.df[,"NEW_CAR"])
dados.df[,"USED_CAR"]     <-as.factor(dados.df[,"USED_CAR"])
dados.df[,"FURNITURE"]    <-as.factor(dados.df[,"FURNITURE"])
dados.df[,"RADIO.TV"]     <-as.factor(dados.df[,"RADIO.TV"])
dados.df[,"EDUCATION"]    <-as.factor(dados.df[,"EDUCATION"])
dados.df[,"RETRAINING"]   <-as.factor(dados.df[,"RETRAINING"])
dados.df[,"SAV_ACCT"]     <-as.factor(dados.df[,"SAV_ACCT"])
dados.df[,"EMPLOYMENT"]   <-as.factor(dados.df[,"EMPLOYMENT"])
dados.df[,"MALE_DIV"]     <-as.factor(dados.df[,"MALE_DIV"])
dados.df[,"MALE_SINGLE"]  <-as.factor(dados.df[,"MALE_SINGLE"])
dados.df[,"MALE_MAR"]     <-as.factor(dados.df[,"MALE_MAR"])
dados.df[,"CO.APPLICANT"] <-as.factor(dados.df[,"CO.APPLICANT"])
dados.df[,"GUARANTOR"]    <-as.factor(dados.df[,"GUARANTOR"])
dados.df[,"TIME_RES"]     <-as.factor(dados.df[,"TIME_RES"])
dados.df[,"REAL_ESTATE"]  <-as.factor(dados.df[,"REAL_ESTATE"])
dados.df[,"PROP_NONE"]    <-as.factor(dados.df[,"PROP_NONE"])
dados.df[,"OTHER_INSTALL"]<-as.factor(dados.df[,"OTHER_INSTALL"])
dados.df[,"RENT"]         <-as.factor(dados.df[,"RENT"])
dados.df[,"OWN_RES"]      <-as.factor(dados.df[,"OWN_RES"])
dados.df[,"NUM_CREDITS"]  <-as.factor(dados.df[,"NUM_CREDITS"])
dados.df[,"JOB"]          <-as.factor(dados.df[,"JOB"])
dados.df[,"NUM_DEPEND"]   <-as.factor(dados.df[,"NUM_DEPEND"])
dados.df[,"TELEPHONE"]    <-as.factor(dados.df[,"TELEPHONE"])
dados.df[,"FOREIGN"]      <-as.factor(dados.df[,"FOREIGN"])

#Variável dependente
dados.df[,"RESPONSE"]     <-as.factor(dados.df[,"RESPONSE"])

#Transforma em numeric
dados.df[,"AMOUNT"]       <-as.numeric(dados.df[,"AMOUNT"])
dados.df[,"INSTALL_RATE"] <-as.numeric(dados.df[,"INSTALL_RATE"])
dados.df[,"AGE"]          <-as.numeric(dados.df[,"AGE"])
dados.df[,"NUM_DEPEND"]   <-as.numeric(dados.df[,"NUM_DEPEND"])
dados.df[,"DURATION"]     <-as.numeric(dados.df[,"DURATION"])
O próximo passo é separar os dados em dois grupos:
  1. Dados para estimação. (Treinamento)
  2. Dados para teste. (Validação)
Uma sugestão é separar a base de dados (aleatoriamente) da seguinte forma: 60% das observações deverão compor a base de treinamento e 40% a base de validação:
#Índices obtidos após a aleatorização
ordena <- sort(sample(nrow(dados.df), nrow(dados.df)*.6))

#Dados para o treinamento
treinamento<-dados.df[ordena,]

#Dados para a validação
validacao<-dados.df[-ordena,]
A ideia é construir o(s) modelo(s) de Credit Scoring com o DataFrame "treinamento" e em seguida avaliar o ajuste com o DataFrame "validacao". Uma das formas mais simples para modelar os dados de crédito é por meio da regressão logística:
#Índices obtidos após a aleatorização
ordena = sort(sample(nrow(dados.df), nrow(dados.df)*.6))

#Dados para o treinamento
treinamento<-dados.df[ordena,]

#Dados para a validação
validacao<-dados.df[-ordena,]
Como há muitas possíveis variáveis para o modelo, podemos proceder com a abordagem de Stepwise para selecionar o modelo com a "melhor" combinação de variáveis explicativas:
#Regressão Logística
modelo.completo <- glm(RESPONSE ~ . ,family=binomial,data=treinamento)

#Abordagem Stepwise para seleção de variáveis
stepwise <- step(modelo.completo,direction="both") 
Após algumas iterações, observa-se que o conjunto de variáveis com o menor valor para o Critério de Informação de Akaike é:
#Modelo com as variáveis indicadas pelo Stepwise
stepwise <- glm(RESPONSE ~  JOB+NUM_CREDITS+EMPLOYMENT+RETRAINING+NEW_CAR+TELEPHONE+MALE_DIV+
                  FURNITURE+PROP_NONE+MALE_MAR+RENT+NUM_DEPEND+REAL_ESTATE+EDUCATION+FOREIGN+
                  TIME_RES, family=binomial,data=treinamento)

#Resume os resultados do modelo
summary(stepwise)
Percebemos que nem todas as variáveis são significantes:
Uma medida interessante para interpretar o modelo é a medida de Razão de chances (Odds Ratio):
#Calcula a razão de chances
exp(cbind(OR = coef(stepwise), confint(stepwise)))
Obtemos como resultados as seguintes Razões de chance:
E como dito anteriormente, a interpretação é bem interessante. Veja por exemplo a variável NUM_DEPEND, nesse caso, para cada dependente a mais que um proponente possui isso aumenta a sua chance de ser considerado inadimplente em aproximadamente 12%. Finalmente, vamos testar a qualidade do modelo aplicando o modelo estimado na base de validação para termos uma ideia do grau de acerto desse modelo:
#Faz a previsão para a base de validação (probabilidade)
predito<-predict(stepwise,validacao,type="response")

#Escolhe quem vai ser "1" e quem vai ser "0"
predito<-ifelse(predito>=0.8,1,0)
  
#Compara os resultados
table(predito,validacao$RESPONSE)
Obtemos assim a seguinte matriz de confusão:

Logo a nossa taxa de acerto (acurácia) nesse modelo é dada por:

$
Tx.Acerto=\frac{75+166}{400}=\frac{241}{400} \approx 60\%
$

Podemos melhorar a taxa de acerto refinando o modelo por meio da exclusão de variáveis não significantes e pela inclusão de componentes estatisticamente significantes.

sexta-feira, 14 de agosto de 2015

Escolha de métodos na Análise Estatística.

Escolher o principal método de análise estatística é sempre um desafio, nesse sentido, alguns fluxogramas podem auxiliar nesse processo. No caso de métodos tradicionais de Inferência Estatística uma boa proposta é:


No caso de Métodos Multivariados ou Teoria de Aprendizagem Estatística uma proposta gráfica que o analista pode seguir é:


Gráficos retirados das seguintes fontes: http://abacus.bates.edu/~ganderso/biology/resources/ e http://scikit-learn.org/stable/tutorial/machine_learning_map/

quarta-feira, 15 de julho de 2015

Revolution R - Parte 3.


Vimos anteriormente no post Revolution R - Parte 2 como manipular bases de dados importadas para o ambiente Revolution R. Nesse post veremos como importar dados brutos das pesquisas do IBGE. Considere por exemplo, os dados da Pesquisa Nacional por Amostra de Domicílios - PNAD 2011.

Após baixar os microdados podemos criar o dataset no formato nativo *.XDF do Revolution R da seguinte forma:

#Define o working directory
setwd("C:/PNAD2011/Dados")

#Define as variáveis e seu tipo
colList <- list(
"V0101"=list(type="factor", start=1, width=4,description="Ano de referência"),
"UF"=list(type="factor", start=5, width=2,description="Unidade da Federação"),
"V0102"=list(type="factor", start=5, width=8,description="Número de controle"),
"V0103"=list(type="factor", start=13, width=3,description="Número de série"),
"V0104"=list(type="factor", start=16, width=2,description="Tipo de entrevista"),
"V0105"=list(type="factor", start=18, width=2,description="Total de moradores"),
"V0106"=list(type="factor", start=20, width=2,description="Total de moradores de 10 anos ou mais"),
"V0201"=list(type="factor", start=22, width=1,description="Espécie do domicílio"),
"V0202"=list(type="factor", start=23, width=1,description="Tipo do domicílio"),
"V0203"=list(type="factor", start=24, width=1,description="Material predominante na construção das paredes externas do prédio"),
"V0204"=list(type="factor", start=25, width=1,description="Material predominante na cobertura (telhado) do domicílio"),
"V0205"=list(type="integer", start=26, width=2,description="Número de cômodos do domicílio"),
"V0206"=list(type="integer", start=28, width=2,description="Número de cômodos servindo de dormitório"),
"V0207"=list(type="factor", start=30, width=1,description="Condição de ocupação do domicílio"),
"V0208"=list(type="numeric", start=31, width=12,description="Aluguel mensal pago no mês de referência"),
"V0209"=list(type="numeric", start=43, width=12,description="Prestação mensal paga no mês de referência"),
"V0210"=list(type="factor", start=55, width=1,description="Terreno onde está localizado o domicílio é próprio"),
"V0211"=list(type="factor", start=56, width=1,description="Tem água canalizada em pelo menos um cômodo do domicílio"),
"V0212"=list(type="factor", start=57, width=1,description="Proveniência da água canalizada utilizada no domicílio"),
"V0213"=list(type="factor", start=58, width=1,description="Água utilizada no domicílio é canalizada de rede geral de distribuição para a propriedade"),
"V0214"=list(type="factor", start=59, width=1,description="Água utilizada no domicílio é de poço ou nascente localizado na propriedade"),
"V0215"=list(type="factor", start=60, width=1,description="Tem banheiro ou sanitário no domicílio ou na propriedade"),
"V0216"=list(type="factor", start=61, width=1,description="Uso do banheiro ou sanitário"),
"V2016"=list(type="integer", start=62, width=2,description="Número de banheiros ou sanitários"),
"V0217"=list(type="factor", start=64, width=1,description="Forma de escoadouro do banheiro ou sanitário"),
"V0218"=list(type="factor", start=65, width=1,description="Destino do lixo domiciliar"),
"V0219"=list(type="factor", start=66, width=1,description="Forma de iluminação do domicílio"),
"V0220"=list(type="factor", start=67, width=1,description="Tem telefone móvel celular"),
"V2020"=list(type="factor", start=68, width=1,description="Tem telefone fixo convencional"),
"V0221"=list(type="factor", start=69, width=1,description="Tem fogão de duas ou mais bocas"),
"V0222"=list(type="factor", start=70, width=1,description="Tem fogão de uma boca"),
"V0223"=list(type="factor", start=71, width=1,description="Tipo de combustível utilizado no fogão"),
"V0224"=list(type="factor", start=72, width=1,description="Tem filtro d’água"),
"V0225"=list(type="factor", start=73, width=1,description="Tem rádio"),
"V0226"=list(type="factor", start=74, width=1,description="Tem televisão em cores"),
"V0227"=list(type="factor", start=75, width=1,description="Tem televisão em preto e branco"),
"V2027"=list(type="factor", start=76, width=1,description="Tem aparelho de DVD"),
"V0228"=list(type="factor", start=77, width=1,description="Tem geladeira"),
"V0229"=list(type="factor", start=78, width=1,description="Tem freezer"),
"V0230"=list(type="factor", start=79, width=1,description="Tem máquina de lavar roupa"),
"V0231"=list(type="factor", start=80, width=1,description="Tem microcomputador"),
"V0232"=list(type="factor", start=81, width=1,description="Microcomputador é utilizado para acessar a Internet"),
"V2032"=list(type="factor", start=82, width=1,description="Tem carro ou motocicleta de uso pessoal"),
"V4105"=list(type="factor", start=83, width=1,description="Código de situação censitária"),
"V4107"=list(type="factor", start=84, width=1,description="Código de área censitária"),
"V4600"=list(type="factor", start=85, width=2,description="Dia de referência"),
"V4601"=list(type="factor", start=87, width=2,description="Mês de referência"),
"V4602"=list(type="factor", start=89, width=4,description="Estrato"),
"V4604"=list(type="integer", start=93, width=2,description="Número de municípios selecionados no estrato"),
"V4605"=list(type="numeric", start=95, width=12,description="Probabilidade do município"),
"V4606"=list(type="integer", start=107, width=3,description="Número de setores selecionados no município"),
"V4607"=list(type="numeric", start=110, width=12,description="Probabilidade do setor"),
"V4608"=list(type="factor", start=122, width=6,description="Intervalo de seleção do domicílio"),
"V4609"=list(type="numeric", start=128, width=9,description="Projeção de população"),
"V4610"=list(type="numeric", start=137, width=3,description="Inverso da fração"),
"V4611"=list(type="numeric", start=140, width=5,description="Peso do domicílio"),
"V4614"=list(type="numeric", start=145, width=12,description="Rendimento mensal domiciliar para todas as unidades domiciliares (exclusive o rendimento das pessoas cuja condição na unidade domiciliar era pensionista, empregado doméstico ou parente do empregado doméstico e das pessoas de menos de 10 anos de idade)"),
"UPA"=list(type="numeric", start=157, width=4,description="Delimitação do município"),
"V4617"=list(type="factor", start=161, width=7,description="STRAT - Identificação de estrato de município auto-representativo e não auto-representativo"),
"V4618"=list(type="factor", start=168, width=7,description="PSU - Unidade primária de amostragem"),
"V4620"=list(type="integer", start=175, width=2,description="Número de componentes do domícilio (exclusive as pessoas cuja condição na unidade domiciliar era pensionista, empregado doméstico ou parente do empregado doméstico)"),
"V4621"=list(type="numeric", start=177, width=12,description="Rendimento mensal domiciliar per capita"),
"V4622"=list(type="factor", start=189, width=2,description="Faixa do rendimento mensal domiciliar per capita"),
"V4624"=list(type="factor", start=191, width=1,description="Forma de abastecimento de água"),
"V9992"=list(type="character", start=192, width=8,description="Data de geração do arquivo de microdados")
)
#Endereço dos microdados
pnadFile <- file.path("C:/PNAD2011/Dados", "DOM2011.DAT")

#Especficação para a importação
sourceData <- RxTextData(pnadFile, colInfo=colList)
outputData <- RxXdfData("DOM2011.xdf")
rxImport(sourceData, outputData, overwrite = TRUE)
Note que cada variável possui um tipo definido (integer, numeric ou factor), sua posição incial e final no arquivo texto DOM2011.DAT além de um descritor. Essas informações estão disponíveis no dicionário da PNAD 2011. Outro comentário importante, é o fato da linha número 70 especificar o local em que os dados da PNAD foram extraídos. Uma vez importados esses dados, podemos realizar algumas operações de interesse, por exemplo:
#Confere o formato das variáveis
rxGetInfoXdf("DOM2011.xdf", getVarInfo=TRUE)

#Remove as observações com peso missing ou negativo
rxDataStep(inData = "DOM2011.xdf", outFile = "DOM2011.xdf",rowSelection=(V4611>0),overwrite=TRUE)

#Estatísticas por UF usando somente a média
mun<-rxSummary(~V2016:UF, data = "DOM2011.xdf",
 fweights = "V4611",summaryStats ="Mean")
print(mun)
Nos próximos posts sobre o Revolution R utilizaremos a base de dados criada e citada aqui, qual seja:DOM2011.xdf.

segunda-feira, 15 de junho de 2015

Cálculo do tamanho de amostras: médias.






Frequentemente, desejamos controlar o erro relativo associado a um determinado tamanho amostral, de fato, quando deseja-se trabalhar com médias populacionais ao invés de simplesmente a proporção populacional, podemos considerar o erro relativo como argumento para se determinar o tamanho amostral que o satisfaz (veja por exemplo: Cochran (1977); Schaeffer, et. al (2011), Särndal, et.al (2003)).

Quando desejamos trabalhar com variáveis qualitativas como Gênero, Raça, Satisfação calculamos o tamanho amostral por meio de proporções, já quando a variável de interesse é contínua, como Renda, Altura, Peso, podemos usar a fórmula para o tamanho amostral por meio da média (Cochran (1977) página 77.).

Para a determinação do tamanho amostral mínimo no caso da média, precisamos:

  1. Tamanho da população: é o tamanho da população alvo. Representado usualmente por $N$
  2. Erro relativo: é o erro relativo percentual admitido para o estimador da média. Por exemplo, um erro permissível de 10% indica que a média amostral (obtida após a coleta dos dados) pode diferir no máximo em 10% da média populacional (para mais ou para menos). O erro relativo é dado por $r=0,1$ nesse exemplo.
  3. Confiabilidade: como amostragem é um processo probabilístico, existe uma probabilidade desse erro relativo (ou seja o erro máximo aceitável) não ser satisfeito. Definimos como nível de confiança (confiabilidade) a probabilidade do erro máximo relativo ser satisfeito. Usualmente, trabalha-se com probabilidades como 90%, 95%, 99% ou ainda 99.9% dependendo do tipo de estudo. O nível de confiança (representado por $1-\alpha$ onde $\alpha$ é o nível crítico.) varia entre 0 e 1 (varia entre 0% a 100%)
  4. Coeficiente de variação:o último ingrediente para o cálculo do tamanho amostral necessário é o valor do coeficiente de variação representado por $\hat{CV}=\frac{\hat{S}}{\overline{Y}}$, para a variável de interesse, isso é, é a razão entre o desvio-padrão e a média populacional para a variável de interesse. Frequentemente essa informação não está disponível previamente, nesse caso há no entanto algumas sugestões:
    • Utilize $\hat{CV}$ como sendo aproximadamente a razão entre desvio-padrão e média de uma distribuição uniforme com $\hat{S}=\sqrt{\frac{(Max-Min)^{2}}{12}}$,e $\overline{Y}=\frac{Min+Max}{2}$ onde $Max$ e $Min$ são os valores máximos e mínimos conhecidos para a variável contínua de interesse. Nesse caso, o "pior dos casos" é construído e o tamanho amostral máximo é obtido.
    • Encontre o valor de $\hat{CV}$ utilizando outro estudo. Por exemplo, procure alguma estatística ou algum artigo que indique "mais ou menos" qual deve ser o valor do coeficiente de variação para a variável de interesse. Suponha que o objetivo seja calcular a média do IRA (Índice de Rendimento Acadêmico) dos alunos na UnB, podemos usar algum texto ou artigo que estime o coeficiente de variação do IRA para os alunos de outra instituição de pesquisa e então essa estimativa é utilizada como proxy.
    • Faça uma amostra piloto utilizando um tamanho amostral arbitrário. Com base nessa amostra piloto calcule o valor do desvio-padrão e média para o IRA dos alunos amostrados e então estime o tamanho amostral necessário. Por exemplo, suponha que o objetivo seja, novamente, calcular a média do IRA dos alunos na UnB. Podemos amostrar dez alunos e calcular o valor do coeficiente de variação do IRA desses alunos. Então utilizamos esse valor para estimar o tamanho da amostra. Suponha que encontremos o valor de 234 para o tamanho amostral. Isso significa que temos que amostrar 234 alunos... Mas, como dez alunos já foram amostrados, o tamanho da amostra deverá ser 224=234-10.

Com base nessa breve explicação, considere o seguinte exemplo: o total de alunos matriculados na UnB em 2011 foi de 30757 alunos. Suponha que desejamos fazer uma pesquisa sobre o IRA médio dos alunos da UnB, nesse caso temos que definir algumas informações.

Podemos definir o erro relativo como 0.05, ou seja, admite-se que a média do IRA dos alunos da UnB pode variar em até 5 pontos percentuais para mais ou para menos do IRA médio populacional (isto é, o IRA médio calculado com base no censo de alunos), o nível de confiança mais utilizado é de 95% isso significa que se o processo amostral for repetido muitas vezes espera-se que a margem de erro ±5% seja satisfeita em 95% das vezes. Por fim, como não conhecemos a priori nenhuma informação sobre o desvio-padrão do IRA dos alunos, podemos usar o fato do IRA variar somente de 0 a 5, logo, pela abordagem da distribuição uniforme temos fazer $\hat{S}=\sqrt{\frac{(5-0)^{2}}{12}}=1.443376$ e $\overline{Y}=\frac{5+0}{2}=2.5$ como abordagem conservadora. Fornecendo assim, $CV=\frac{1.443376}{2.5}\approx 0.5773503$

Abaixo segue o programa para o cálculo do tamanho amostral:


Tamanho da população:


Erro relativo:


Coeficiente de variação:


Confiabilidade:


Tamanho da amostra:




Note que ao executar o programa o tamanho da amostra estimado foi de 505. Esse valor poderia ser reduzido se fizéssemos uma amostra piloto ou se tivéssemos uma estimativa menos rígida para o coeficiente de variação.

sexta-feira, 15 de maio de 2015

Revolution R - Parte 2.


Dando continuidade ao nosso estudo sobre o Revolution R, vimos no post Revolution R - Parte 1, como trabalhar inicialmente com a IDE fornecida pelo grupo Revolution Analytics.

Nesse post veremos como manipular grandes bases de dados usando para isso o Revolution R. Para isso, considere o arquivo Pobreza.csv. É claro que essa base não é uma "base grande", mas vamos fazer os exercícios com ela, e, sem perda de generalidade, podemos posteriormente utilizar algum conjunto de dados com maior massa de informações.

O primeiro comando para manipulação de banco de dados no formato *.XDF será o comando rxDataStep. O primeiro passo é ler a base de dados Pobreza.csv:


#Importa os dados no formato XDF
pobreza.xdf <- rxImport(inData = "Pobreza.csv",outFile = "pobreza.xdf", overwrite=TRUE)
Esse comando importa o arquivo Pobreza.csv e cria outro denominado pobreza.xdf o qual estará no formato nativo do Revolution R, em seguida, podemos fazer consultas simples nessa base de dados "grande". Para utilizar o comando rxDataStep, podemos escrevê-lo diretamente ou usando o Snippet:
Aqui a base Pobres.xdf será composta pelos municípios com Índice de Gini (Var09) maior ou igual a 0.40 e Incidência de Pobreza (Var03) superior a 50. Já a base Ricos.xdf será composta pelos municípios com Índice de Gini (Var09) menores do que 0.40 e Incidência de Pobreza (Var03) inferior ou igual a 50. O código é dado por:
#Cria as bases Pobres.xdf e Ricos.xdf
Pobres.xdf<- rxDataStep(inData = "pobreza.xdf", 
 outFile = "Pobres.xdf", rowSelection = Var09>=0.4 & Var03>50)

#Mostra as primeiras linhas
head(Pobres.xdf)

Ricos.xdf<- rxDataStep(inData = "pobreza.xdf", 
 outFile = "Ricos.xdf", rowSelection = Var09<0.4 & Var03<=50)

#Mostra as primeiras linhas
head(Ricos.xdf)
Para facilitar a alocação de memória podemos manter somente algumas variáveis de interesse, por exemplo, suponha que gostaríamos de ter na base somente as variáveis Var01, Var03 e Var09. Para isso, basta fazer:
#Mantêm somente as variáveis de interesse
rxDataStep(inData = "pobreza.xdf", 
 outFile = "Pobres.xdf", rowSelection = Var09>=0.4 & Var03>50,
 overwrite=TRUE,varsToKeep=c("Var01","Var03","Var09"))


#Retira as variáveis que não são necessárias
rxDataStep(inData = "pobreza.xdf", 
 outFile = "Ricos.xdf", rowSelection = Var09<0.4 & Var03<=50,
 overwrite=TRUE,varsToDrop=c("Var02","Var04","Var05",
 "Var06","Var07","Var08","Var10","Var11"))
O comando rxDataStep também pode ser utilizado para se criar novas variáveis na base de dados, por exemplo, suponha que desejamos criar uma variável Tipo que assume valor igual a 0, quando a base é Pobres.xdf e 1 quando a base é Ricos.xdf:
#Cria a variavel Tipo para a base Pobres
rxDataStep(inData = "Pobres.xdf",outFile = "Pobres.xdf",
 transforms = list(Tipo = Var01 < 0),overwrite=TRUE)

#Mostra as primeiras observações
head(Pobres.xdf)

#Cria a variavel Tipo para a base Ricos
rxDataStep(inData = "Ricos.xdf",outFile = "Ricos.xdf",
 transforms = list(Tipo = Var01 > 0),overwrite=TRUE)

#Mostra as primeiras observações
head(Ricos.xdf)
Seguiremos com mais dicas sobre o Revolution R em outros posts.

quarta-feira, 15 de abril de 2015

Text Mining com kernels string no R - Parte 3.


Nesse post daremos continuidade ao que já vimos em posts anteriores como Text Mining com kernels string no R - Parte 1 e Text Mining com kernels string no R - Parte 2. Nessa etapa, vamos tentar reconhecer o padrão das informações prestadas pela Petrobras de 2006 até o primeiro trimestre 2014, para isso, considere os seguintes pacotes:

#Habilita as bibliotecas necessárias
library(tm)
library(wordcloud)
library(Rstem)

Vamos trabalhar com os textos informados pela Petrobras, os quais podem ser obtidos nesse link. Salve esses arquivos em alguma pasta do seu computador, como por exemplo, C:\Text Mining. Em seguida, precisamos transformar os arquivos em PDF para o formato TXT, o qual é possível de ser lido pelo R:

#Seleciona o arquivo de texto a ser minerado
dest<- "C:\\Text Mining\\2006 - 1T.pdf"

#Executa o programa que converte pdf em txt
exe<-"C:\\Program Files\\xpdfbin-win-3.03\\bin64\\pdftotext.exe"
system(paste("\"", exe, "\" \"", dest,"\"", sep= ""), wait= F)

#Cria o arquivo no formato txt
filetxt<- sub(".pdf", ".txt", dest)

#Lê os dados do arquivo
txt<- readLines(filetxt, warn=FALSE)

#Deixa tudo minúsculo
txt<- tolower(txt)

#Remove algumas expressões regulares
txt<-removeWords(txt, c("\\f", stopwords("portuguese")))

#Cria os objetos no formato corpus
corpus<-Corpus(VectorSource(txt))

#Remove a pontuação
corpus<-tm_map(corpus,removePunctuation)

#Remove os números
corpus<-tm_map(corpus,removeNumbers)
Note que o código acima só é possível de ser executado se você tiver instalado o arquivo que converte PDF para TXT, qual seja, Xpdf. Na linha número 5 deve estar o endereço do arquivo executável pdftotext.exe o qual é obtido ao se instalar o programa Xpdf. Cada texto em PDF no arquivo é transformado em TXT e pode ser lido. Para cada texto, podemos associar um retorno obtido no trimestre após a divulgação da informação, por exemplo, para o primeiro trimestre de 2006 temos:
#Coloca os textos na forma de Matriz
tdm<- TermDocumentMatrix(corpus)
m<- as.matrix(tdm)

#Cria um data.frame com a frequência das palavras
d<- data.frame(freq= sort(rowSums(m), decreasing =TRUE))

#Coloca as palavras nas colunas
d$word <- row.names(d)
d$stem<- wordStem(d$word, language="portuguese")

#Remove as palavras muito grandes
d<- d[nchar(row.names(d))<20,]

#Agrega as frequências das palavras
agg_freq<- aggregate(freq ~stem, data = d, sum)
agg_word<- aggregate(word ~stem, data =d, function(x) x[1])

#Une os dados de frequência e palavras
d<- cbind(freq = agg_freq[,2], agg_word)

#Ordena as palavras pela frequência
d<- d[order(d$freq,decreasing =TRUE),]

#1 Trimestre 2006
d$Trimestre<-"1t2006"
d$Retorno<-0.03789927 #Retorno do 2t2006
Uma vez criado um DataFrame para cada texto com o seu respectivo retorno trimestral, podemos usar as frequências das palavras como uma covariável na previsão do retorno trimestral. Essa proposta não é nova e tem sido usada em finanças na busca por padrões. Um bom texto para se começar a estudar esse assunto é o artigo Text mining for market prediction: A systematic review.

segunda-feira, 16 de março de 2015

Text Mining com kernels string no R - Parte 2.


Inicialmente vimos no post Text Mining com kernels string no R - Parte 1 como gerar Wordclouds as quais são úteis na análise descritiva de informação textual, entretanto, podemos fazer algo mais complexo do que simplesmente uma análise descritiva. Podemos reconhecer padrões associados aos textos.

Considere uma lista de 40 textos obtidos da reuters e armazenados no pacote kernlab:

#Habilita o pacote kernlab
library(kernalb)

#Base de notícias da reuters
data(reuters)

#Armazena as palavras chave
y <- rlabels

#Armazena os textos
x <- reuters
Em seguida podemos fazer um Kernel Principal Component Analysis para estudar o padrão das palavras chave:
#Cria o Kernel String
sk <- stringdot(type="spectrum", length=2, normalized=TRUE)

#Faz o Kernel Principal Component Analysis com base nos textos
kpc <- kpca(x,kernel=sk,scale=c())

#Plota o Kernel Principal Component Analysis
plot(rotated(kpc),col=ifelse(y==levels(y)[1],1,2))
Obtendo assim:
O Kernel Principal Component Analysis é uma boa ferramenta para estudar os "tipos" de textos existentes, por exemplo, podemos encontrar textos mais "otimistas", "pessimistas" e atribuir a esses textos um valor quantitativo que pode ser utilizado nos campos de marketing e finanças.

segunda-feira, 16 de fevereiro de 2015

Utilizando RStudio na nuvem com Amazon AMI.


Frequentemente, temos códigos que rodam por horas ou até dias... Impossibilitando o computador para a realização de outras tarefas. Uma solução interessante é "alugar" máquinas que façam essas análises para você.

Hoje está cada vez mais fácil trabalhar com computação em nuvem. Vou mostrar nesse post como utilizar o serviço da Amazon para computação em nuvem com RStudio.

1) O primeiro passo é criar uma conta AWS (Amazon Web Services).

2) Em seguida, vá ao site RStudio Server Amazon Machine Image (AMI). Você pode escolher alguma das regiões disponíveis, para as quais o RStudio já está instalado, como por exemplo, São Paulo:


Nessa etapa você já deve ter sua conta AWS (Amazon Web Services), após preencher com o nome de usário e senha a seguinte tela surge:


Note que há uma lista de possíveis máquinas e configurações disponíveis. Cada máquina possui um valor diferente para o aluguel do serviço. No caso de programação em paralelo deve-se atentar ao número de núcleos (cores) disponíveis e no caso de grandes bases de dados, também na memória RAM disponível.

3) Após escolher a máquina clique no botão Next: Configure Instance Details. Pode aceitar as configurações default e seguir para o próximo passo clicando em Next: Add Storage.

4) Também, nessa etapa (caso não necessite de nenhuma das opções apresentadas) pode seguir para o passo 5) clicando em Next: Tag Instance.

5) Essa é mais uma etapa optativa, não é necessário preencher os campos Key e Value.

6) Essa é mais uma etapa optativa, não é necessário preencher os campos Key e Value. Clique em Next: Configure Security Group.

7) Nessa etapa você deve escolher Type igual a HTTP e nos campos Security group name e Description pode preencher com o nome RStudio:


8) A última etapa para a criação do RStudio Server é clicar no botão Review and Launch. A Amazon irá apresentar as configurações escolhidas e então basta clicar em Launch para criar a instância. Obs: Caso surja uma janela com os dizeres Select an existing key pair or create a new key pair escolha a opção Proceed without a key pair.

Abrindo o RStudio Server na nuvem.

Se as etapas anteriores foram executadas corretamente uma tela como a apresentada abaixo deve surgir:


1)Para abrir o RStudio Server, clique no nome da instância destacada na imagem anterior no retângulo vermelho.

2)Na tela que surgir, procure pelo endereço explicitado em Public DNS. O endereço estará no canto inferior direito da tela. Ao copiar e colar nesse endereço no navegador a seguinte tela aparecerá:


3) Escolha como Username: rstudio e Password: rstudio e então clique em Sign In. Após essas etapas surge:


4) Você deve alterar a senha, para isso, na linha 18, troque o nome "mypassword" pelo nova senha desejada e então execute todo o programa. Pronto! A senha foi atualizada. A partir de agora todo o login deverá ser feito com a nova senha escolhida.

Importando e exportando objetos no RStudio Server.

Caso você tenha alguma base de dados que deseja trabalhar, basta fazer o Upload para o servidor. Basta clicar na opção Upload da aba Files:


Como exercício, fiz o Upload do arquivo german.csv (a descrição da base também está disponível):


Podemos criar alguns novos objetos:

#Importa a base de dados
german.df<-read.csv("german.csv")

#Gera uma variável N(0,1)
x<-rnorm(100,0,1)
Da mesma maneira que podemos fazer o Upload de dados podemos também fazer Download. Por exemplo, se fizermos um histograma:
#Cria o histograma
hist(x)
Podemos salvá-lo em nossa máquina, basta clicar na opção Export:
Em seguida, salve o gráfico na nuvem:
Por fim, escolhemos os objetos que desejamos exportar e clicamos na opção Export:
. Escolha a opção Download:

Interface entre o RStudio Server da Amazon e Dropbox.

Por facilidade, podemos fazer o Download e o Upload de objetos entre o RStudio Server da Amazon e Dropbox. Para isso, além de uma conta no Dropbox, instalar o pacote rDrop. Os passos são: 1) Instale o pacote devtools no RStudio Server da Amazon. 2) Instale o pacote rDrop:
#Habilita o pacote devtools
library(devtools)

#Instala o pacote rDrop
install_github("duncantl/ROAuth")
install_github("karthik/rDrop")
3) Autorize o uso do RStudio Server da Amazon na sua conta do Dropbox. Para isso você precisa criar um App Dropbox developer site. Clique em Create App:
4) Escolha a opção Dropbox Api app:
5) Uma vez que o App foi criado, certifique-se de guardar sua App key e App secret em algum lugar seguro.
#Habilita o pacote rDrop
library(rDrop)

#Insere a App key necessária para a utilização do dropbox
dropbox_credentials <- dropbox_auth("App key", "App secret")
6) Podemos importar dados do Dropbox e exportar dados, basta usar os comandos:
#Salva arquivos no dropbox
?dropbox_save

#Importa arquivos do dropbox
?dropbox_get
Por fim, por se tratar de um serviço PAGO após terminar o uso do RStudio Server Amazon você precisa encerrar a instância, ou então vai continuar pagando pelo serviço. Para isso, basta ir no site AWS onde o servidor foi configurado e na opção Instances parar a instância:
Também é possível encerrar a instância por meio do código:
#Enccerra a instância do RStudio Server na Amazon:
system("sudo shutdown -h now", wait = FALSE)

quinta-feira, 15 de janeiro de 2015

Processamento e paralelo no R.


A computação paralela é uma forma de computação na qual muitos cálculos que são realizados em simultaneamente, aumentando assim a velocidade de execução de determinados códigos. Esse tipo de processamento opera no princípio de que grandes problemas muitas vezes podem ser divididos em partes menores, que são então resolvidos simultaneamente ("em paralelo"). O possível aumento de velocidade máxima de um único programa, como resultado de paralelização é conhecida como Lei de Amdahl.

No R podemos trabalhar com a computação em paralelo por meio de dois pacotes, a saber: foreach e doParallel.

Inicialmente é necessário instalar e carregar os pacotes de interesse:

#Carrega os pacotes necessários para realizar o paralelismo
library(foreach)
library(doParallel)

Cada núcleo existente na sua máquina (ou uma parte deles) pode ser utilizado para se dividir as tarefas e, consequentemente, os cálculos, para isso, precisamos saber quantos núcleos temos disponíveis:

#Checa quantos núcleos existem
ncl<-detectCores()
ncl

#Registra os clusters a serem utilizados
cl <- makeCluster(ncl)
registerDoParallel(cl)
Note que podemos registrar menos clusters do que temos disponíveis, caso você não deseje que sua máquina fique "travada" realizando somente as operações de cálculo exigidas. Para testar a velocidade que o processamento em paralelo trás para o nosso código vamos inicialmente gerar uma base de dados:
#Gera o número de observações
n<-1000

#Variáveis geradas
x<-rnorm(n,0,1)
y<-rnorm(n,1+2*x,2)

#Dataframe com as variáveis geradas
dados<-data.frame(x,y)
Em seguida, vamos realizar o Bootstrap usando a função for usual do R para computar o Erro-Padrão dos parâmetros do modelo de regressão na forma: $y=\alpha+\beta x +\epsilon$ O código é dado por:
#Inicia a contagem do tempo
ptm <- proc.time()

#Cria o vetor para armazenar o parãmetro beta em cada iteração Bootstrap
beta<-rep(0,5000)

#Faz o Bootstrap usando a função for
for(i in 1:5000)
{
  #Gera a amostra Bootstrap
  bdados<-dados[sample(nrow(dados),nrow(dados),replace=T),]
  beta[i]<-unname(lm(y~x,bdados)$coef[2])
}
mean(beta)
sd(beta)

#Para de contar o tempo
proc.time() - ptm
Os resultados para a minha máquina com 4 núcleos e utilizando for foi:
Agora podemos usar a função foreach sem paralelismo para comparar também:
#Inicia a contagem do tempo
ptm <- proc.time()

#Cria o vetor para armazenar o parãmetro beta em cada iteração Bootstrap
beta<-rep(0,5000)

#Faz o Bootstrap usando a função foreach
boot_b <- foreach(i=1:5000, .combine=c) %do% {
  #Gera a amostra Bootstrap
  bdados<-dados[sample(nrow(dados),nrow(dados),replace=T),]
  beta[i]<-unname(lm(y~x,bdados)$coef[2])
}
mean(beta)
sd(beta)

#Para de contar o tempo
proc.time() - ptm
Os resultados foram mais demorados do que utilizando somente o for:
Finalmente, podemos considerar realizar essas tarefas em paralelo, distribuindo os cálculos entre os núcleos registrados:
#Inicia a contagem do tempo
ptm <- proc.time()

#Cria o vetor para armazenar o parãmetro beta em cada iteração Bootstrap
beta<-rep(0,5000)

#Faz o Bootstrap usando a função foreach
boot_b <- foreach(i=1:5000, .combine=c) %dopar% {
  #Gera a amostra Bootstrap
  bdados<-dados[sample(nrow(dados),nrow(dados),replace=T),]
  beta[i]<-unname(lm(y~x,bdados)$coef[2])
}
mean(boot_b)
sd(boot_b)

#Para de contar o tempo
proc.time() - ptm

#Stop clusters
stopCluster(cl)

Usando o paralelismo, o Bootstrap foi muito mais rápido:
Nota-se que o processamento em paralelo realmente é vantajoso, mas outras funções como a família apply, sapply, lapply, mapply, etc. também são muito boas quando deseja-se que o código rode o mais rapidamente possível.