[3] Simple Linear Regression

 

Visualization of multiple regression with two predictor variables

 

fi05 <- read.table("C:/Users/eznay/OneDrive/바탕 화면/회귀분석 데이터/CH06FI05.txt"); names(fi05) <- c("X1", "X2","Y"); library(car)

 

### 3D Scatter plot matrix in R : scatter3D( )  # 3D 산점도 그래프 생성 코드 

install.packages("plot3D") ; library(plot3D)

 

par(mfrow=c(2,2))

with(data=fi05, scatter3D(X1, Y, X2, phi=20, theta=20, xlab="X1", ylab="X2", zlab="Y", pch=20, col=1))

with(data=fi05, scatter3D(X1, Y, X2, phi=20, theta=40, xlab="X1", ylab="X2", zlab="Y", pch=20, col=1))

with(data=fi05, scatter3D(X1, Y, X2, phi=20, theta=60, xlab="X1", ylab="X2", zlab="Y", pch=20, col=1))

with(data=fi05, scatter3D(X1, Y, X2, phi=20, theta=80, xlab="X1", ylab="X2", zlab="Y", pch=20, col=1))

# fi05 <- read.table("CH06FI05.TXT"); names(fi05) <- c("X1", "X2","Y"): "CH06FI05.TXT" 파일을 읽어와서 데이터프레임 fi05를 생성하고 열 이름을 설정합니다. 이는 데이터를 불러오고 변수 이름을 설정하는 과정입니다. 데이터프레임에는 X1, X2, Y 열이 있다

# library(car): car 패키지를 로드합니다. 이 패키지는 산점도 행렬을 그리기 위해 사용

# install.packages("plot3D"); library(plot3D): plot3D 패키지를 설치하고 로드합니다. 이 패키지는 3D 산점도를 그리기 위해 사용

# par(mfrow=c(2,2)): 그래프를 2x2의 4개 패널로 나누어 출력하기 위해 그래프 레이아웃을 설정

# with(data=fi05, scatter3D(X1, Y, X2, phi=20, theta=20, xlab="X1", ylab="X2", zlab="Y", pch=20, col=1)): fi05 데이터프레임을 사용하여 3D 산점도를 그립니다. X1을 x축, Y를 y축, X2를 z축으로 설정하고, phi와 theta 매개변수를 이용하여 그래프의 시점과 각도를 조절합니다. xlab, ylab, zlab은 각각 x축, y축, z축의 레이블을 설정하며, pch와 col은 점의 모양과 색깔을 설정합니다. 이렇게 4개의 그래프를 그린다

 

 

### Scatter plot matrix in R : pairs( ) 

# 산점도 행렬 그리기 → 변수들 간 관계 확인 (fi05 변수들간의 모든 조합에 대한 산점도 생성)

pairs(fi05)

pairs(as.matrix(fi05))

Note that a matrix can be used as the argument of the function pairs

 

# pairs(fi05): fi05 데이터프레임을 사용하여 산점도 행렬을 그립니다. 각 변수들 간의 쌍별로 산점도를 그려줍니다. 예를 들어, X1과 X2의 관계, X1과 Y의 관계, X2와 Y의 관계를 확인

# pairs(as.matrix(fi05)): fi05 데이터프레임을 행렬로 변환한 후, 이를 사용하여 산점도 행렬을 그립니다. as.matrix 함수를 사용하여 데이터프레임을 행렬로 변환합니다. 이렇게 하면 데이터프레임의 열들이 변수로 간주되어 산점도 행렬을 그린다

---

Fitting multiple linear regression

### multiple linear regression fit with lm function 

# Im 함수를 사용해 다중 선형회귀 모델 적합, 결과 요약

→ 회귀 모델을 통해 독립 변수들이 종속 변수에 미치는 영향, 모델 유의성 평가 

fi05.lm <- lm(Y~X1+X2, data=fi05); summary(fi05.lm) ; anova(fi05.lm)

