Segunda Prova de Tópicos Especiais em Econometria

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

1 ##Primeira questão

#Microdados: O que é uma amostra complexa, e qual é a implicação de “complexa” para a análise econométrica? Questões relacionadas à análise estatística de dados na utilização de informações coletadas em inquéritos com amostras complexas, uma das preocupações refere-se à análise apropriada dos dados, que não podem ser tratados como se fossem observações independentes e identicamente distribuídas (amostragem aleatória simples). Na análise, os seguintes aspectos devem ser considerados: • probabilidades distintas de seleção das unidades; • estratificação; • efeito de conglomeração das unidades; • calibração para ajuste de não resposta e outros ajustes para totais conhecidos da população. Elementos dos desenhos complexos de amostragem Na amostragem probabilística, cada indivíduo da população tem uma probabilidade conhecida e não nula de ser selecionado. São vários os métodos de amostragem probabilística, sendo os mais utilizados a amostra aleatória simples, a amostragem estratificada, e a amostragem por conglomerados. A amostragem estratificada consiste em dividir a população em subgrupos homogêneos para determinadas características e selecionar uma amostra em cada um deles, separadamente. Estes subgrupos não se interceptam e totalizam a população. Cada uma das subdivisões populacionais é denominada de estrato. Este tipo de amostragem é recomendado quando se deseja obter estimativas com certa precisão para cada uma das subdivisões. Cochran demonstra que o uso de estratos homogêneos pode resultar em um ganho na precisão de estimativas de características da totalidade da população. Intuitivamente, se cada estrato for homogêneo, a estimativa de uma determinada característica pode ser obtida dentro do estrato com um número relativamente pequeno de observações. Estas estimativas são então ponderadas, adequadamente, para obtenção da estimativa populacional. Na amostragem por conglomerados, a unidade amostral é um conjunto de elementos da população. Pelas dificuldades de custo e de obtenção de listas dos elementos da população, a aplicação da amostragem por conglomerados é muito frequente nos inquéritos de base populacional. Geralmente, este tipo de amostragem é repetido em múltiplos estágios, quando uma amostra de conglomerados é selecionada em cada etapa. Os conglomerados são unidades compostas de subunidades, que vão sendo selecionadas em cada estágio até chegar ao objeto de interesse da pesquisa. Por exemplo, no Brasil podem-se selecionar no primeiro estágio os setores censitários, denominados de unidades primárias de seleção. No segundo estágio, dentro dos setores selecionados, os domicílios, e no último o indivíduo a ser pesquisado. Particularmente, um método frequentemente utilizado para a seleção dos conglomerados é a amostragem com probabilidade proporcional ao tamanho (PPT). Neste caso, as unidades primárias são selecionadas com probabilidades proporcionais ao seu tamanho. Este método, além da facilidade de aplicação, tem a vantagem de contribuir para a redução da variância entre as unidades de seleção. Além dos tipos de amostragem e da combinação entre eles, constituindo desenhos complexos de amostragem, outro aspecto merecedor de atenção são as probabilidades desiguais de seleção das unidades em cada um dos estágios. Para compensar as probabilidades desiguais de seleção, são atribuídas ponderações diferenciadas aos elementos da amostra, chamadas de fatores naturais de expansão, correspondentes ao inverso do produto das probabilidades de inclusão nos diversos estágios de seleção. Finalmente, para ajustar os pesos naturais do desenho e/ou corrigir os problemas originados pela ausência ou recusa de resposta é necessário fazer a calibração para totais conhecidos da população. A calibração dos fatores naturais de expansão consiste em estimar novos pesos para cada elemento da amostra, através de ajuste dos pesos naturais do desenho segundo informações de variáveis auxiliares da amostra. Por exemplo, um dos propósitos da calibração é obter distribuições amostrais das variáveis auxiliares semelhantes às distribuições populacionais obtidas pelo censo. O plano de amostragem, processo usado para selecionar a amostra da população, é importante para a estimação dos parâmetros, já que a estimação baseada no plano amostral é fundamentada na ponderação de cada elemento da amostra pelo inverso da sua probabilidade de seleção. Quando os dados são utilizados para fins descritivos, as estimativas pontuais dos parâmetros populacionais (como totais, médias, proporções) são influenciadas apenas pelos pesos amostrais determinados pelo plano de amostragem. Na maioria das vezes, porém, são usados com fins analíticos, requerendo a estimação de medidas de variabilidade (como a variância e o desvio padrão) para o cálculo dos intervalos de confiança e a realização de testes de hipóteses. Nos estudos analíticos, as estimativas pontuais são associadas à precisão, dada pelo tamanho do intervalo de confiança. Quanto maior o intervalo, ou quanto maior o erro padrão, menor a precisão da estimativa. A estimação do erro padrão é, pois, essencial para a inferência estatística. Em planos de amostragem complexa, a estimação da variância é influenciada não somente pelos pesos, mas também pela estratificação e conglomeração, conjuntamente. No que se refere à estratificação, se esta for adequadamente utilizada, a variância da média estimada é, geralmente, inferior à variância se a amostragem fosse aleatória simples, ou seja, são obtidos ganhos de precisão nas estimativas médias. Já o efeito da conglomeração é, em geral, inverso, produzindo variâncias maiores do que as calculadas com observações independentes. Na amostragem por conglomerados, a variância da média é função direta da variância total e do coeficiente de correlação intra-classe, que mede a correlação entre as observações dentro dos conglomerados. Significa dizer que quanto maior a homogeneidade dentro dos conglomerados, maior a variabilidade entre conglomerados, e maior a variância da média. Portanto, a amostragem por conglomerados acarreta, em geral, a perda de precisão das estimativas médias. Para medir o efeito do plano amostral (EPA) sobre a variância das estimativas médias, utiliza-se o chamado “efeito de desenho”, calculado pela razão entre a estimativa da variância determinada pelo plano amostral e a estimativa da variância obtida por uma amostra aleatória simples de mesmo tamanho. Além de ser utilizado para verificar a perda de precisão da estimativa, o efeito de desenho é utilizado também para o planejamento de estudos futuros, no cálculo do tamanho de amostra.

2 ##Segunda questão

