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?