resd <- fi05.lm$resid; fitted <- fi05.lm$fitted

mse <- summary(fi05.lm)$sigma^2

 

[실습]

fi05.lm <- lm(Y ~ X1 + X2, data = fi05)
summary(fi05.lm)
anova(fi05.lm)

resd <- residuals(fi05.lm)
fitted <- fi05.lm$fitted.values
mse <- deviance(fi05.lm) / df.residual(fi05.lm)

 

# fi05.lm <- lm(Y~X1+X2, data=fi05): fi05 데이터프레임을 사용하여 종속 변수 Y와 독립 변수 X1, X2를 활용하여 다중 선형 회귀 모델을 적합합니다. lm 함수를 사용하여 회귀 모델을 적합하고, fi05.lm에 결과를 할당

# summary(fi05.lm): fi05.lm 모델의 요약 통계를 출력합니다. 회귀 계수, 회귀식의 설명력, 통계적 유의성 등을 확인

# anova(fi05.lm): fi05.lm 모델의 분산 분석(ANOVA) 결과를 출력합니다. 독립 변수들의 효과 및 모델 전체의 유의성을 평가하는 F-검정 결과가 제공

# resd <- fi05.lm$resid: 회귀 모델의 잔차(Residuals)를 resd에 할당

# fitted <- fi05.lm$fitted: 회귀 모델의 예측값(Fitted values)을 fitted에 할당

# mse <- summary(fi05.lm)$sigma^2: 회귀 모델의 평균 제곱 오차(Mean Squared Error, MSE)를 계산하여 mse에 할당합니다. summary 함수의 결과 중 잔차의 표준 편차를 제곱하여 MSE를 계산

 

## estimating beta through design matrix 

# 디자인 행렬을 사용해 베타 추정, 공분산 행렬 계산 → 베타의 표준 오차 산출 

Dx <- cbind(rep(1,21), fi05$X1, fi05$X2)

beta <- solve(t(Dx)%*%Dx)%*%t(Dx)%*%fi05$Y

cov.beta <- solve(t(Dx)%*%Dx)*mse sqrt(cov.beta[1,1]);

sqrt(cov.beta[2,2]); sqrt(cov.beta[3,3])

 

[실습]

cov.beta <- solve(t(Dx) %*% Dx) * mse
sqrt(cov.beta[1, 1])
sqrt(cov.beta[2, 2])
sqrt(cov.beta[3, 3])

# Dx <- cbind(rep(1,21), fi05$X1, fi05$X2): 디자인 행렬(Dx)을 생성합니다. fi05 데이터프레임의 X1과 X2 열을 이용하여 디자인 행렬을 구성합니다. 여기서 첫 번째 열은 1로 구성된 상수항

# beta <- solve(t(Dx)%*%Dx)%*%t(Dx)%*%fi05$Y: 최소제곱법을 사용하여 베타(beta)를 추정합니다. 디자인 행렬(Dx)과 종속 변수(fi05$Y)를 활용하여 베타를 추정합니다. solve 함수를 사용하여 역행렬을 계산하고, 행렬 연산을 수행하여 베타를 추정

# cov.beta <- solve(t(Dx)%*%Dx)*mse: 베타(beta)의 공분산 행렬(covariance matrix)을 계산합니다. 디자인 행렬(Dx)의 전치행렬과 자기 자신의 행렬 곱을 역행렬로 계산하고, 이를 MSE(Mean Squared Error)와 곱하여 공분산 행렬을 계산

# sqrt(cov.beta[1,1]): 베타(beta)의 첫 번째 요소의 표준 오차(standard error)를 계산합니다. 공분산 행렬의 첫 번째 요소를 제곱근하여 표준 오차를 계산합니다. 이는 상수항에 대한 베타의 표준 오차를 나타냄

# sqrt(cov.beta[2,2]): 베타(beta)의 두 번째 요소의 표준 오차를 계산합니다. 공분산 행렬의 두 번째 요소를 제곱근하여 표준 오차를 계산합니다. 이는 X1에 대한 베타의 표준 오차를 나타냄