#Painel: explique as vantagens e desvantagens de se trabalhar com dados em painel comparando os modelos de efeitos fixos, o de efeitos aleatórios e o de dados empilhados. A principal vantagem da utilização de modelos de dados em painel refere-se ao controle da heterogeneidade individual, ou seja, à possibilidade de se medirem separadamente os efeitos gerados por conta de diferenças existentes entre cada observação em cada cross-section, além de ser possível avaliar a evolução, para um dado indivíduo, das variáveis em estudo ao longo do tempo. Os dados em painel providenciam uma maior quantidade de informação, maior variabilidade dos dados, menor colinearidade entre as variáveis, maior número de graus de liberdade e maior eficiência na estimação. A inclusão da dimensão em cross-section, num estudo temporal, confere uma maior variabilidade aos dados, na medida em que a utilização de dados agregados resulta em séries mais suaves do que as séries individuais que lhes servem de base. Esse aumento na variabilidade dos dados contribui para a redução de uma eventual colinearidade existente entre variáveis. Efeitos Fixos Os modelos de efeitos fixos apresentam a complicação adicional de que os regressores sejam correlacionados com os efeitos do nível do indivíduo e, portanto, uma estimação consistente dos parâmetros do modelo requer uma eliminação ou controle dos efeitos fixos. Assim, um modelo que leva em conta os efeitos específicos do indivíduo para uma variável dependente Yit. A vantagem do modelo de efeitos fixos é que pode ser obtido um estimador consistente do efeito marginal do j-ésimo regressor de E(yit|β0i, xit), dado que xj,it varia no tempo. Efeitos Aleatórios No modelo de efeitos aleatórios, por outro lado, assume-se que β0i é puramente aleatório, ou seja, não é correlacionado com os regressores. A estimação, portanto, é elaborada com um estimador FGLS (feasible generalized least squares). A vantagem do modelo de efeitos aleatórios é que este estima todos os coeficientes, mesmo dos regressores invariantes no tempo, e, portanto, os efeitos marginais. Ademais, E(yit|xit) pode ser estimado. Porém a grande desvantagem é que estes estimadores são inconsistentes se o modelo de efeitos fixos for mais apropriado. Conforme já discutido, a variável dependente e os regressores podem potencialmente variar simultaneamente ao longo do tempo e entre indivíduos. Enquanto a variação, ao longo do tempo ou para um dado indivíduo, é conhecida por within variance, a variação entre indivíduos é chamada de between variance. De acordo com Wooldridge (2010), no modelo de efeitos fixos, o coeficiente de um regressor com baixa variação within será imprecisamente estimado e não será identificado se não houver qualquer within variance. Assim, é de fundamental importância a distinção entre essas variações para a definição do melhor modelo de dados em painel. O estimador de efeitos aleatórios será consistente e completamente eficiente se o modelo de efeitos aleatórios for apropriado, porém será inconsistente se o modelo de efeitos fixos for apropriado, uma vez que a correlação entre xit e β0i resulta numa correlação entre os regressores e o termo do erro na equação. Da mesma forma, ainda segundo Cameron & Trivedi (2009), se não ocorrerem efeitos fixos, então o estimador de efeitos aleatórios será consistente, mas ineficiente e, portanto, uma estimação com erros-padrão robustos clusterizados deveria ser obtida. A expressão da estimativa por mínimos quadrados generalizados factíveis de um coeficiente de regressão do modelo, supondo-se efeitos aleatórios, torna-se igual à da estimativa do mesmo coeficiente em um modelo de efeitos fixos (estimação within) se θ = 1. Empilhado (Pooled) Se não ocorrerem efeitos fixos, mas os erros demonstrarem correlação dentro do painel, então o estimador de efeitos aleatórios será consistente, mas ineficiente e, portanto, uma estimação com erros-padrão robustos clusterizados deverá ser obtida. Dessa forma, para um painel curto, em que T < N, uma estimação com erros-padrão robustos clusterizados pode ser obtida considerando-se a premissa de que os erros são independentes entre indivíduos e que N = ∞, ou seja, que ((epsilon)it, (epsilon)js) = 0 para i ≠ j, que E(epsilon)it, (epsilon)js não seja restrita e que (epsilon)it seja heterocedástico. De acordo com Cameron e Trivedi (2009), o passo inicial para a aplicação de um modelo com dados em painel é a aplicação de um modelo POLS (pooled ordinary least squares), que assume que os regressores sejam exógenos e que o termo de erro seja (Mu)it, em vez da decomposição (alfa)i + (epsilon)i. Os parâmetros desse modelo são estimados por meio de OLS, mas a inferência requer que haja controle da correlação within do erro (Mu)it para um dado indivíduo, a ser elaborado utilizando-se erros-padrão robustos com agrupamento no nível do indivíduo. Para dados em painel longo, ou seja, com muitos períodos para um número relativamente menor de indivíduos, os efeitos individuais β0i podem ser incorporados em xit como variáveis dummy para cada período, de modo a haver muitos efeitos de tempo yt (efeitos mensais, trimestrais ou anuais, por exemplo). Um modelo pooled, para T > N, em que os regressores xit contemplam o intercepto, o efeito temporal e, possivelmente, um vetor de variáveis de indivíduo, como T é maior do que N passa a ser necessária a especificação de um modelo que considere a existência de correlação serial do erro (Beck e Katz, 1995). Dessa forma, para dados em painel longo, os modelos pooled com métodos de estimação OLS (POLS) e FGLS passam a ser mais adequados, já que permitem a utilização de um modelo AR(1) para (Mu)it ao longo do tempo, em que (Mu)it seja heterocedástico (Hoechle, 2007).

3 ##Terceira Questão: Microdados:

Para esta questão, baixe os dados do ano de 2019 (dados anuais da primeira visita) da PNADC, microdados consolidados para Brasil e variáveis selecionadas. Pede-se: baixar os dados; organizar e estabelecer o plano amostral complexo; investigar os determinantes do rendimento mensal efetivo de todos os trabalhos. Avaliar o ajuste em termos de significância, impactos das variáveis preditoras, presença de heterocedasticidade e correlação serial. Vocês utilizarão os dados anuais acumulados da 1ª visita. Para baixar os dados anuais, é preciso colocar o número da entrevista (interview), entre 1 e 5 visitas, conforme expresso em https://www.ibge.gov.br/estatisticas/sociais/habitacao/17270-pnad-continua.html?=&t=o-que-e. Questionam-se temas específicos em trimestres distintos (entrevistas de 1 a 5), de modo que nem todas as variáveis estão em todas as entrevistas, conforme o tema.

# Temas e tópicos pesquisados ao longo do ano em determinada visita:
# Habitação (1a visita); 
# Características gerais dos moradores (1a visita); 
# Informações adicionais da força de trabalho (1a visita);
# Rendimentos de outras fontes (1a e 5a visitas). 
# UF Unidade da Federação 
# V2007 Sexo
# V2009 Idade do morador na data de referência
# V2010 Cor ou raça 
# V3007 já concluiu algum outro curso de graduação? 
# VD3004 Nível de instrução mais elevado alcançado (pessoas de 5 anos ou mais de idade) padronizado para o Ensino fundamental sistema de 9 anos
# VD4001 Condição em relação à força de trabalho na semana de referência para pessoas de 14 anos ou mais de idade 
# VD4002 Condição de ocupação na semana de referência para pessoas de 14 anos ou mais de idade # VD4020 Rendimento mensal efetivo de todos os trabalhos para pessoas de 14 anos ou mais de idade (apenas para pessoas que receberam em dinheiro, produtos ou mercadorias em qualquer trabalho) # VD4035 Horas efetivamente trabalhadas na semana de referência em todos os trabalhos para pessoas de 14 anos ou mais de idade

4 Puxando os dados

library(PNADcIBGE) httr::timeout(1e+12) # variaveis selecionadas variaveis_selecionadas <- c(“UF”, “V2007”, “V2009”, “V2010”, “V3007”, “VD3004”, “VD4001”, “VD4002”, “VD4020”, “VD4035”) dadosPNADc <- get_pnadc(year = 2019, interview = 1, vars = variaveis_selecionadas) class(dadosPNADc) # Salva objeto final saveRDS(dadosPNADc, “dadosPNADc.rds”)

5 Salvando os dados

library(readr)
dadosPNADc<- read_rds(path = "~/TEE 2020/Prova2/dadosPNADc.rds")
class(dadosPNADc)

[1] "survey.design2" "survey.design" 

library(survey)
variaveis_selecionadas <- c("UF", "V2007", "V2009", "V2010",
                            "V3007", "VD3004", "VD4001", "VD4002", "VD4020", "VD4035")

#Pode ser utilizada para a estimação do total de uma variável numérica, além da estimativa do total, o comando também retorna o erro padrão (SE) dessa estimativa.

# contagem por sexo
# 1 Homem e 2   Mulher
totalsexo <- svytotal(~V2007, dadosPNADc, na.rm = T)
totalsexo

                total     SE
V2007Homem  101073920 167095
V2007Mulher 108422543 167095

plot(totalsexo)

A opção na.rm = T, que remove as observações onde a variável é não-aplicável. O vetor v2007 utilizou pela svytotal valores binários para identificar as variaveis homem e mulher, fazendo a contagem do total no survey pela amostra da primeira entrevista da PNAD, a amostra demonstra que o numero de homens brasileiros são: 101073920 e mulheres são: 108422543

# contagem por sexo e raca sem interacao
# 1 Branca
# 2 Preta
# 3 Amarela
# 4 Parda 
# 5 Indígena
# 9 Ignorado

totalsexoraca <- svytotal(~V2007 + V2010, dadosPNADc, na.rm = T)
totalsexoraca

                  total       SE
