Category

Hackeando o R

Hackeando o R: Modelos lineares mistos

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar como realizar a estimação de um modelo linear misto. O método deriva seu nome devido à possibilidade de um modelo apresentar tanto efeitos fixos como aleatórios, porém é bom fazer um aviso: a nomenclatura utilizada nesse modelo é proveniente da bioestatística, que chama de efeitos fixos aqueles que são invariantes na população, e efeitos aleatórios aqueles que possuem uma distribuição de probabilidade que assume um valor para cada indivíduo.

Para compreender essas definições, vamos utilizar um exemplo do pacote lme4. Abaixo, fazemos uma regressão simples sobre os dados do dataset sleepstudy. Ele é proveniente de um estudo onde indivíduos foram restritos a 3 horas de sono por noite, e a cada dia seu tempo de reação foi observado.

library(lme4)
library(ggplot2)

ggplot(sleepstudy, aes(Days, Reaction, colour = Subject)) +
geom_point() +
geom_smooth(method = "lm", colour = "black", se = FALSE)

Como podemos ver, há indivíduos que se mantém consistentemente abaixo da reta de regressão, pois devem possuir tempo de reação menor de modo generalizado, enquanto há outros que estão sempre acima. Isso reflete o fato de que a amostra não possui observações independentes (pois é feita para cada indivíduo ao longo de 10 dias), e não podemos utilizar o modelo linear clássico. A solução oferecida por um modelo linear misto é gerar um componente aleatório que é específico a cada indivíduo, de modo que podemos ter, para cada um, um intercepto e uma inclinação distintas. Para fazer isso, basta utilizar a função lmer:

m4 <- lmer(Reaction ~ Days + (1|Subject) + (0 + Days|Subject),
data = sleepstudy)

summary(m4, correlation = FALSE)

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
Data: sleepstudy

REML criterion at convergence: 1743.7

Scaled residuals:
Min 1Q Median 3Q Max
-3.9626 -0.4625 0.0204 0.4653 5.1860

Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 627.57 25.051
Subject.1 Days 35.86 5.988
Residual 653.58 25.565
Number of obs: 180, groups: Subject, 18

Fixed effects:
Estimate Std. Error t value
(Intercept) 251.405 6.885 36.513
Days 10.467 1.560 6.712

O modelo acima possui tanto intercepto e inclinação populacionais, como também coeficientes correspondentes para cada um dos indivíduos. No formato acima, supomos que os coeficientes individuais não são correlacionados (note que são escritos separadamente), o que pode ser testado, porém está fora do escopo do post de hoje. A partir desse modelo, podemos visualizar as retas de regressão, e compará-las com a reta que seria estimada ao fazer uma regressão para cada indivíduo:

mixed_mod <- coef(m4)$Subject
mixed_mod$Subject <- row.names(mixed_mod)

ggplot(sleepstudy, aes(Days, Reaction)) +
geom_point() +
theme(legend.position = "none") +
facet_wrap(~ Subject, nrow = 3) +
geom_smooth(method = "lm", colour = "cyan3", se = FALSE,
size = 0.8) +
geom_abline(aes(intercept = `(Intercept)`, slope = Days,
color = "magenta"),
data = mixed_mod, size = 0.8) +
theme_bw()

Como podemos ver, as retas são diferentes (sendo a azul a regressão OLS), o que ocorre pois as regressões do modelo misto são aproximadas em direção à média de todas elas.


							
						

Hackeando o R: Análise de variáveis latentes

By | Hackeando o R

No Hackeando o R de hoje, vamos fazer uma expansão sobre o conteúdo do nosso post de semana passada, sobre análise de fatores. Dentro do contexto de fatores, a classificação com um único fator é equivalente ao clustering, e, considerando que as observações dentro de um mesmo cluster não são correlacionadas, entramos em um modelo chamado de análise de perfil latente. No post de hoje, consideraremos o caso em que as variáveis conhecidas são categóricas, chamado de análise de classes latentes. A ideia é podermos classificar indivíduos em grupos mutuamente exclusivos, a partir de suas respostas a perguntas categóricas. Nesse sentido, podemos ver que esse modelo é útil para a análise de dados de censos, como veremos no exemplo a seguir.

