Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
231 views
Tweet
Share
Share

Построение и анализ ЛММ с помощью R

Подготовка данных

# преобразуем файл к "читабельному виду"
fin = open('doctors.txt', 'r')
fout = open('doc.tsv', 'w')
data = fin.readlines()
fin.close()
head = [s[:-1] for s in data[:4]]
fout.write('\t'.join(head)+'\n')
for row in data[4:] :
    row = row.split()
    fout.write('\t'.join(row)+'\n')
fout.close()

%r
doctors <- read.table('doc.tsv', header=TRUE, sep='\t')
str(doctors) # yep!
'data.frame': 53 obs. of 4 variables: $ Index : int 1 2 3 4 5 6 7 8 9 10 ... $ Retardation.Index : num 2.8 3.1 2.59 3.36 2.8 3.35 2.99 2.99 2.92 3.23 ... $ Distrust.Index : num 6.1 5.1 6 6.9 7 5.6 6.3 7.2 6.9 6.5 ... $ Degree.of.Illness.after.6.Months: int 44 25 10 28 25 72 45 25 12 24 ...

Построение ЛММ

%r
y <- doctors$Degree.of.Illness.after.6.Months
x1 <- doctors$Retardation.Index
x2 <- doctors$Distrust.Index
doc.lm <- lm(y ~ x1 + x2)
summary(doc.lm)
Call: lm(formula = y ~ x1 + x2) Residuals: Min 1Q Median 3Q Max -22.166 -9.484 -2.016 6.412 40.366 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.8253 21.0440 -0.039 0.96887 x1 23.4913 6.8632 3.423 0.00124 ** x2 -7.0591 3.0218 -2.336 0.02354 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 14.91 on 50 degrees of freedom Multiple R-squared: 0.1989, Adjusted R-squared: 0.1668 F-statistic: 6.206 on 2 and 50 DF, p-value: 0.003913

Проверка нормальности

# проверим признаки на нормальность
%r
shapiro.test(y)
shapiro.test(x1)
shapiro.test(x2)
Shapiro-Wilk normality test data: y W = 0.84579, p-value = 7.176e-06 Shapiro-Wilk normality test data: x1 W = 0.96178, p-value = 0.08797 Shapiro-Wilk normality test data: x2 W = 0.95222, p-value = 0.03355 
# проверим остатки на нормальность
%r
u <- doc.lm$residuals # так можно "извлечь" сами остатки
summary(u) # можем убедиться, что M(u)=0 (Mean - среднее)
shapiro.test(u) # Ok
ks.test(u, "pnorm") # Ok
Min. 1st Qu. Median Mean 3rd Qu. Max. -22.170 -9.484 -2.016 0.000 6.412 40.370 Shapiro-Wilk normality test data: u W = 0.93689, p-value = 0.007625 One-sample Kolmogorov-Smirnov test data: u D = 0.49473, p-value = 1.678e-12 alternative hypothesis: two-sided 
# графическая проверка нормальности
%r
library(lattice)
histogram(u)

# и ещё один способ
%r
qqnorm(u)
qqline(u, col="red")

Оценка качества модели

# R-квадрат вычисляется сам по себе
# F-статистика вычисляется сама по себе
%r
yhat <- doc.lm$fitted.values # предсказанные значения отклика
sigma.u <- var(u)
print(paste("Дисперсия остатков: ", sigma.u))
mse <- mean(u^2)
print(paste("MSE =", mse))
mape <- mean(abs(u/y)*100)
print(paste("MAPE =", mape))
[1] "Дисперсия остатков: 213.789820845966" [1] "MSE = 209.756050641326" [1] "MAPE = 75.8067348907951"