Пример построения и анализа регрессионной модели средствами SageMath

# read data
f = open('data.tsv', 'r')
head = f.readline() # skip the header
Y = []
X = []
for line in f :
    line = line.replace(',', '.').split() # replace commas by points and separate the cells
    row = [float(cell) for cell in line[1:]]
    Y.append(row[0])
    X.append([1] + row[1:]) # add a column of '1's
f.close()
Y = vector(RDF, Y)
X = matrix(RDF, X)
show(Y, X) # check if everything is correct

# computing 'beta_hat' coefficients
beta_hat = (X.transpose()*X)^(-1)*X.transpose()*Y
show(r'$\widehat{\beta}='+latex(beta_hat)+'$')

# computing 'u_hat' residuals
Y_hat = X*beta_hat # response (y) predicted by a model
u_hat = Y - Y_hat # just residuals
show(r'$\widehat{u}='+latex(u_hat)+'$')

# графическая проверка нормальности распределения остатков
a, b = min(u_hat), max(u_hat)
h = (b-a)/10
q = []
for i in xrange(10) :
    m = len([p for p in u_hat if a+i*h < p < a+(i+1)*h])
    q.append(float(m)/len(u_hat))
list_plot(q, plotjoined=True)
bar_chart(q)

print u_hat

# тест Шапиро - Уилка, если вы всё ещё не уверены :)
# UPD 16.04.17: p-value для этого теста должно быть БОЛЬШЕ 0.05 (нулевая гипотеза - распределение нормальное)
# Значит, тест показывает, что нормального распределения у нас нет, хотя было очень похоже (не дотянули совсем немного до 0.05)
%r
u = c(-0.5468890348882383, -0.6394825124745029, 1.1914692360582286, 0.029242763986664144, 0.8001097571719047, 3.105107767131045, -0.22834817910063343, 1.0849806355832943, 1.1862089694514628, -0.06558655308077999, -0.07230801203933712, -0.12910997423303794, -1.2140881897151417, -0.1916352895978788, -1.7888224594156386, -0.9023464418541947, -0.2799854876885579, 0.05888578451796178, 0.05464520892173752, 0.6828861286625809, -0.2819356441609493, 0.11921467245509199, -2.084379376131719, 0.025517761740225353, 0.08664846870165377)
shapiro.test(u)

# model quality - R^2, MSE, MAPE etc.
y_mean = mean(Y) # y average
T = X.nrows() # length of the sample
N = X.ncols() # number of predictors (with constant)
TSS = sum([(y-y_mean)^2 for y in Y]) # total sum of squares
ESS = sum([(yh-y_mean)^2 for yh in Y_hat]) # explained sum of squares
RSS = sum([uh^2 for uh in u_hat]) # residual sum of squares
R2 = ESS/TSS # R^2
show(r'$R^2='+latex(R2)+'$')
R2T = 1 - (1-R2)*(T-1)/(T-N) # adjusted R^2
show(r'$R^2_T='+latex(R2T)+'$')
R2A = 1 - (1-R2)*(T+N)/(T-N) # adjusted R^2
show(r'$R^2_A='+latex(R2A)+'$')
sigma_u2 = RSS/(T-N) # variance of residuals
show(r'$\widehat{\sigma}_u^2='+latex(sigma_u2)+'$')
MSE = RSS/T # MSE
show('$MSE='+latex(MSE)+'$')
MAPE = 100/T*sum([abs(u_hat[i])/Y[i] for i in range(T)]) # MAPE
show('$MAPE='+latex(MAPE)+'\,\%$')
F_stat = R2*(T-N)/(1-R2)/(N-1) # F-statistics
import scipy.stats
F_crit = scipy.stats.f.isf(0.1, N-1, T-N) # F-critical (or use F-distribution table directly)
if F_stat > F_crit :
    show('Model is significant')
else :
    show('Model is NOT significant')