V2007Homem    101073920 167095.4
V2007Mulher   108422543 167095.4
V2010Branca    89399699 488330.5
V2010Preta     19787999 243951.2
V2010Amarela    1395667  84689.2
V2010Parda     98121582 441478.4
V2010Indígena    765305  38645.6
V2010Ignorado     26211   8638.4

plot(totalsexoraca)

A contagem da amostra para grupos raciais sem interação. Em sua maioria são brancos: 89399699 e pardos: 98121582

# contagem por sexo, raca e com interacao
totalsexoEraca <- svytotal(~ interaction(V2007, V2010), 
                           dadosPNADc, na.rm = T)
ftable(totalsexoEraca)  # coluna A com valor e B com SE

                                                    A            B
interaction(V2007, V2010)Homem.Branca    42228563.133   258889.995
interaction(V2007, V2010)Mulher.Branca   47171136.144   287783.674
interaction(V2007, V2010)Homem.Preta      9822064.302   138013.525
interaction(V2007, V2010)Mulher.Preta     9965934.264   135156.874
interaction(V2007, V2010)Homem.Amarela     649805.359    47442.374
interaction(V2007, V2010)Mulher.Amarela    745861.666    46043.315
interaction(V2007, V2010)Homem.Parda     47981018.273   244559.350
interaction(V2007, V2010)Mulher.Parda    50140563.440   257646.403
interaction(V2007, V2010)Homem.Indígena    378649.756    22717.207
interaction(V2007, V2010)Mulher.Indígena   386655.577    22585.454
interaction(V2007, V2010)Homem.Ignorado     13819.218     6033.188
interaction(V2007, V2010)Mulher.Ignorado    12391.869     4340.246

plot(totalsexoEraca)

A interação serve para demonstrar os grupos raciais correlacionados com homens e mulheres. Em sua maioria, são homens e mulheres brancos, e, homens e mulheres pardos.

6 Estimando Médias

#Médias

# média da renda
mediarenda <- svymean(~VD4020, dadosPNADc, na.rm = T)
mediarenda

         mean     SE
VD4020 2386.1 33.662

cv(mediarenda)

           VD4020
VD4020 0.01410757

confint(mediarenda, level = .99)

          0.5 %   99.5 %
VD4020 2299.398 2472.814

A média de uma variável numérica é estimada através da função svymean, que possui uma sintaxe condizentemente idêntica ao svytotal. O exemplo do total da renda efetiva pode ser reescrito para médias: 2386.1 SE: 33.662. Observando estimativa do coeficiente de variação (CV) é 0.01410757, com intervalo de confiança de 99%. A renda média varia entre $2299,39 e $2472,81 no ano de 2019.

# proporção em variável categórica sexo
propsexo <- svymean(~V2007, dadosPNADc, na.rm = T)
propsexo

               mean    SE
V2007Homem  0.48246 8e-04
V2007Mulher 0.51754 8e-04

plot(propsexo)

# proporção em variável categórica sexo e raca sem interacao
propsexoraca <- svymean(~V2007 + V2010, dadosPNADc, na.rm = T)
propsexoraca

                    mean     SE
V2007Homem    0.48246122 0.0008
V2007Mulher   0.51753878 0.0008
V2010Branca   0.42673608 0.0023
V2010Preta    0.09445505 0.0012
V2010Amarela  0.00666201 0.0004
V2010Parda    0.46836868 0.0021
V2010Indígena 0.00365307 0.0002
V2010Ignorado 0.00012511 0.0000

plot(propsexoraca)

# proporção em variável categórica sexo e raca com cruzamento (interaction)
propsexoEraca <- svymean(~ interaction(V2007, V2010), dadosPNADc, na.rm = T)
ftable(propsexoEraca)

                                                    A            B
interaction(V2007, V2010)Homem.Branca    2.015717e-01 1.235773e-03
interaction(V2007, V2010)Mulher.Branca   2.251644e-01 1.373692e-03
interaction(V2007, V2010)Homem.Preta     4.688415e-02 6.587869e-04
interaction(V2007, V2010)Mulher.Preta    4.757090e-02 6.451511e-04
interaction(V2007, V2010)Homem.Amarela   3.101749e-03 2.264591e-04
interaction(V2007, V2010)Mulher.Amarela  3.560259e-03 2.197809e-04
interaction(V2007, V2010)Homem.Parda     2.290302e-01 1.167367e-03
interaction(V2007, V2010)Mulher.Parda    2.393385e-01 1.229837e-03
interaction(V2007, V2010)Homem.Indígena  1.807428e-03 1.084372e-04
interaction(V2007, V2010)Mulher.Indígena 1.845643e-03 1.078083e-04
interaction(V2007, V2010)Homem.Ignorado  6.596397e-05 2.879852e-05
interaction(V2007, V2010)Mulher.Ignorado 5.915073e-05 2.071752e-05

plot(propsexoEraca)

Utilizando variáveis categóricas, é possível estimar a frequência relativa de cada categoria. Isso pode ser feito também através da função svymean, com uma sintaxe análoga à utilizada para estimar os totais das categorias. Pelas amostras, é evidente e perceptível que os valores proporcioanais são proximos da realidade, pois a estimação de medidas de variabilidade (como a variância e o desvio padrão) para o cálculo dos intervalos de confiança e a realização de testes de hipóteses. As estimativas pontuais são associadas à precisão, dada pelo tamanho do intervalo de confiança. Quanto maior o intervalo, ou quanto maior o erro padrão, menor a precisão da estimativa. O SE em mairoia das variaveis tem valores contudentes,logo a baixo de 0,1%, ou seja, o nível de precisão da variaveis é alta.

7 Teste entre médias de renda para grupos de sexo

tt<-svyttest(VD4020 ~ V2007, dadosPNADc)
tt

    Design-based t-test

data:  VD4020 ~ V2007
t = -18.418, df = 11506, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -603.0318 -487.0295
sample estimates:
difference in mean 
         -545.0306 

confint(tt, level=0.99)

[1] -621.2617 -468.7996
attr(,"conf.level")
[1] 0.99

Os testes de hipóteses incluídos no pacote survey incluem teste-t para médias, teste qui-quadrado e teste de postos.O teste t para diferenças de médias de rendimentos efetivos entre sexos: Hipótese alternativa aceita, há diferenças entre rendas.

Também pode ser utilizado para estimar a renda média efetiva por unidade da federação

mediaRendaUF <- svyby(~VD4020, ~UF, dadosPNADc, svymean, na.rm = T)
mediaRendaUF

Sendo a renda efetiva mais alta é Distrito Federal $4235.22, porém seu SE é 258.40795, então os valores podem ser relativos, tendo diferanças do dado real a amostra complexa. A renda mais baixa é Maranhão, com $1322.92, seu SE é 53.75046.

confint(mediaRendaUF)

                       2.5 %   97.5 %
Rondônia            1859.002 2247.424
Acre                1722.538 2029.252
Amazonas            1520.849 1990.366
Roraima             1857.077 2481.246
Pará                1404.737 1639.242
Amapá               1646.825 2206.526
Tocantins           1672.235 2148.635
Maranhão            1217.573 1428.271
Piauí               1206.238 1535.658
Ceará               1472.307 1799.699
Rio Grande do Norte 1516.271 2103.177
Paraíba             1465.428 1912.317
Pernambuco          1556.890 2004.400
Alagoas             1379.361 1622.599
Sergipe             1409.613 1904.953
Bahia               1443.533 1727.837
Minas Gerais        1953.415 2161.383
Espírito Santo      2122.959 2519.779
Rio de Janeiro      2784.726 3137.633
São Paulo           2802.972 3280.010
Paraná              2454.935 2745.748
Santa Catarina      2543.014 2773.045
Rio Grande do Sul   2494.941 2800.663
Mato Grosso do Sul  2175.724 2552.926
Mato Grosso         2184.335 2483.387
Goiás               1985.374 2225.952
Distrito Federal    3728.751 4741.691