# sqrt(cov.beta[3,3]): 베타(beta)의 세 번째 요소의 표준 오차를 계산합니다. 공분산 행렬의 세 번째 요소를 제곱근하여 표준 오차를 계산합니다. 이는 X2에 대한 베타의 표준 오차를 나타냄

 

### Model diagnosis # 다중 회귀 모델 진단 

par(mfrow=c(2,2)); plot(fi05.lm)

ncvTest(fi05.lm); bc <- powerTransform(fi05.lm); summary(bc); boxCox(fi05.lm)

# par(mfrow=c(2,2)); plot(fi05.lm): 다중 회귀 모델(fi05.lm)의 진단을 위해 그래프를 그립니다. par 함수를 사용하여 그래프를 2x2의 모양으로 나열하고, plot 함수를 사용하여 모델의 진단 그래프를 생성

# ncvTest(fi05.lm): 다중 회귀 모델(fi05.lm)의 비선형성을 검정하는 Brown-Forsythe 검정을 수행합니다. ncvTest 함수를 사용하여 검정을 수행하고, 모델의 비선형성 여부를 평가

# bc <- powerTransform(fi05.lm); summary(bc); boxCox(fi05.lm): 다중 회귀 모델(fi05.lm)의 종속 변수에 대한 Box-Cox 변환을 수행합니다. powerTransform 함수를 사용하여 Box-Cox 변환을 수행하고, summary 함수를 통해 변환 결과를 요약합니다. 또한, boxCox 함수를 사용하여 Box-Cox 변환의 람다(lambda) 값을 확인

 

### Interval estimation and prediction for multiple regression

## Estimation of all data # 예측  

esti.95<-cbind(fi05, predict.lm(fi05.lm, newdata=fi05, se.fit=TRUE, interval="confidence"))

 

X.h <- data.frame(1, X1=65.4, X2=17.6)

# esti.95<-cbind(fi05, predict.lm(fi05.lm, newdata=fi05, se.fit=TRUE, interval="confidence")): 다중 회귀 모델(fi05.lm)을 사용하여 모든 데이터에 대한 추정값과 신뢰 구간을 계산합니다. predict.lm 함수를 사용하여 새로운 데이터(fi05)에 대한 추정값과 신뢰 구간을 계산하고, cbind 함수를 사용하여 원래 데이터(fi05)와 결합

# X.h <- data.frame(1, X1=65.4, X2=17.6): 다중 회귀 모델(fi05.lm)을 사용하여 새로운 데이터(X1=65.4, X2=17.6)에 대한 예측값을 계산하기 위해 새로운 데이터 프레임(X.h)을 생성

 

## Estimation for X.h # 추정값, 신뢰구간 계산 

predict.lm(fi05.lm, newdata=X.h, se.fit=TRUE, interval="confidence")

 

## by directly using the estimated covariance matrix

se.X.h <- c(sqrt (as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h)) ))

EY <- c(as.matrix(X.h)%*% fi05.lm$coef)

c(EY- qt(0.975, 18)*se.X.h, EY+qt(0.975, 18)*se.X.h)

# predict.lm(fi05.lm, newdata=X.h, se.fit=TRUE, interval="confidence"): 다중 회귀 모델(fi05.lm)을 사용하여 새로운 데이터(X.h)에 대한 추정값과 신뢰 구간을 계산합니다. predict.lm 함수를 사용하여 새로운 데이터에 대한 추정값과 신뢰 구간을 계산합니다. se.fit=TRUE는 추정값의 표준 오차를 함께 계산하고, interval="confidence"는 신뢰 구간을 계산하도록 지

# se.X.h <- c(sqrt (as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h)) )): 추정된 공분산 행렬을 직접 사용하여 X.h에 대한 표준 오차를 계산합니다. as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h))를 통해 X.h에 대한 공분산 행렬을 계산하고, sqrt 함수를 사용하여 표준 오차를 계산

