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

f = open('data.tsv', 'r')
s = f.read().replace(',', '.')
f.close()
fnew = open('d.tsv', 'w')
fnew.write(s)
fnew.close()

%r
d <- read.table('d.tsv', header=TRUE, sep='\t')
str(d) # Ok!

%r
d.lm <- lm(Y ~ X2 + X3 + X4 + X5, data=d)
summary(d.lm)
yhat <- d.lm$fitted.values
u <- d.lm$residuals
mse <- mean(u^2)
mape <- 100*mean(abs(u/d$Y))
print(paste("MSE =", mse))
print(paste("MAPE =", mape))

%r
par(mfrow = c(2, 2))
plot(d.lm)
par(mfrow = c(1, 1))

%r
u <- d.lm$residuals
# 1
shapiro.test(u)
# 2
library(lattice)
histogram(u)
# 3
qqnorm(u)
qqline(u, col="green")

%r
# тест Рэмси
library(lmtest)
resettest(d.lm, power=2:4) # большое p-значение --> хорошо!
# кол-во регрессоров + константа
N = ncol(d) - 1
# длина выборки
T = nrow(d)
# функция Амемии
af <- sum(u^2) * (T+N) / (T-N)
print(paste("AF =", af))

%r
sd.lm <- lm(Y ~ X2 + X3 + X4 + X5, data=data.frame(scale(d)))
summary(sd.lm)

%r
eps <- 0
beta <- d.lm$coefficients
for (k in 2:5) {
    eps[k] <- beta[k] * mean(d[,k+1]) / mean(d[,2])
}
eps

%r
sigma.beta <- summary(d.lm)$cov.unscaled
sigma.beta

%r
r.beta <- summary(d.lm, correlation=TRUE)$correlation
r.beta

# Значимость проверяется автоматически; см. результаты функции lm()
%r
# резюме по коэффициентам
coef(summary(d.lm))
# доверительные интервалы
confint(d.lm, level=0.9)

%r
beta <- d.lm$coefficients
R2 <- summary(d.lm)$r.squared
delta.R2 <- 0
for (k in 2:5) {
    delta.R2[k] <- ((beta[k]/sqrt(sigma.beta[k,k]))^2)*(1-R2)/(T-N)
}
delta.R2

# выберем для построения модели 20 случайных наблюдений
%r
train <- sample(1:25, 20)
d.train <- lm(Y ~ X2 + X3 + X4 + X5, data=d, subset=train)
summary(d.train)

%r
new <- d[-train,]
d.predict <- predict(d.train, new)
d.predict
d$Y[-train]
predict(d.train, new, interval="prediction")
predict(d.train, new, interval="confidence")

# прогнозы для средних
%r
d[26,] <- c(26, 0, mean(d$X2), mean(d$X3), mean(d$X4), mean(d$X5))
new <- d[26,]
predict(d.lm, new)
predict(d.lm, new, interval="prediction")
predict(d.lm, new, interval="confidence")