Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Eu sei que com o Python, você pode fazer praticamente tudo... Mas uma pessoa que não tem muito traquejo com computação, consegue rodar modelos de regressões em Python?

+1 voto
1,674 visitas
perguntada Fev 19, 2016 em Estatística por danielcajueiro (5,226 pontos)  
Compartilhe

1 Resposta

+1 voto
respondida Fev 19, 2016 por danielcajueiro (5,226 pontos)  

Sim. O Python "imitou" muito bem muitas ferramentas boas de outras linguagens incluindo o R. Já mencionamos aqui o Pandas uma ferramenta bastante adequada para lidar com base de dados inspirada nas estruturas de dados do R. Veja o exemplo abaixo como rodar regressões em Python.

Vamos usar a amostra de dados "MEAP93.raw" que é bem conhecida e discutida no livro "Introductory Econometrics - Jeffrey Wooldridge", que pode ser baixada no site da Cengage. Esse conjunto de dados pode ser usado para estudar o desempenho de estudantes em matemática usando a variável math10, que é a porcentagem de estudantes do décimo ano que tiram a nota mínima para passar no exame de matemática padrão americano, em função do tamanho da escola. Estamos interessados em estimar o seguinte modelo de regressão linear:

\[math10 = \beta_0 + \beta_1 \log(totcomp) + \beta_2 \log(staff) + \beta_3 \log(enroll),\]

onde totcomp é uma variável que mede o salário e benefícios dos professores e staff o número de funcionários por cada 1000 estudantes. A pergunta relevante e de interesse aqui é: Alunos tem melhores desempenho em escolas menores?

Rodar essa regressão em Python é super simples... A semenhança com o R provavelmente não deve ser mera coincidência... Veja no código abaixo, que já apresenta um sumário dos resultados:

# Load packages
import numpy as np
import statsmodels.api as sm
import pandas as pd
import patsy as ps


if __name__ == '__main__':
    # Load data    
    df = pd.read_csv('/home/daniel/Documents/Projetos/Prorum/Python For Econometrics/DataSets/Txt/MEAP93.raw',delim_whitespace=True,header=None)
    df.columns=['lnchprg','enroll','staff','expend','salary','benefits','droprate','gradrate','math10','sci11','totcomp','ltotcomp','lexpend','lenroll','lstaff','bensal','lsalary']

    # Multiple regression
    y,X = ps.dmatrices('math10 ~ ltotcomp + lstaff + lenroll',data=df, return_type='dataframe')
    model = sm.OLS(y,X) # Describe Model
    results = model.fit() # Fit model
    print results.summary()

Essa regressão apresenta os seguintes resultados:

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 math10   R-squared:                       0.065
Model:                            OLS   Adj. R-squared:                  0.058
Method:                 Least Squares   F-statistic:                     9.420
Date:                Fri, 19 Feb 2016   Prob (F-statistic):           4.97e-06
Time:                        17:17:38   Log-Likelihood:                -1523.7
No. Observations:                 408   AIC:                             3055.
Df Residuals:                     404   BIC:                             3072.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept   -207.6647     48.703     -4.264      0.000      -303.408  -111.922
ltotcomp      21.1550      4.056      5.216      0.000        13.182    29.128
lstaff         3.9800      4.190      0.950      0.343        -4.256    12.216
lenroll       -1.2680      0.693     -1.829      0.068        -2.631     0.095
==============================================================================
Omnibus:                       27.703   Durbin-Watson:                   1.666
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               34.495
Skew:                           0.568   Prob(JB):                     3.23e-08
Kurtosis:                       3.860   Cond. No.                     1.34e+03
==============================================================================

A regressão acima mostra que a variável totcomp é significante, ou seja, o salário dos professores tem efeito positivo no desempenho dos estudantes. Por outro lado, o número de funcionários não é significante.

Nesse problema específico estamos interessados em testar a hipótese nula \(H_0: \beta_{lenroll}\ge 0\) contra a hipótese alternativa \(H_a: \beta_{lenroll}\lt 0\). Para fazer isso, precisamos acrescentar algumas poucas linhas de código (ou usar a \(t\) estatística apresentada no sumário acima e a tabela de valores críticos). Para calcular apenas o p-valor, você precisa apenas de uma linha de programação:

pvalue=scipy.stats.t.cdf(results.tvalues['lenroll'],results.df_resid) 

Se você quiser a função completa para calcular qualquer caso e também os intervalos de confiança, você pode implementar o código abaixo:

def tTestOneTail(coefficient,referenceCoefficient,results,significanceLevel,oneCaudalTestType):
    coefficientValue=results.params[coefficient]

    standardErrorEstimator=results.bse[coefficient]
    degreesOfFreedom=results.df_resid

    tstat=(coefficientValue-referenceCoefficient)/standardErrorEstimator
    if(oneCaudalTestType=='smaller'):
        pvalue=scipy.stats.t.cdf(tstat,degreesOfFreedom)
    else:
        pvalue=1-scipy.stats.t.cdf(tstat,degreesOfFreedom)        

    conf1=coefficientValue-scipy.stats.t.ppf(1.0-significanceLevel, degreesOfFreedom)*standardErrorEstimator
    conf2=coefficientValue+scipy.stats.t.ppf(1.0-significanceLevel, degreesOfFreedom)*standardErrorEstimator    


    print "=================================================================="
    print "coef"+"     "+"std err"+"     "+"  t  "+"     "+"pvalue"+"     "+"[95%   conf. int]"
    print "{0:.4f}".format(coefficientValue),"  ","{0:.4f}".format(standardErrorEstimator),"   ","{0:.4f}".format(tstat),"   ","{0:.4f}".format(pvalue),"   ","{0:.4f}".format(conf1)," ","{0:.4f}".format(conf2)

Para chamar essa função, você precisa chamar as seguintes linhas que definem o coeficiente de interesse, o valor de referência do coeficiente, o nível de significância desejado e o tipo de teste unicaudal ('larger' ou 'smaller'):

coefficient='lenroll'
referenceCoefficient=0
significanceLevel=0.05
oneCaudalTestType='smaller' # Alternative Hypothesis
tTestOneTail(coefficient,referenceCoefficient,results,significanceLevel,oneCaudalTestType)

Essas linhas de código apenas recalculam os p-valores e os intervalos de confiança para testes de uma cauda e imprime os resultados:

coef     std err       t       pvalue     [95%   conf. int]
-1.2680    0.6932     -1.8292     0.0340     -2.4109   -0.1252

Esse resultado rejeita a hipótese nula em favor da hipótese alternativa que \(\beta_{enroll}\lt 0\). Ou seja, concluímos de acordo com esse modelo estatístico que escolas maiores minam o desempenho dos estudantes.

De fato, vários testes usuais de estatística ou econometria básico podem ser feitos em apenas um linha e de forma bem próxima a outros programas econométricos especializados em estatística e econometria. Por exemplo, imagine que você gostaria de testar a hipótese conjunta que \(\beta_{log(totcomp)}=0\) e \(\beta_{log(staff)}=0\), você precisa adicionar apenas as seguintes linhas para fazer o teste \(F\):

hypotheses = '(ltotcomp = 0),(lstaff=0)'
f_test = results.f_test(hypotheses)    
print f_test

Essas linhas apresentam o sumário do teste:

<F test: F=array([[ 13.62024751]]), p=1.88664201305e-06, df_denom=404, df_num=2>

onde o primeiro elemento é a estatística \(F\), o segundo é o \(p\)-valor e os últimos números são respectivamente os graus de liberdade do denominador da estatística e do numerador.

Se você tem interesse no assunto, veja esse caderno "Python for Introductory Econometrics - Wooldridge" que apresenta os exemplos do texto do Wooldridge em Python e inclui exemplos de OLS, WLS, FGLS, testes de hipótese, testes de adequação do modelo linear (como o de White) e vários outros testes estatísticos.

Se você deseja aumentar seu traquejo com programação, considere dar uma olhada nessas perguntas:

Como um programador iniciante pode se tornar um programador intermediário?

Como um programador em nível intermediário pode se tornar um programador em nível avançado?

Se você deseja aprender mais sobre python, dê uma olhada nessa pergunta:

Como iniciar em Python?

...