# EY <- c(as.matrix(X.h)%*% fi05.lm$coef): X.h에 대한 예측값을 계산하기 위해 X.h와 회귀 모델의 계수를 곱하여 예측값을 계산

# c(EY- qt(0.975, 18)*se.X.h, EY+qt(0.975, 18)*se.X.h): X.h에 대한 예측값과 신뢰 구간을 계산합니다. EY에서 qt(0.975, 18)*se.X.h를 빼고 더하여 신뢰 구간을 계산합니다. qt(0.975, 18)는 자유도 18인 t-분포에서 상위 2.5%에 해당하는 값으로, 95% 신뢰 수준

 

## Prediction for X.h # 예측값, 예측 구간 계산 

predict.lm(fi05.lm, newdata=X.h, se.fit=TRUE, interval="prediction")

 

## by directly using the estimated covariance matrix

se.X.h.pred <- c(sqrt (mse+as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h)) ))

Y.pred <- c(as.matrix(X.h)%*% fi05.lm$coef)

c(Y.pred- qt(0.975, 18)*se.X.h.pred, EY+qt(0.975, 18)*se.X.h.pred)

# predict.lm(fi05.lm, newdata=X.h, se.fit=TRUE, interval="prediction"): 다중 회귀 모델(fi05.lm)을 사용하여 새로운 데이터(X.h)에 대한 예측값과 예측 구간을 계산합니다. predict.lm 함수를 사용하여 새로운 데이터에 대한 예측값과 예측 구간을 계산합니다. se.fit=TRUE는 예측값의 표준 오차를 함께 계산하고, interval="prediction"은 예측 구간을 계산하도록 지정

# se.X.h.pred <- c(sqrt (mse+as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h)) )): 추정된 공분산 행렬을 직접 사용하여 X.h에 대한 예측 구간의 표준 오차를 계산합니다. mse에 as.matrix(X.h)%*%cov.beta%*% t(as.matrix(X.h))를 더하여 예측 구간의 표준 오차를 계산

# Y.pred <- c(as.matrix(X.h)%*% fi05.lm$coef): X.h에 대한 예측값을 계산하기 위해 X.h와 회귀 모델의 계수를 곱하여 예측값을 계산

# c(Y.pred- qt(0.975, 18)*se.X.h.pred, EY+qt(0.975, 18)*se.X.h.pred): X.h에 대한 예측값과 예측 구간을 계산합니다. Y.pred에서 qt(0.975, 18)*se.X.h.pred를 빼고 더하여 예측 구간을 계산합니다. qt(0.975, 18)는 자유도 18인 t-분포에서 상위 2.5%에 해당하는 값으로, 95% 예측 구간을 가진다

---

General Linear Test # 일반 선형 검정 수행 

ta701<-read.table("C:/Users/eznay/OneDrive/바탕 화면/회귀분석 데이터/CH07TA01.txt"); names(ta701) <- c("X1","X2", "X3", "Y")

 

## testing beta1=0 with predictors (X1, X2, X3)

H0.1<- lm(Y~X2+X3, data=ta701); H1.1<- lm(Y~X1+X2+X3, data=ta701); anova(H0.1, H1.1)

# H0.1 <- lm(Y~X2+X3, data=ta701); H1.1 <- lm(Y~X1+X2+X3, data=ta701); anova(H0.1, H1.1): 예측 변수(X1, X2, X3) 중에서 beta1=0인지를 검정합니다. H0.1은 beta1=0인 모델을, H1.1은 beta1이 0이 아닌 모델을 나타냅니다. anova 함수를 사용하여 두 모델 간의 F-통계량과 p-값을 계산

 

## testing beta1=beta2=beta3

X123 <- ta701$X1+ta701$X2+ta701$X3

H0.2<- lm(Y~X123, data=ta701); H1.2<- lm(Y~X1+X2+X3, data=ta701); anova(H0.2, H1.2)