Para apresentarmos o modelo na prática, utilizaremos a base de dados do censo de instituições de saúde mental N-HMSS, disponível aqui. O modelo em si está disponível no pacote poLCA, logo iremos carregá-lo em conjunto com os dados abaixo:

library(poLCA)

seed = 1

data = read.csv('nmhss-puf-2019.csv')

No exemplo, vamos nos limitar aos tratamentos que estão disponíveis em cada instituição, identificados nas colunas 17 a 30. Além disso, precisamos fazer duas mudanças nos dados: números negativos querem dizer NA, e variáveis binárias precisam ser modificadas de 0 e 1 para 1 e 2 - uma limitação do pacote poLCA.

tratamento = data[, names(data)[17:30]]
tratamento[tratamento < 0] <- NA

tratamento <- tratamento + 1

summary(tratamento)

 

Para começar a análise, vamos procurar classes se baseando na oferta de 5 tratamentos: psicoterapia individual (TREATPSYCHOTHRPY), terapia familiar/de casal (TREATFAMTHRPY), terapia em grupo (TREATGRPTHRPY), terapia cognitivo-comportamental (TREATCOGTHRPY) e medicação psicotrópica (TREATPSYCHOMED).

A função poLCA requer 3 inputs: uma fórmula, um dataframe, e o número de classes a procurar. Abaixo, rodamos ela com 2 classes. O resultado é reportado automaticamente. Como podemos ver, para a primeira variável, temos que a probabilidade da segunda classe ser 2 é de 99.27%, ou seja, uma alta probabilidade de que há o tratamento, enquanto que a primeira classe apresenta apenas 33.72% de probabilidade de fornecer ele. Para todas as variáveis temos uma probabilidade maior para a classe 1 do que para a classe 2, indicando que a segunda classe abrange clínicas com maior variedade de tratamentos. Ademais, temos que é estimado que 88.33% das clínicas pertencem à segunda classe.

m <- poLCA(cbind(TREATPSYCHOTHRPY, TREATFAMTHRPY,
TREATGRPTHRPY, TREATCOGTHRPY,
TREATPSYCHOMED) ~ 1,
data = tratamento, nclass = 2)

Podemos visualizar o resultado de modo mais intuitivo com a função plot():


plot(m)

Para finalizar, fazemos algumas adições ao modelo para mostrar outros pontos importantes. Abaixo, trabalhamos com 3 classes, e, para garantir a convergência da estimação (algo que nem sempre ocorre com a função), aumentamos o número de iterações de 1000 para 2500, e rodamos ele 5 vezes. Ademais, mudamos a fórmula para incluir a variável assistência ao pagamento (PAYASST), de modo a identificar como as classes que encontramos se relacionam com ela.


tratamento$PAYASST <- data$PAYASST
tratamento$PAYASST[tratamento$PAYASST < 0] <- NA

m <- poLCA(cbind(TREATPSYCHOTHRPY, TREATFAMTHRPY,
TREATGRPTHRPY, TREATCOGTHRPY,
TREATPSYCHOMED) ~ PAYASST,
data = tratamento, nclass = 3,
maxiter = 2500, nrep = 5)

Como vemos no resultado, temos 3 classes com características distintas: a classe 1 (cerca de 12% da população) tem diversos tratamentos, exceto terapia familiar; a 2 e a 3 (79 e 9% da população, respectivamente), incluem todos os tratamentos, exceto pelo fato de que a 2 tem probabilidade muito menor de incluir medicação. Ademais, podemos verificar a relação que cada classe tem com a assistência ao pagamento: em comparação com a classe 1, clínicas da classe 2 tendem a oferecer menos assistência, enquanto clínicas da classe 3 tendem a oferecer mais.

Exportando data frames do R para planilhas

By | Hackeando o R