# normalizing the data
Y_norm = vector([(y-y_mean)/std(list(Y)) for y in Y])
X_norm = matrix([[(X[i,j]-mean(X.column(j)))/std(list(X.column(j))) for j in xrange(1, N)] for i in xrange(T)])
# standartized equation
beta_hat_st = (X_norm.transpose()*X_norm)^(-1)*X_norm.transpose()*Y_norm
show(r'$\widehat{\beta}^{\text{st}}='+latex(beta_hat_st)+'$')

# elasticity
for k in xrange(1, N) :
    elast = beta_hat[k]*mean(X.column(k))/mean(Y)
    show(r'$\varepsilon_'+str(k+1)+'='+latex(elast)+'$')

# 'beta' variance-covariance matrix
Sigma_beta_hat = sigma_u2*(X.transpose()*X)^(-1)
show(r'$\widehat{\Sigma}_{\beta}='+latex(Sigma_beta_hat)+'$')

# 'beta' correlation matrix
r = matrix([[Sigma_beta_hat[i,j]/sqrt(Sigma_beta_hat[i,i])/sqrt(Sigma_beta_hat[j,j]) for j in xrange(N)] \
            for i in xrange(N)])
show('$r='+latex(r)+'$')

# t-test for 'beta' coefficients
for k in xrange(N) :
    t_stat = beta_hat[k]/sqrt(Sigma_beta_hat[k,k])
    import scipy.stats
    t_crit = scipy.stats.t.isf(0.1, T-N)
    if abs(t_stat) > t_crit :
        low = beta_hat[k] - t_crit*sqrt(Sigma_beta_hat[k,k])
        up = beta_hat[k] + t_crit*sqrt(Sigma_beta_hat[k,k])
        show(r'$\widehat{\beta}_'+str(k+1)+'$'+'is significant: '+\
             r'$\widehat{\beta}_'+str(k+1)+r'\in'+latex((low, up))+'$')
    else :
        show(r'$\widehat{\beta}_'+str(k+1)+'$'+'is NOT significant')

# prediction point
x_pred = vector([mean(x) for x in X.columns()])
Y_ppred = beta_hat * x_pred
show('Prediction point:', Y_ppred)

# prediction interval for mathematical expectation
low = Y_ppred - t_crit * sqrt(x_pred*(Sigma_beta_hat*x_pred))
up = Y_ppred + t_crit * sqrt(x_pred*(Sigma_beta_hat*x_pred))
show('Prediction interval for mathematical expectation: '+'$'+latex((low, up))+'$')

# prediction interval for individual value
low = Y_ppred - t_crit * sqrt(sigma_u2+x_pred*(Sigma_beta_hat*x_pred))
up = Y_ppred + t_crit * sqrt(sigma_u2+x_pred*(Sigma_beta_hat*x_pred))
show('Prediction interval for individual value: '+'$'+latex((low, up))+'$')

# partial determination coefficients
for k in xrange(1, N) :
    delta_R2 = (1-R2)/(T-N)*(beta_hat[k]/sqrt(Sigma_beta_hat[k,k]))^2
    show(r'$\Delta R^2_'+str(k+1)+'='+latex(delta_R2)+'$')

# model specification
# Ramsey test
RSS1 = RSS # RSS for model (1)
X_new = X
for p in xrange(2, 5) :
    X_new = X_new.augment(vector([yh^p for yh in Y_hat]))
beta_new = (X_new.transpose()*X_new)^(-1)*X_new.transpose()*Y
Y_hat_new = X_new * beta_new
RSS2 = sum([(Y[i]-Y_hat_new[i])^2 for i in xrange(T)]) # RSS for model (2)
F_stat = (RSS1-RSS2)*(T-3)/3/RSS2
F_crit = scipy.stats.f.isf(0.1, 3, T-3)
if F_stat > F_crit :
    show('Model is badly-specified')
else :
    show('Model is well-specified')
# Amemiya test
AF = RSS * (T+N) / (T-N)
show('$AF='+latex(AF)+'$')