# X123 <- ta701$X1+ta701$X2+ta701$X3; H0.2 <- lm(Y~X123, data=ta701); H1.2 <- lm(Y~X1+X2+X3, data=ta701); anova(H0.2, H1.2): 예측 변수(X1, X2, X3)들의 계수들이 모두 같은지를 검정합니다. X123는 X1, X2, X3의 합을 나타내는 변수입니다. H0.2는 beta1=beta2=beta3인 모델을, H1.2는 적어도 하나의 beta가 다른 모델을 나타냅니다. anova 함수를 사용하여 두 모델 간의 F-통계량과 p-값을 계산

 

## testing beta1=beta2 with predictors (X1, X2, X3)

H0.3<- lm(Y~X12+X3, data=ta701); H1.3<- lm(Y~X1+X2+X3, data=ta701); anova(H0.3, H1.3)

 

[실습]

H0.3 <- lm(Y ~ X1:X2 + X3, data = ta701)
H1.3 <- lm(Y ~ X1 + X2 + X3, data = ta701)
anova(H0.3, H1.3)

# H0.3 <- lm(Y~X12+X3, data=ta701); H1.3 <- lm(Y~X1+X2+X3, data=ta701); anova(H0.3, H1.3): 예측 변수(X1, X2, X3) 중에서 beta1=beta2인지를 검정합니다. H0.3은 beta1=beta2인 모델을, H1.3은 beta1과 beta2가 다른 모델을 나타냅니다. anova 함수를 사용하여 두 모델 간의 F-통계량과 p-값을 계산

 

## testing beta1=beta2=0 with predictors (X1, X2, X3)

H0.4<- lm(Y~X3, data=ta701); H1.4<- lm(Y~X1+X2+X3, data=ta701); anova(H0.4, H1.4)

# H0.4 <- lm(Y~X3, data=ta701); H1.4 <- lm(Y~X1+X2+X3, data=ta701); anova(H0.4, H1.4): 예측 변수(X1, X2, X3) 중에서 beta1=beta2=0인지를 검정합니다. H0.4는 beta1=beta2=0인 모델을, H1.4는 적어도 하나의 beta가 0이 아닌 모델을 나타냅니다. anova 함수를 사용하여 두 모델 간의 F-통계량과 p-값을 계산

Standardized Regression # 변수 표준화 

## standardizing all variables

ta701.s <- as.data.frame(scale(ta701))

## checking zero means and 1 standard deviation

apply(ta701.s, 2, mean) ; apply(ta701.s, 2, sd) ;

## fitting

stan.reg <- lm(Y~X1+X2+X3, data=ta701.s); summary(stan.reg)

# ta701.s <- as.data.frame(scale(ta701)): scale 함수를 사용하여 데이터셋 ta701의 모든 변수를 표준화합니다. 결과는 ta701.s라는 새로운 데이터프레임으로 저장

# apply(ta701.s, 2, mean): 표준화된 데이터셋 ta701.s의 각 변수에 대해 평균을 계산합니다. apply 함수를 사용하여 각 열(axis=2)에 대해 평균을 계산합니다. 이를 통해 각 변수의 평균이 거의 0에 가까운지 확인

# apply(ta701.s, 2, sd): 표준화된 데이터셋 ta701.s의 각 변수에 대해 표준 편차를 계산합니다. apply 함수를 사용하여 각 열(axis=2)에 대해 표준 편차를 계산합니다. 이를 통해 각 변수의 표준 편차가 거의 1에 가까운지 확인

# stan.reg <- lm(Y~X1+X2+X3, data=ta701.s): 표준화된 데이터셋 ta701.s를 사용하여 선형 회귀 모델을 적합시킵니다. 종속 변수 Y와 독립 변수 X1, X2, X3 사이의 관계를 모델링하는 선형 회귀 모델이 생성됩니다. 결과는 stan.reg에 저장되며, summary(stan.reg)를 통해 모델의 요약 정보를 확인