Imagine-se na seguinte situação: o seu chefe lhe pede uma análise de um conjunto de dados "sujos", que requerem um enorme trabalho de limpeza, além de que se faça cálculos complexos. Você, inteligente e inovador, pretende realizar todo o trabalho com o R, já que a linguagem permite que você faça o trabalho rapidamente. Há um porém, o seu chefe ainda não conhece muito bem a linguagem, e também necessita passar as informações para outros setores da empresa, e por isso, pede que o conjunto de dados seja entregue em uma planilha. Neste texto, iremos mostrar como é possível resolver esse problema, exportando data frames para planilhas utilizando os pacote {XLConnect} e {googlesheets4}.

# Carrega os pacotes necessários
library(googlesheets4)
library(XLConnect)
library(tidyverse)

Primeiro, iremos mostrar como exportar uma data frame para um arquivo .xlsx. O primeiro passo é carregar o arquivo
de formato .xlsx, para que seja feita a conexão com o R. Veja que é possível conectar a uma planilha já existente, utilizando a função loadWorkbook(), passando nome do arquivo em questão no primeiro argumento, bem como criar um arquivo direto do R, passando o argumento create = TRUE. Fazemos isso tudo salvando em um objeto.

Logo em seguida, realizamos a criação de páginas dentro da planilha, no qual iremos colocar nossos dados. Podemos fazer isso com a função createSheet(), em que passamos o objeto que conecta com o arquivo .xlsx, e o nome da página. Após isso, utilizamos a função writeWorksheet() para enfim escrever os dados na planilha. Veja que podemos escolher as coordenadas em que os dados estarão guardados usando os argumentos startRow e startCol, que por padrão começam na linha 1, coluna 1.

Por fim, salvamos o arquivo utilizando a função saveWorkbook().


# Cria um arquivo .xlsx
flowers <- loadWorkbook("iris.xlsx", create = TRUE)

# Cria uma página no arquivo
createSheet(flowers, "iris1")

# Escreve na página criada
writeWorksheet(flowers, iris, sheet = "iris1", startRow = 1, startCol = 1)

# Salva o arquivo
saveWorkbook(flowers)

Outro pacote muito útil para transformar data frames em planilhas é o {googlesheets4], que justamente conecta os dados em uma planilha no Google Planilhas.

Para isso, é necessário que o pacote tenha permissão da sua conta Google para realizar as mudanças. Para isso, utiliza-se a função gs4_auth(), colocando seu e-mail de uso.

Em seguida, cria-se a planilha utilizando a função gs4_create(), passando o nome da planilha como primeiro argumento, e o data frame que podem ser utilizados.

Em caso de alterações, podemos escrever por cima da planilha utilizando a função write_sheet(). O argumento ss é usado para conectar a planilha, podemos usar tanto a url, quanto o id salvo dentro do argumento iris_gg.

## Google Sheets

# Conecta a conta do Google, para garantir a permissão
gs4_auth(email = "luiz@exemplo.com")

# Cria uma planilha dentro do Google Planilhas
iris_gg <- gs4_create("flowers", sheets = iris)

# Escreve por cima do arquivo criado anteriormente
write_sheet(data = filter(iris, Species == "setosa"),
            ss = "1V2722XxqQHbi_9V_bi5LAxwEHrQTXHCW4OnpyjjxiQg",
            sheet = "iris")

Serve como um bom quebra-galho, não?

Hackeando o R: Fazendo análise de fatores no R

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar como fazer a análise de fatores no R. A ideia por trás dessa metodologia é tentar identificar a relações não-observáveis em amostras de variáveis observáveis (como uma correlação). Esse método é comum em psicometria, logo para o exemplo utilizaremos o pacote psych, que possui ferramentas específicas para isso. Também utilizaremos o pacote GPArotation, que é útil para a análise de fatores. Diferentemente do PCA, onde é pressuposto que toda a variância é comum entre as variáveis, aqui supomos que há alguma variância em cada observação que é única a ela, através de fatores que são específicos, ou erros. Com isso, vamos fazer um exemplo usando o dataset attitude, que já vem no R base.

Primeiramente, vamos carregar os pacotes que iremos utilizar:

library(psych)
library(GPArotation)