Intervalo de confiança de 95%

8 Estimando domínios

Muitas vezes queremos fazer uma análise para um domínio específico da população. Com o survey, esse domínio pode ser selecionado utilizando a função subset no objeto do plano amostral, aplicando a condição que define aquele domínio. É necessário que a seleção desse domínio nos dados seja feita somente após a criação do objeto que define o plano amostral, caso contrário podem ser obtidos resultados incorretos. Para a seleção dos domínios podem ser utilizadas condicionais com igualdade, como por exemplo para estimar a renda efetiva média de mulheres, outro exemplo para testar diferenças entre médias de horas trabalhadas.

mediarendaM <- svymean(~VD4020, subset(dadosPNADc, V2007 == "Mulher")  , na.rm = T)
mediarendaM

         mean     SE
VD4020 2076.8 29.009

mediarendaH <- svymean(~VD4020, subset(dadosPNADc, V2007 == "Homem")  , na.rm = T)
mediarendaH

         mean     SE
VD4020 2621.8 41.661

Pelos domínio é evidente a diferença entre médias salariais, sendo a média de renda efetiva da mulher: $2076.8 e homem: $2621.8, ou seja, em 2019, o homem aufere $545 reais a mais, do que a mulher.

#Por nivel de instrução (educação) Além disso, é possível utilizar múltiplas condições com os operados lógicos & (“e”) e |(“ou”), para estimar a frequência relativa de cada nível de instrução para mulheres pretas com mais de 25 anos, por exemplo.

nivelinstrHP25 <- svymean(~VD3004, subset(dadosPNADc, V2007 == "Mulher" & V2010 == "Preta" & V2009 > 25), na.rm = T)
nivelinstrHP25

                                                   mean     SE
VD3004Sem instrução e menos de 1 ano de estudo 0.082053 0.0032
VD3004Fundamental incompleto ou equivalente    0.321993 0.0056
VD3004Fundamental completo ou equivalente      0.075289 0.0030
VD3004Médio incompleto ou equivalente          0.052541 0.0026
VD3004Médio completo ou equivalente            0.300936 0.0056
VD3004Superior incompleto ou equivalente       0.040062 0.0029
VD3004Superior completo                        0.127125 0.0044

O estimador demonstra que mulheres pretas com 25 anos tem grau de instrução: Fundamental incompleto ou equivalente igual a 32,19%, sendo somente 12,71% das mulheres pretas tem superior completo.

Em outros casos há interesse em estimar quantidades de interesse para diversos domínios mutuamente exclusivos, a fim de possibilitar comparações. Para isto, podemos utilizar a função svyby.

# Se desejamos estimar a frequência relativa de homens e mulheres em cada nível de instrução, 
# usamos o seguinte código:

freqSexoInstr <- svyby(~V2007, ~VD3004, dadosPNADc, svymean, na.rm = T)
freqSexoInstr

Neste caso, foi utilizado para estimar frequência relativa de homens e mulheres em cada nível de instrução, na qual, SE esta abaixo de 1%. Logo, os valores identifcam mulheres superio completo são 59,48%, enquanto de homens é 40,51%, ou seja, quase 10% a mais de mulheres tem ensino superior, porém, aufere menos renda.

9 Estimando Razões

Além das proporções, também é possível estimar razões entre duas variáveis. Um exemplo de razão é a taxa de desocupação divulgada pelo IBGE: ela é a razão entre o total de pessoas desocupadas pelo total de pessoas na força de trabalho. A função para estimar razões é a svyratio. Sua sintaxe utiliza quatro argumentos: a variável cujo total estará no numerador, a variável cujo total estará no denominador, o objeto do plano amostral e a opção de remover valores não aplicáveis. Um exemplo para a estimativa dessa taxa é

# taxa de desocupação
txdesocup <- svyratio(~VD4002 == "Pessoas desocupadas",~VD4001 == "Pessoas na força de trabalho", dadosPNADc, na.rm = T)
txdesocup

Ratio estimator: svyratio.survey.design2(~VD4002 == "Pessoas desocupadas", ~VD4001 == 
    "Pessoas na força de trabalho", dadosPNADc, na.rm = T)
Ratios=
                                VD4001 == "Pessoas na força de trabalho"
VD4002 == "Pessoas desocupadas"                                0.1165897
SEs=
                                VD4001 == "Pessoas na força de trabalho"
VD4002 == "Pessoas desocupadas"                              0.001237404

O estimador relata que a taxa de desocupação é de 11,65% no Brasil em 2019

Cálculos de coeficiente de variação e intervalos de confiança para essa taxa:

cv(txdesocup)

                                VD4001 == "Pessoas na força de trabalho"
VD4002 == "Pessoas desocupadas"                               0.01061333

confint(txdesocup)

                                                                             2.5 %
VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.1141644
                                                                           97.5 %
VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.119015

O coeficiente de variação ficou em 1% e o intervalo de confiança de 95%, ajustou-se entre 11,41% e 11,19% de taxa de desocupação

Taxa de desocupação com desigualdade (V2009 idade >= 25 anos), condicionais com desigualdade, como por exemplo estimar a taxa de desocupação para pessoas com idade igual ou superior a 25 anos:

txdesocup25 <- svyratio(~VD4002 == "Pessoas desocupadas",~VD4001 == "Pessoas na força de trabalho", subset(dadosPNADc, V2009>=25) , na.rm = T)
txdesocup25

Ratio estimator: svyratio.survey.design2(~VD4002 == "Pessoas desocupadas", ~VD4001 == 
    "Pessoas na força de trabalho", subset(dadosPNADc, V2009 >= 
    25), na.rm = T)
Ratios=
                                VD4001 == "Pessoas na força de trabalho"
VD4002 == "Pessoas desocupadas"                               0.08491103
SEs=
                                VD4001 == "Pessoas na força de trabalho"
VD4002 == "Pessoas desocupadas"                              0.001112636

O estimador demonstrou que a taxa de desocupação de pessoas com 25 anos ou mais é de 8,4%

É possível definir domínios que sejam cruzamentos de variáveis categóricas com a função interaction. Na svyby também é possível utilizar vartype=“cv” caso desejemos que no output apareça o coeficiente de variação ao invés do erro padrão da estimativa. Para estimar a taxa de desocupação por sexo e cor/raça utilizando svyratioe svyby e o respectivo cv:

txdesocupSexoRaca <- svyby(~VD4002 == "Pessoas desocupadas", ~interaction(V2007,V2010), dadosPNADc, svyratio, denominator = ~VD4001 == "Pessoas na força de trabalho", na.rm = T, vartype = "cv")
txdesocupSexoRaca

De acordo com a taxa de desocupação por raça e genêro, a mulher indígena fica em primeiro lugar com 19,12% de desocupação, e, em seguida, a mulher preta em segundo lugar, com 18,25% de desocupação. O homem pardo é o que menor índice de desocupação, com 4,70%. Analisando veementemente, as mulheres índigenas e pretas sofrem mais descriminação.

#Test T para horas trabalhadas Testando diferenças entre médias de horas trabalhadas, entre concluintes e não concluintes de graduação:

svyttest(as.numeric(VD4035) ~ V3007, dadosPNADc)

    Design-based t-test

data:  as.numeric(VD4035) ~ V3007
t = -1.0594, df = 5209, p-value = 0.2895
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -2.0537607  0.6125467
sample estimates:
difference in mean 
         -0.720607 

Hipotese alternativa: tem diferenças.

10 Histograma

A função svyhist permite a criação de histogramas para variáveis numéricas que consideram os pesos amostrais para computar as frequências ou densidades das barras. A sintaxe básica é semelhante as de estimação apresentadas anteriormente. Além disto, possui os mesmos parâmetros que a função hist: breaks, xlab, main, etc… . Para estimativas das frequências absolutas, deve-se utilizar freq = TRUE. Abaixo, um exemplo do histograma do número de horas trabalhadas para cada um dos casos:

svyhist(~ as.numeric(VD4035), dadosPNADc, main = "Histograma", xlab = "Número de Horas Trabalhadas")

11 Boxplot

Para a construção de boxplots que considerem os pesos amostrais, a função é svyboxplot. A sintaxe dela difere um pouco das demais. É necessário declarar a variável numérica, sucedida por um ~ e a variável de agrupamento do boxplot. Caso não deseje usar grupos, basta colocar o número 1. Além disto, há a opção de plotar apenas os outliers mais extremos ou todos, através da opção all.outliers. Exemplos:

svyboxplot(VD4035 ~ V2007, dadosPNADc, main = "Boxplot do Número de Horas Trabalhadas por Sexo")

12 Modelos lineares

Além de testes de hipóteses, é possível estimar modelos lineares generalizados para dados amostrais complexos através da função svyglm. Para regressão linear simples ou múltipla, basta descrever a fórmula do modelo desejado na função svyglm. No exemplo abaixo, um modelo de regressão com rendimento efetivo como variável dependente e escolaridade, cor/raça e idade como variáveis explicativas. Para obter as estatísticas do modelo, pode ser utilizada as função summary, assim como feito para a regressão linear convencional.

modeloLin <- svyglm(VD4020 ~ VD3004 + V2010 + V2009, dadosPNADc)
summary(modeloLin)

Call:
svyglm(formula = VD4020 ~ VD3004 + V2010 + V2009, design = dadosPNADc)

Survey design:
survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)

Coefficients:
                                            Estimate Std. Error t value
(Intercept)                                 -809.079     83.050  -9.742
VD3004Fundamental incompleto ou equivalente  523.235     41.004  12.761
VD3004Fundamental completo ou equivalente    972.831     45.787  21.247
VD3004Médio incompleto ou equivalente       1203.046     56.892  21.146
VD3004Médio completo ou equivalente         1469.039     48.394  30.356
VD3004Superior incompleto ou equivalente    2083.572     86.041  24.216
VD3004Superior completo                     4714.437    113.937  41.377
V2010Preta                                  -669.560     40.297 -16.616
V2010Amarela                                 934.026    574.212   1.627
V2010Parda                                  -646.928     35.232 -18.362
V2010Indígena                               -494.432    150.009  -3.296
V2010Ignorado                               -796.822    876.541  -0.909
V2009                                         41.692      1.488  28.015
                                            Pr(>|t|)    
(Intercept)                                  < 2e-16 ***
VD3004Fundamental incompleto ou equivalente  < 2e-16 ***
VD3004Fundamental completo ou equivalente    < 2e-16 ***
VD3004Médio incompleto ou equivalente        < 2e-16 ***
VD3004Médio completo ou equivalente          < 2e-16 ***
VD3004Superior incompleto ou equivalente     < 2e-16 ***
VD3004Superior completo                      < 2e-16 ***
V2010Preta                                   < 2e-16 ***
V2010Amarela                                0.103845    
V2010Parda                                   < 2e-16 ***
V2010Indígena                               0.000984 ***
V2010Ignorado                               0.363341    
V2009                                        < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 6149977)

Number of Fisher Scoring iterations: 2

De acordo com a regressão linear,todos os anos estudados e raças explicam a varíavel: Rendimento mensal efetivo de todos os trabalhos para pessoas de 14 anos ou mais de idade (apenas para pessoas que receberam em dinheiro, produtos ou mercadorias em qualquer trabalho).

13 Regressão Logística

Outro modelo pertencente a classe de modelos generalizados é a Regressão Logística. Para utilizar este modelo, basta colocar o argumento family = “binomial”. A partir dos valores binários, podem ser melhor demonstrados os valores de saída para cada varíavel. No exemplo abaixo modelamos a não-conclusão de graduação pelo sexo, raça e idade

modeloLog <- svyglm(V3007 ~ VD3004 + V2010 + V2009, dadosPNADc, family = "binomial")
summary(modeloLog)

Call:
svyglm(formula = V3007 ~ VD3004 + V2010 + V2009, design = dadosPNADc, 
    family = "binomial")

Survey design:
survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)

Coefficients:
                                           Estimate Std. Error   t value
(Intercept)                               2.671e+01  4.197e-02   636.365
VD3004Superior incompleto ou equivalente -7.600e-03  2.579e-02    -0.295
VD3004Superior completo                  -5.341e+01  4.633e-02 -1152.941
V2010Preta                               -2.910e-02  4.336e-02    -0.671
V2010Amarela                             -4.728e-02  1.399e-01    -0.338
V2010Parda                               -6.353e-02  2.517e-02    -2.524
V2010Indígena                            -8.374e-02  1.783e-01    -0.470
V2010Ignorado                            -5.426e-02  8.816e-01    -0.062
V2009                                     8.415e-04  1.412e-03     0.596
                                         Pr(>|t|)    
(Intercept)                                <2e-16 ***
VD3004Superior incompleto ou equivalente   0.7682    
VD3004Superior completo                    <2e-16 ***
V2010Preta                                 0.5022    
V2010Amarela                               0.7354    
V2010Parda                                 0.0116 *  
V2010Indígena                              0.6387    
V2010Ignorado                              0.9509    
V2009                                      0.5513    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 3.609054e-14)

Number of Fisher Scoring iterations: 25

Nitidamente, a varíavel superior completo demonstra-se correlacionada com Rendimento mensal efetivo, ou seja, pela regressão logistica, é de fato, o individuo com maior anos de estudos, na qual, ensino superior completo, aufere maior renda. Em quesito de raça, o pardo teve p-valor de 1%.

14 ##Quarta questão: Painel:

Sejam os dados em painel enviados no e-mail, analisar os determinantes do crescimento do PIB municipal com os modelos por: Pooled OLS, Between, Fixed effects (within), First differences, e Random effects. São dados em painel para análise do Crescimento econômico e comércio externo nos municípios do estado de Mato Grosso do Sul entre 2003 e 2015, excluindo Figueirão e Paraíso das Águas. As variáveis foram deflacionadas para reais de 2015 e são listadas abaixo: • Taxa de Crescimento do Produto Interno Bruto Municipal (TXCRESC já com logs) • Produto Interno Bruto Municipal (PIBM) • Composição Setorial da Economia Municipal (SAGRO, SIND, SSERV e SADM já com logs) • Grau de Abertura Comercial Municipal (GA) • Imposto Líquido de Subsídios Municipal (IMPLIQ) • População Estimada Municipal (POPESTIM)

library(readr)
Painel<- read_rds(path = "~/TEE 2020/Prova2/dados3.rds")
class(Painel)

[1] "tbl_df"     "tbl"        "data.frame"

library(plm)
Painel.set <- pdata.frame(Painel, index = c("i", "t","ANO"))
pdim(Painel.set)

Balanced Panel: n = 77, T = 13, N = 1001

attach(Painel.set)

n = 77 (cross section), T = 13 (periodos de tempo), N = 1001 (observações)

15 Pooled

pool.1<-lm(txcresc ~ icms+sadm+sserv+sind+sagro+ypercap+impliq, data = Painel.set)
pool.2<-plm(txcresc ~ sadm+sserv+sind+sagro,data = Painel.set,
            model = "pooling")
library(stargazer)
stargazer(pool.1,pool.2, type = "text", title = "Pooled: Taxa de crescimento econômico" )

Pooled: Taxa de crescimento econômico
====================================================================
                                  Dependent variable:               
                    ------------------------------------------------
                                        txcresc                     
                              OLS                    panel          
                                                     linear         
                              (1)                     (2)           
--------------------------------------------------------------------
icms                         0.000                                  
                            (0.000)                                 
                                                                    
sadm                       -0.109***               -0.108***        
                            (0.026)                 (0.025)         
                                                                    
