All Posts By

Diogo Wolff

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.


							
						

Relatório #21 - Expectativa

By | Indicadores

Para o Relatório AM dessa semana, decidimos comentar um pouco sobre a trajetória das expectativas do mercado para as principais variáveis macroeconômicas. Para iniciar, deixamos abaixo as expectativas iniciais para o IPCA:

A história é aquela já contada muitas vezes nos últimos meses: choques que já eram vistos no final do ano passado custam a se desfazer, como a escassez de componentes eletrônicos para automóveis e os preços de commodities, e a perspectiva é de que a política monetária tem se mantido um passo atrás da trajetória de inflação, aumentando a taxa de juros porém sem o fazer a um nível satisfatoriamente impactante para controlar os preços. Essa aparente desconfiança sobre o fim de tais choques - e sobre o comprometimento do BC - já contamina as expectativas do ano que vem, se elevando continuamente nas últimas semanas.

Simultaneamente, podemos ver abaixo que o juro real, que foi mantido em níveis baixos no pico da pandemia para tentar aquecer a economia, está agora se aproximando cada vez mais do juro neutro, apontando para uma perspectiva mais neutra do BC e em favor da estabilização dos preços.

O resultado dessa aproximação, em conjunto com uma série de fatores de instabilidade no país, é uma perspectiva negativa para o ano que vem. Enquanto o PIB para esse ano possui expectativas favoráveis - muito por conta de uma comparação com um 2020 péssimo -, 2022 já aponta para apenas cerca de 1.5% de crescimento.

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.

Relatório #20 - PMS e IBC-Br

By | Indicadores

No final dessa semana, são anunciados os valores da PMS e do IBC-Br referentes ao mês de agosto, completando as principais pesquisas sobre o nível de atividade da economia. A PMS é particularmente importante pois o setor de serviços teve alta variância durante a pandemia, e há grande expectativa sobre a manutenção de sua recuperação. Como podermos ver abaixo, após a queda causada pelas restrições do começo desse ano, houve uma trajetória de crescimento, porém ela já está diminuindo, indicando que estamos convergindo para um nível estável de comércio.

Podemos ver que, em comparação com as outras pesquisas, a PMS apresenta uma variação anual menor, indicando que a recuperação está atrasada para o setor de serviços. É importante notar, porém, que esse setor ficou abaixo dos valores anteriores por um período prolongado, o que explica seu vale mais profundo no começo de 2021.

Mudando o foco agora para o IBC-Br, que sai nessa sexta-feira, podemos ver no gráfico abaixo que a oferta da economia como um agregado já parece estar recuperada. Observando seu valor em nível, tanto os valores em nível como dessazonalizados aparentam estar levemente abaixo da tendência pré-pandemia, com a variação anual indicando que os impactos do lockdown de 2020 já foram superados pela recuperação que ocorreu na segunda metade do ano.

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.

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":""}