Para criarmos uma análise simples, podemos iniciar identificando o possível número de fatores que explicam a variação comum das observações, utilizando um gráfico scree. Sua interpretação é parecida com a do PCA: estamos buscando uma dobra brusca no gráfico (como um cotovelo), indicando que a passagem de n fatores para n+1 traz pouca informação nova. Isso pode ser feito com a função fa.parallel, que gera o gráfico e algumas informações adicionais. Como podemos ver abaixo, utilizar 2 fatores parece adequado.

fa.parallel(attitude, fm = "ml", fa = "fa")

Com isso, podemos rodar o modelo (cujos detalhes serão omitidos por hoje) usando a função fa(). Dois argumentos merecem atenção: o rotate = "oblimin" indica que permitimos que os fatores possuam correlação, enquanto que fm = "ml" faz a estimação através de máxima verossimilhança.

attitude_fa <- fa(attitude, nfactors = 2,
rotate = "oblimin", fm = "ml")

attitude_fa
fa.diagram(attitude_fa, simple = FALSE)

O diagrama gerado pela fa.diagram() mostra que os dois fatores gerados possuem correlação com características observáveis distintas. O ML1 é correlacionado com learning, raises e advance, o que pode ser interpretado como as oportunidades de carreira para funcionários da empresa. O ML2 tem correlação com todas as variáveis exceto advance e critical, o que pode ser visto como o bem-estar geral dos funcionários. É importante notar que essa é apenas uma interpretação: o objetivo da análise é justamente identificar fatores que explicam variações comuns, e tentar descrevê-los de uma forma causal.

Como retirar dados de contas públicas municipais via API do SICONFI

By | Hackeando o R

Neste texto, iremos mostrar como podemos retirar os dados dos demonstrativos contábeis de entes federativos do Brasil pelo Sistema de Informações Contábeis e Fiscais do Setor Público Brasileiro (SICONFI) via API no R.

Primeiramente, é necessário estar a par dos parâmetros que devem ser colocados como entradas para obter os dados dos diversos tipos de demonstrativos. O site http://apidatalake.tesouro.gov.br/docs/siconfi/ fornece detalhadamente quais parâmetros devem ser fornecidos para cada tipo de demonstrativos, bem como a url base para realizar a requisição do API.

Aqui iremos trabalhar como exemplo a Declaração de Contas Anuais (DCA) Anexo I-D do município de Varginha - Minas Gerais, no qual nos fornecerá as Despesas Orçamentárias por Natureza.

Para o DCA, há 3 parâmetros que devem ser inseridos: an_exercicio (Ano de exercício do demonstrativo); no_anexo (Qual anexo do relatório deseja obter) e id_ente (O código IBGE do ente em questão). Sendo an_exercicio e id_ente obrigatórios para esse demonstrativo em questão.

É fundamental a utilização dos pacotes a seguir.


library(httr)
library(jsonlite)
library(magrittr)
library(tibble)

Em seguida vamos realizar a chamada da API criando uma URL.

# URL da DCA de Varginha no ano de 2020 
url_dca <- paste("https://apidatalake.tesouro.gov.br/ords/siconfi/tt/dca?", # URL base para a chamada 
"an_exercicio=", 2020, "&", # Insere o parâmetro de Ano do exercício 
"no_anexo", "DCA-Anexo+I-D","&", # Insere o parâmetro do Anexo que se deseja obter 
"id_ente=", "3170701", sep = "") # Insere o parâmetro do Ente de acordo com o código IBGE do mesmo

Após isso, devemos realizar a requisição da API usando a função GET do pacote httr, bem como realizar a extração do conteúdo com as funções content e fromJSON dos pacotes httr e jsonlite, respectivamente.


api_dca <- GET(url_dca) 
 
# A chamada irá nos retornar os dados requisitados. Agora só precisamos extrair o conteúdo que nos interessa

json_dca <- api_dca %>%

content(as = "text", encoding = "UTF-8") %>%

fromJSON(flatten = FALSE)

# E após isso transforma-los em um tibble

dca_tb <- as.tibble(json_dca[["items"]])

Desta forma podemos obter os dados do DCA Anexo I-D do município de Varginha. O método pode ser replicado para outros anexos e demonstrativos, bem como para qualquer outro Ente do Brasil.

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente
{"cart_token":"","hash":"","cart_data":""}