sserv                      0.091***                 0.091***        
                            (0.021)                 (0.020)         
                                                                    
sind                       0.104***                 0.105***        
                            (0.021)                 (0.021)         
                                                                    
sagro                      0.164***                 0.163***        
                            (0.009)                 (0.009)         
                                                                    
ypercap                    -0.00000                                 
                           (0.00000)                                
                                                                    
impliq                      -0.000                                  
                            (0.000)                                 
                                                                    
Constant                   0.036***                 0.035***        
                            (0.003)                 (0.002)         
                                                                    
--------------------------------------------------------------------
Observations                 1,001                   1,001          
R2                           0.341                   0.341          
Adjusted R2                  0.337                   0.339          
Residual Std. Error    0.044 (df = 993)                             
F Statistic         73.539*** (df = 7; 993) 128.931*** (df = 4; 996)
====================================================================
Note:                                    *p<0.1; **p<0.05; ***p<0.01

O primeiro modelo é o empilhado (pooled) a primeira tabela, as varíaveis: “ypercap” e “impliq” não foram significativas, sendo retiradas.

Teste Arellano

summary(pool.2,
        vcov = function(x) vcovHC(x, method = "arellano"))

Pooling Model

Note: Coefficient variance-covariance matrix supplied: function(x) vcovHC(x, method = "arellano")

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel.set, 
    model = "pooling")

Balanced Panel: n = 77, T = 13, N = 1001

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1124684 -0.0183431 -0.0043853  0.0083392  0.6175562 

Coefficients:
              Estimate Std. Error t-value  Pr(>|t|)    
(Intercept)  0.0353192  0.0039444  8.9543 < 2.2e-16 ***
sadm        -0.1077600  0.0404025 -2.6672 0.0077736 ** 
sserv        0.0908422  0.0248246  3.6594 0.0002661 ***
sind         0.1046267  0.0302309  3.4609 0.0005612 ***
sagro        0.1634202  0.0161186 10.1386 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.935
Residual Sum of Squares: 1.9337
R-Squared:      0.34115
Adj. R-Squared: 0.3385
F-statistic: 81.4853 on 4 and 76 DF, p-value: < 2.22e-16

Teste White

summary(pool.2,
        vcov = function(x) vcovHC(x, method = "white1"))

Pooling Model

Note: Coefficient variance-covariance matrix supplied: function(x) vcovHC(x, method = "white1")

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel.set, 
    model = "pooling")

Balanced Panel: n = 77, T = 13, N = 1001

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1124684 -0.0183431 -0.0043853  0.0083392  0.6175562 

Coefficients:
              Estimate Std. Error t-value  Pr(>|t|)    
(Intercept)  0.0353192  0.0034851 10.1345 < 2.2e-16 ***
sadm        -0.1077600  0.0302611 -3.5610 0.0003869 ***
sserv        0.0908422  0.0161833  5.6133 2.572e-08 ***
sind         0.1046267  0.0190369  5.4960 4.935e-08 ***
sagro        0.1634202  0.0160829 10.1611 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.935
Residual Sum of Squares: 1.9337
R-Squared:      0.34115
Adj. R-Squared: 0.3385
F-statistic: 118.391 on 4 and 76 DF, p-value: < 2.22e-16

16 First difference

firstdiff.1 <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel, model= "fd")
summary(firstdiff.1)

Oneway (individual) effect First-Difference Model

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel, 
    model = "fd")

Balanced Panel: n = 77, T = 13, N = 1001
Observations used in estimation: 924

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.3693565 -0.0088450  0.0028261  0.0135207  0.6539874 

Coefficients:
              Estimate Std. Error  t-value  Pr(>|t|)    
(Intercept) -0.0190844  0.0016701 -11.4269 < 2.2e-16 ***
sadm         0.4786508  0.0825123   5.8010 9.067e-09 ***
sserv        0.4050347  0.0378556  10.6995 < 2.2e-16 ***
sind         0.1665020  0.0356718   4.6676 3.500e-06 ***
sagro        0.2411619  0.0118826  20.2953 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.4796
Residual Sum of Squares: 1.3506
R-Squared:      0.45532
Adj. R-Squared: 0.45295
F-statistic: 192.054 on 4 and 919 DF, p-value: < 2.22e-16

17 Between

bet<- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set, model = "between")
summary(bet)

Oneway (individual) effect Between Model

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel.set, 
    model = "between")

Balanced Panel: n = 77, T = 13, N = 1001
Observations used in estimation: 77

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.0199496 -0.0069039 -0.0033663  0.0044138  0.0595513 

Coefficients:
             Estimate Std. Error t-value  Pr(>|t|)    
(Intercept) 0.0063090  0.0042439  1.4866  0.141489    
sadm        0.1032248  0.0362061  2.8510  0.005682 ** 
sserv       0.1830302  0.0291765  6.2732 2.345e-08 ***
sind        0.1505469  0.0259232  5.8074 1.590e-07 ***
sagro       0.1378101  0.0152975  9.0087 2.012e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    0.043937
Residual Sum of Squares: 0.010753
R-Squared:      0.75526
Adj. R-Squared: 0.74166
F-statistic: 55.547 on 4 and 72 DF, p-value: < 2.22e-16

18 Efeitos Fixos (Normais, individuais (one way), dois caminhos(two way) e no tempo(time))

fixed <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set,
             index=c("i", "t","ANO"), model="within")
###################################################################
fixed.ind <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set, 
                 index=c("i", "t","ANO"), effect = "individual",
                 model="within")
###################################################################

fixed.two <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set,
                 index=c("i", "t","ANO"), effect = "twoways",
                 model="within")
###################################################################
fixed.time <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set,
                  index=c("i", "t","ANO"), effect = "time",
                  model="within")
###################################################################
stargazer(fixed.ind,fixed.two,fixed.time, type = "text",
          column.labels = c("FE.Indiv", "FE.Two", "FE.Time"),
          align=TRUE,
          style = "all",
          keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

===============================================
                    Dependent variable:        
             ----------------------------------
                          txcresc              
              FE.Indiv     FE.Two     FE.Time  
                 (1)         (2)        (3)    
-----------------------------------------------
sadm          -0.243***    0.132**    0.107*** 
               (0.038)     (0.056)    (0.027)  
             t = -6.482   t = 2.362  t = 3.927 
              p = 0.000   p = 0.019  p = 0.0001
sserv         0.122***    0.195***    0.187*** 
               (0.029)     (0.028)    (0.019)  
              t = 4.160   t = 6.910  t = 10.019
             p = 0.00004  p = 0.000  p = 0.000 
sind            0.045     0.131***    0.145*** 
               (0.033)     (0.031)    (0.018)  
              t = 1.336   t = 4.175  t = 7.972 
              p = 0.182  p = 0.00004 p = 0.000 
sagro         0.187***    0.176***    0.157*** 
               (0.011)     (0.012)    (0.009)  
             t = 17.008  t = 14.116  t = 17.247
              p = 0.000   p = 0.000  p = 0.000 
-----------------------------------------------
Observations    1,001       1,001      1,001   
R2              0.305       0.265      0.383   
Adjusted R2     0.245       0.191      0.373   
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

Os valores foram próximos, com maior R(quadrado ajustado) o Efeito Fixo no tempo, porém somente os testes para averiguar melhor modelo.

#Fixed effects usando Least squares dummy variable model
fixed.dum <-lm(txcresc ~ sadm+sserv+sind+sagro+factor(i) - 1, data=Painel.set)
summary(fixed.dum)

Call:
lm(formula = txcresc ~ sadm + sserv + sind + sagro + factor(i) - 
    1, data = Painel.set)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.10224 -0.01848  0.00104  0.01259  0.55858 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
sadm        -0.243247   0.037525  -6.482 1.47e-10 ***
sserv        0.122079   0.029347   4.160 3.48e-05 ***
sind         0.044537   0.033335   1.336 0.181857    
sagro        0.187406   0.011019  17.008  < 2e-16 ***
factor(i)1   0.096259   0.011838   8.131 1.37e-15 ***
factor(i)2   0.040218   0.011941   3.368 0.000788 ***
factor(i)3   0.040270   0.011977   3.362 0.000804 ***
factor(i)4   0.066305   0.012147   5.459 6.17e-08 ***
factor(i)5   0.019718   0.011835   1.666 0.096043 .  
factor(i)6   0.039202   0.012328   3.180 0.001523 ** 
factor(i)7   0.060790   0.012128   5.012 6.45e-07 ***
factor(i)8   0.041225   0.012783   3.225 0.001304 ** 
factor(i)9   0.040431   0.011964   3.379 0.000757 ***
factor(i)10  0.031522   0.011933   2.642 0.008393 ** 
factor(i)11  0.015798   0.011806   1.338 0.181168    
factor(i)12  0.060547   0.013181   4.593 4.97e-06 ***
factor(i)13  0.069366   0.012376   5.605 2.75e-08 ***
factor(i)14  0.039707   0.012009   3.306 0.000982 ***
factor(i)15  0.053121   0.012399   4.284 2.03e-05 ***
factor(i)16  0.032564   0.012009   2.712 0.006823 ** 
 [ reached getOption("max.print") -- omitted 61 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.04225 on 920 degrees of freedom
Multiple R-squared:  0.6611,    Adjusted R-squared:  0.6313 
F-statistic: 22.16 on 81 and 920 DF,  p-value: < 2.2e-16

(fixed.dum$AIC <- AIC(fixed.dum))

[1] -3414.512

(fixed.dum$BIC <- BIC(fixed.dum))

[1] -3011.994

19 Efeitos aleatórios

random.set <- plm(txcresc ~ icms+sadm+sserv+sind+sagro+ypercap+impliq, 
                  data = Painel.set, model="random")
summary(random.set)

Oneway (individual) effect Random Effect Model 
   (Swamy-Arora's transformation)

Call:
plm(formula = txcresc ~ icms + sadm + sserv + sind + sagro + 
    ypercap + impliq, data = Painel.set, model = "random")

Balanced Panel: n = 77, T = 13, N = 1001

Effects:
                    var   std.dev share
idiosyncratic 1.779e-03 4.218e-02 0.994
individual    1.108e-05 3.329e-03 0.006
theta: 0.03817

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1125836 -0.0182762 -0.0042869  0.0087773  0.6155412 

Coefficients:
               Estimate  Std. Error z-value  Pr(>|z|)    
(Intercept)  3.6407e-02  2.8249e-03 12.8876 < 2.2e-16 ***
icms         6.6873e-12  1.2980e-11  0.5152    0.6064    
sadm        -1.1449e-01  2.6229e-02 -4.3649 1.272e-05 ***
sserv        9.0579e-02  2.1699e-02  4.1744 2.988e-05 ***
sind         1.0234e-01  2.1084e-02  4.8538 1.211e-06 ***
sagro        1.6529e-01  9.4586e-03 17.4747 < 2.2e-16 ***
ypercap     -1.7709e-08  5.6395e-08 -0.3140    0.7535    
impliq      -3.2536e-09  8.0966e-09 -0.4019    0.6878    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.8922
Residual Sum of Squares: 1.9145
R-Squared:      0.33806
Adj. R-Squared: 0.33339
Chisq: 507.136 on 7 DF, p-value: < 2.22e-16

Teste para fixed effects, hipótese: OLS (pool) é melho que fixed effects.Teste F para efeitos fixos (Chow test),‘pFtest’:

pFtest(fixed, pool.2)

    F test for individual effects

data:  txcresc ~ sadm + sserv + sind + sagro
F = 2.1502, df1 = 76, df2 = 920, p-value = 1.675e-07
alternative hypothesis: significant effects

pFtest(fixed.ind, pool.2) 

    F test for individual effects

data:  txcresc ~ sadm + sserv + sind + sagro
F = 2.1502, df1 = 76, df2 = 920, p-value = 1.675e-07
alternative hypothesis: significant effects

pFtest(fixed.two, pool.2)

    F test for twoways effects

data:  txcresc ~ sadm + sserv + sind + sagro
F = 5.3642, df1 = 88, df2 = 908, p-value < 2.2e-16
alternative hypothesis: significant effects

pFtest(fixed.time, pool.2)

    F test for time effects

data:  txcresc ~ sadm + sserv + sind + sagro
F = 29.669, df1 = 12, df2 = 984, p-value < 2.2e-16
alternative hypothesis: significant effects

O efeito fixo é melhor

O teste de Gourieroux et al. (1982) aparece numa alternativa:

g <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set, model = "pooling")
plmtest(g, effect = "twoways", type = "ghm")

    Lagrange Multiplier Test - two-ways effects (Gourieroux, Holly and
    Monfort) for balanced panels

data:  txcresc ~ sadm + sserv + sind + sagro
chibarsq = 1083.4, df0 = 0.00, df1 = 1.00, df2 = 2.00, w0 = 0.25, w1 =
0.50, w2 = 0.25, p-value < 2.2e-16
alternative hypothesis: significant effects

O efeito fixo é melhor

Teste de Hausman, H0: preferir random effects, H1: preferir efeitos fixos.

phtest(fixed, random.set)

    Hausman Test

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 152.94, df = 4, p-value < 2.2e-16
alternative hypothesis: one model is inconsistent

phtest(fixed.ind, random.set)

    Hausman Test

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 152.94, df = 4, p-value < 2.2e-16
alternative hypothesis: one model is inconsistent

phtest(fixed.two, random.set)

    Hausman Test

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 78.679, df = 4, p-value = 3.317e-16
alternative hypothesis: one model is inconsistent

phtest(fixed.time, random.set)

    Hausman Test

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 64.617, df = 4, p-value = 3.098e-13
alternative hypothesis: one model is inconsistent

O efeito fixo é melhor

Teste time-fixed effects. É nulo a função time-fixed effects se for:

#p-value  is < 0.05 then use . 
# F test for individual effects
#EFEITOS FIXOS NO TEMPO E RECOMENDAVEL
# Lagrange Multiplier Test - time effects (Breusch-Pagan)
plmtest(fixed, c("time"), type=("bp"))

    Lagrange Multiplier Test - time effects (Breusch-Pagan) for balanced
    panels

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 1066, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

plmtest(fixed.ind, c("time"), type=("bp"))

    Lagrange Multiplier Test - time effects (Breusch-Pagan) for balanced
    panels

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 1066, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

plmtest(fixed.two, c("time"), type=("bp"))

    Lagrange Multiplier Test - time effects (Breusch-Pagan) for balanced
    panels

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 1066, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

plmtest(fixed.time, c("time"), type=("bp"))

    Lagrange Multiplier Test - time effects (Breusch-Pagan) for balanced
    panels

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 1066, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

Os efeitos fixos no tempo foi rejeitado

Testing para correlação serial: Breusch–Godfrey teste para modelos de painel, H0: não há correlação serial.

pbgtest(fixed, data = Painel.set)

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 150.55, df = 13, p-value < 2.2e-16
alternative hypothesis: serial correlation in idiosyncratic errors

pbgtest(fixed.ind, data = Painel.set)

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 150.55, df = 13, p-value < 2.2e-16
alternative hypothesis: serial correlation in idiosyncratic errors

pbgtest(fixed.two, data = Painel.set)

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 55.132, df = 13, p-value = 3.832e-07
alternative hypothesis: serial correlation in idiosyncratic errors

pbgtest(fixed.time, data = Painel.set)

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  txcresc ~ sadm + sserv + sind + sagro
chisq = 61.175, df = 13, p-value = 3.235e-08
alternative hypothesis: serial correlation in idiosyncratic errors

Tem correlação serial nos erros idiosincráticos

Teste de Durbin Watson para correlação serial no modelo em diferenças, (H0: não tem correl.serial x H1: tem).

g2 <- plm(txcresc ~ sadm+sserv+sind+sagro+factor(ANO), data=Painel.set,
          index=c("i", "t","ANO"), model = "within")
pdwtest(g2)

    Durbin-Watson test for serial correlation in panel models

data:  txcresc ~ sadm + sserv + sind + sagro + factor(ANO)
DW = 2.0422, p-value = 0.628
alternative hypothesis: serial correlation in idiosyncratic errors

No modelo em diferenças, tem correlação diferencial.

Teste de raiz unitária, (H0: tem x H1: não tem).

library(tseries) 
adf.test(Painel.set$txcresc, k=2) 

    Augmented Dickey-Fuller Test

data:  Painel.set$txcresc
Dickey-Fuller = -11.768, Lag order = 2, p-value = 0.01
alternative hypothesis: stationary

Rejeita a 1%

Teste de heterocedasticidade, (H0: homocedástico x H1: heterocedástico)

library(lmtest) 
bptest(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set, studentize=F)

    Breusch-Pagan test

data:  txcresc ~ sadm + sserv + sind + sagro
BP = 989.92, df = 4, p-value < 2.2e-16

É homocedástico

Os testes indicaram problemas de correlacao serial e homocedasticidade.No caso,fazer o teste estimação de covariância de matrizes robustas (Sandwich estimator)

coeftest(fixed)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.037525 -6.4823 1.471e-10 ***
sserv  0.122079   0.029347  4.1599 3.482e-05 ***
sind   0.044537   0.033335  1.3361    0.1819    
sagro  0.187406   0.011019 17.0076 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fixed.ind)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.037525 -6.4823 1.471e-10 ***
sserv  0.122079   0.029347  4.1599 3.482e-05 ***
sind   0.044537   0.033335  1.3361    0.1819    
sagro  0.187406   0.011019 17.0076 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fixed.two)

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
sadm  0.132037   0.055893  2.3623   0.01837 *  
sserv 0.195444   0.028284  6.9100 9.117e-12 ***
sind  0.130894   0.031355  4.1745 3.273e-05 ***
sagro 0.175932   0.012463 14.1163 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fixed.time)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  0.1069961  0.0272443  3.9273 9.190e-05 ***
sserv 0.1867332  0.0186388 10.0185 < 2.2e-16 ***
sind  0.1452135  0.0182149  7.9722 4.301e-15 ***
sagro 0.1569249  0.0090988 17.2467 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No teste, efeitos fixos e efeitos fixos individuais, o p-valor da varíavel “sind” não foi significativa.

Segundo, Driscoll and Kraay (1998), teste estimação de covariância de matrizes robustas.

library(car)
zz <- plm(txcresc ~ sadm+sserv+sind+sagro, data=Painel.set,
          index=c("i", "t","ANO"), model="within")

#Teste de significância do coeficiente padrão

coeftest(zz)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.037525 -6.4823 1.471e-10 ***
sserv  0.122079   0.029347  4.1599 3.482e-05 ***
sind   0.044537   0.033335  1.3361    0.1819    
sagro  0.187406   0.011019 17.0076 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

SCC teste de significancia no modelo robusto padrão

coeftest(zz, vcov.=vcovSCC)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.125472 -1.9387   0.05285 .  
sserv  0.122079   0.050587  2.4133   0.01600 *  
sind   0.044537   0.038667  1.1518   0.24969    
sagro  0.187406   0.041385  4.5283 6.727e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Todos os parâmetros, passe vcov como uma função de argumento

coeftest(zz, vcov.=function(x) vcovSCC(x, type="HC1", maxlag=4))

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.098091 -2.4798   0.01332 *  
sserv  0.122079   0.050345  2.4248   0.01551 *  
sind   0.044537   0.032474  1.3715   0.17056    
sagro  0.187406   0.047561  3.9403 8.753e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(zz, vcov.=function(x) vcovSCC(x, type="HC1", maxlag=4))

Oneway (individual) effect Within Model

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel.set, 
    model = "within", index = c("i", "t", "ANO"))

Balanced Panel: n = 77, T = 13, N = 1001

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1022373 -0.0184782  0.0010377  0.0125917  0.5585828 

Coefficients:
       Estimate Std. Error t-value  Pr(>|t|)    
sadm  -0.243247   0.037525 -6.4823 1.471e-10 ***
sserv  0.122079   0.029347  4.1599 3.482e-05 ***
sind   0.044537   0.033335  1.3361    0.1819    
sagro  0.187406   0.011019 17.0076 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.3638
Residual Sum of Squares: 1.642
R-Squared:      0.30534
Adj. R-Squared: 0.24493
F-statistic: 101.097 on 4 and 920 DF, p-value: < 2.22e-16

#Alternativa de Croissant e Millo (2008,p.31)

coeftest(zz, vcov.=function(x) vcovHC(x, method="arellano", type="HC1"))

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
sadm  -0.243247   0.055723 -4.3653 1.413e-05 ***
sserv  0.122079   0.036994  3.3000  0.001004 ** 
sind   0.044537   0.049409  0.9014  0.367606    
sagro  0.187406   0.021548  8.6973 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(zz, vcov.=function(x) vcovHC(x, method="arellano", type="HC1"))

Oneway (individual) effect Within Model

Call:
plm(formula = txcresc ~ sadm + sserv + sind + sagro, data = Painel.set, 
    model = "within", index = c("i", "t", "ANO"))

Balanced Panel: n = 77, T = 13, N = 1001

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-0.1022373 -0.0184782  0.0010377  0.0125917  0.5585828 

Coefficients:
       Estimate Std. Error t-value  Pr(>|t|)    
sadm  -0.243247   0.037525 -6.4823 1.471e-10 ***
sserv  0.122079   0.029347  4.1599 3.482e-05 ***
sind   0.044537   0.033335  1.3361    0.1819    
sagro  0.187406   0.011019 17.0076 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    2.3638
Residual Sum of Squares: 1.642
R-Squared:      0.30534
Adj. R-Squared: 0.24493
F-statistic: 101.097 on 4 and 920 DF, p-value: < 2.22e-16

### metodologia parecida o teste estimação de covariância de matrizes robustas (Sandwich estimator)
fixef(zz,
      vcov.=function(x) vcovSCC(x, type="HC1", maxlag=4),
      effect = "individual")

         1          2          3          4          5          6          7 
 0.0962586  0.0402183  0.0402701  0.0663051  0.0197185  0.0392018  0.0607895 
         8          9         10         11         12         13         14 
 0.0412247  0.0404311  0.0315216  0.0157982  0.0605470  0.0693662  0.0397074 
        15         16         17         18         19         20         21 
 0.0531206  0.0325636 -0.0029739  0.0468802  0.0288873  0.0360804  0.0387891 
        22         23         24         25         26         27         28 
 0.0471172  0.0318366  0.0543707  0.0466447  0.0638597  0.0314187  0.0387586 
        29         30         31         32         33         34         36 
 0.0388266  0.0392987  0.0537103  0.0425143  0.0583311  0.0295258  0.0451973 
        37         38         39         40         41         42         43 
 0.0864612  0.0701419  0.0269606  0.0483509  0.0345045  0.0284094  0.0946750 
        44         45         46         47         48         49         50 
 0.0450376  0.0370874  0.0499149  0.0364127  0.0841676  0.0328454  0.0574399 
        51         52         53         54         55         56         57 
 0.0491154  0.0827274  0.0502011  0.0361257  0.0393278  0.0581578  0.0225664 
        59         60         61         62         63         64         65 
 0.0429644  0.1002649  0.0518472  0.0277103  0.0499841  0.0287218  0.0439645 
        66         67         68         69         70         71         72 
 0.0404121  0.0242084  0.0388206  0.0125797  0.0351352  0.0424167  0.0370882 
        73         74         75         76         77         78         79 
 0.0489222  0.0304575  0.0882792  0.0393123  0.0597596  0.0791390  0.0546583