R- 시계열 모형

1. 정상성

– 평균이 일정하다. 즉 모든 시점에 대해 일정한 평균을 가진다.

– 분산도 시점에 의존하지 않는다

– 공분산은 단지 시차에만 의존하고 실제 어느 시점 t,s 에는 의존하지 않는다.

대부분의 경우 실제의 데이터는 이러한 정상성을 만족하지 않는다.

그래서 정상성을 만족하지 않는 데이터를 정상 시계열 데이터로 변환해야 하는데

이때, 변환, 차분과 같은 방법을 사용할 수 있다.

일반적으로 평균이 일정하지 않는 경우에는 차분을 이용하고, 분산이 일정하지 않은 경우에는

변환을 이용하여 정상 시계열 데이터로 변환한다.

차분의 경우 바로 전 시점의 데이터를 빼는 것을 일반차분, 계절성을 갖는 데이터의 경우 여러 시점전의

데이터를 빼는 계절차분을 사용한다.

2. 시계열 모형

– 자기 회귀 모형 (AR)

: P 시점전의 자료가 현재 자료에 영향을 주는 모형

Z(현재 시점) = Z (현재-(P+1)) * (영향도 계수) + Z (현재-(P+2)) * (영향도 계수) + …. + 백색잡음

–> AR(1) : Z = @Z(t-1)

AR(2) : Z= @Z(t-1) + @Z(t-2)

– 이동평균 모형 (MA)

: 유한한 개수의 백색잡음의 결합으로 언제나 정상성을 만족한다

그전 시간의 변수가 아니라 백색잡음만 가지고 식을 만든다…

(백색 잡음이란 시계열 분석에서의 오차항을 의미함 )

–> MA(1) : Z = a(핸재잡음) – a(t-1)(현재시간 -1 잡읍)

MA(2) : Z = a(핸재잡음) – a(t-1)(현재시간 -1 잡읍) – a(t-2)(현재시간 -2 잡읍)

– 자기회귀누적이동평균모형 (ARIMA)

: 비정상 시계열 모델 ARIMA 를 정상 시계열 ARMA 로 만들어야 하는데 …

차분 몇번 변환 몇번을 하였는가에 따라 ARIMA(p, d, q) 가 달라진다..

ARIMA(p , d, q)

-> p : AR , q : MA, d : 차분 횟수

d(차분횟수)가 0 이면 ARMA(p,q) 이되고 정상성을 만족함

p(AR) 0이면 IMA(d,q) 모형으 따르고 d 번 차분하면 MA(q) 로 정상성을 만족

나머지 케이스도 동일 ~

3. 분해 시계열

– 시계열에 영향을 주는 일반적인 요인을 시계열에서 분리해 분석하는 방법을 말하면 회귀적인 방법을 주로사용

1) 추세요인

: 데이터 오르거나 내리거나 2차식 형태를 취하거나 하는 형태가 있을 경우~ 추세가 있다고 한다.

2) 계절요인

: 고정주기에 자료가 변화할 경우 (계절처럼 고정 주기가 있는 경우) 계절용인이 있다고 한다.

3) 순환요인

: 명백한 경제적이거나 자연적인 이유 없이 알려지지 않은 주기를 가지고 변화하는 자료 (여튼 주기가 있다!)

4) 불규칙요인

: 위의 세가지 요인으로 설명할수 없는 오차에 해당하는 데이터를 불규칙 요인이라고 한다.

※ 이제 코드를 좀 봐야 할 시간
install.packages(“xts”)
install.packages(“TTR”)
install.packages(“forecast”)
install.packages(“zoo”)
install.packages(“timeDate”)
install.packages(“colorspace”)
install.packages(“fracdiff”)
install.packages(“tseries”)
install.packages(“Rcpp”)
library(TTR)
library(forecast)

R – Decision Tree

install.packages(“rattle”)
install.packages(“RColorBrewer”)
install.packages(“rpart.plot”)

#######################################
#1. 데이터 구성
#######################################
library(rattle)
library(RColorBrewer)
library(rpart.plot)

#일정한 Random 값을 갖도록한다. (언제 실행해도 1426 이라면 같은 Random 값)
set.seed(1426)

#R 에 내장된 샘플 CSV weather 를 읽어 온다
data(weather)

#weather 데이터를 ds 에 담는다.
ds <- get(“weather”)

# 제외할 컬럼 “Date” “Location” “RISK_MM”
id <- c(“Date”, “Location”)
risk <- “RISK_MM”
ignore <- c(id, if (exists(“risk”)) risk )

# vars : RainTommmorw 포함, 제외 컬럼은 제거
(vars <- setdiff(names(ds), ignore)) #B에는 없고 A 에만 있는 값만 리턴

# vars : RainTommmorw 제거
target <- “RainTomorrow”
inputs <- setdiff(vars, target)

# 원본 데이터의 Row 수 366
(nobs <- nrow(ds))

#intersect : 두 데이터가 공통적으로 가지고 있는 데이터만 리턴
#inputs : RainTomorrow 를 제외한 컬럼 리스트
#names(ds) : 데이터의 전체 컬럼 리스트
#vars : RainTomorrow 를 포함하고 제외항목을 뺀 리스트
#ds[vars] : vars 컬럼과 일치하는 데이터 셋
#apply : apply(M, 1, min)
# – M : 데이터
# – 1 : 열 , 2 : 행
# – 처리할 함수 : min 최소, max 최대, mean 평균
#sapply : apply 와 비슷함. 행으로 리턴
#lapply : apply 와 비슷함. 열로 리턴
#names(ds)[which(sapply(ds[vars], is.numeric))] : 값이 숫자인 컬럼만 리턴
(numerics <- intersect(inputs, names(ds)[which(sapply(ds[vars], is.numeric))]))

#데이터가 Fator (Char 포함) / RainTomorrow 제외
(categorics <- intersect(inputs, names(ds)[which(sapply(ds[vars], is.factor))]))

#RainTomorrow ~ . 이라고 formula 생성 / 막다른 변수 같은거 넣어서도 만들 수 있음
(form <- formula(paste(target, “~ .”)))

## nobs(366) 안에서 366*0.7 만큰의 Rnadom 값을 생성
length(train <- sample(nobs, 0.7*nobs))
length(test <- setdiff(seq_len(nobs), train))

## 내일 비가 올 것이냐, 그 확률 두가지 데이터를 뽑았다.
actual <- ds[test, target]
risks <- ds[test, risk]
#######################################
#2. 데이터 확인
#######################################
dim(ds)
names(ds)
head(ds)
tail(ds)
str(ds)
summary(ds)
######################################
#3. Decision Tree 생성
######################################
#rpart 메서드 사용 도움말
str(rpart)

#
model <- rpart(formula=form, data=ds[train, vars])

print(model)
summary(model)
plotcp(model)

# Plot Type 1
plot(model, uniform=TRUE)
text(model)

# Plot Type 2
plot(model, uniform=TRUE)
text(model, use.n=TRUE, all=TRUE, cex=.8)

# Plot Type 3
fancyRpartPlot(model)

# Plot Type 4
prp(model)

# Plot Type 5
prp(model, type=2, extra=104, nn=TRUE, fallen.leaves=TRUE,
faclen=0, varlen=0, shadow.col=”grey”, branch.lty=3)

# Plot Type 6
col <- c(“#FD8D3C”, “#FD8D3C”, “#FD8D3C”, “#BCBDDC”,
“#FDD0A2”, “#FD8D3C”, “#BCBDDC”)
prp(model, type=2, extra=104, nn=TRUE, fallen.leaves=TRUE,
faclen=0, varlen=0, shadow.col=”grey”, branch.lty=3, box.col=col)

# Plot Type 7
prp(model, type=1)
prp(model, type=2)
prp(model, type=3)
prp(model, type=4)

#
install.packages(“partykit”, repos=”http://R-Forge.R-project.org“)
library(partykit)
plot(as.party(model))

R – Nonlinear Regression

######################################################
## 19. Nonlinear Regression
## 비선형 회귀 분석에 사용한다.
## 말그대로 데이터가 선형관계가 아닌 비선형 관계일 경우 사용한다.
######################################################
## (1) One-Compartment Model
## 하나의 exp 만 가지고 증가하는 모델
## 데이터 로딩
setwd(“D:/DEVSource/CSV_DATA/”)
data = read.csv(“one_comp.csv”)
print(data)

## 비선형 분석
## nls 함수 사용
## C0, K 등 초기값 사용
## 결과 = a X
## exp : e^x (e = 20172882)
out = nls(conc~C0*exp(-K*time), start=list(C0=41.3, K=0.64), data=data)
summary(out)

## Grapth 로 표현
install.packages(“NRAIA”)
library(NRAIA)
plotfit(out)

## (2) Tow-Compartment Model
## 두개의 exp 를 가지고 더 급하게 증가하는 모델
data2 = read.csv(“two_comp.csv”)
print(data2)

out2 = nls(conc~SSbiexp(time, A1, lrcl, A2, lrc2), data=data2)
summary(out2)

plotfit(out2)
## (3) One-Compartment Model – 급속증가 – 증가둔화 – 감소
data3 = read.csv(“oral_dose.csv”)
print(data3)
out3 = nls(conc~SSfol(Dose=4.4, time, lKe, lKa, LCl), data=data3)
summary(out3)

plotfit(out3)
## (4)Michaelis-Menten model
data4 = read.csv(“MM.csv”)
out4 = nls(rate~SSmicmen(conc, Vm, K), data=data4)
summary(out4)
plotfit(out4)

## 초기 변수가 필요없는 self-starting function
1. SSasymp()
Description
This selfStart model evaluates the asymptotic regression function and its gradient. It has an initial attribute that will evaluate initial estimates of the parameters Asym, R0, and lrc for a given set of data.

Usage
SSasymp(input, Asym, R0, lrc)

Arguments
input a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the horizontal asymptote on the right side (very large values of input).

R0 a numeric parameter representing the response when input is zero.

lrc a numeric parameter representing the natural logarithm of the rate constant.

Value
a numeric vector of the same length as input. It is the value of the expression Asym+(R0-Asym)*exp(-exp(lrc)*input). If all of the arguments Asym, R0, and lrc are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples

Lob.329 <- Loblolly[ Loblolly$Seed == “329”, ]
SSasymp( Lob.329$age, 100, -8.5, -3.2 ) # response only
Asym <- 100
resp0 <- -8.5
lrc <- -3.2
SSasymp( Lob.329$age, Asym, resp0, lrc ) # response and gradient
getInitial(height ~ SSasymp( age, Asym, resp0, lrc), data = Lob.329)
## Initial values are in fact the converged values
fm1 <- nls(height ~ SSasymp( age, Asym, resp0, lrc), data = Lob.329)
summary(fm1)
2. SSasympOff()
Description
This selfStart model evaluates an alternative parametrization of the asymptotic regression function and the gradient with respect to those parameters. It has an initial attribute that creates initial estimates of the parameters Asym, lrc, and c0.

Usage
SSasympOff(input, Asym, lrc, c0)

Arguments
input a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the horizontal asymptote on the right side (very large values of input).

lrc a numeric parameter representing the natural logarithm of the rate constant.

c0 a numeric parameter representing the input for which the response is zero.

Value
a numeric vector of the same length as input. It is the value of the expression Asym*(1 – exp(-exp(lrc)*(input – c0))). If all of the arguments Asym, lrc, and c0 are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart; example(SSasympOff) gives graph showing the SSasympOff parametrization, where phi_1 is Asymp, phi_3 is c0.

Examples
CO2.Qn1 <- CO2[CO2$Plant == “Qn1”, ]
SSasympOff(CO2.Qn1$conc, 32, -4, 43) # response only
Asym <- 32; lrc <- -4; c0 <- 43
SSasympOff(CO2.Qn1$conc, Asym, lrc, c0) # response and gradient
getInitial(uptake ~ SSasympOff(conc, Asym, lrc, c0), data = CO2.Qn1)
## Initial values are in fact the converged values
fm1 <- nls(uptake ~ SSasympOff(conc, Asym, lrc, c0), data = CO2.Qn1)
summary(fm1)

3. SSasympOrig()
Description
This selfStart model evaluates the asymptotic regression function through the origin and its gradient. It has an initial attribute that will evaluate initial estimates of the parameters Asym and lrc for a given set of data.

Usage
SSasympOrig(input, Asym, lrc)

Arguments
input a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the horizontal asymptote.

lrc a numeric parameter representing the natural logarithm of the rate constant.

Value
a numeric vector of the same length as input. It is the value of the expression Asym*(1 – exp(-exp(lrc)*input)). If all of the arguments Asym and lrc are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
Lob.329 <- Loblolly[ Loblolly$Seed == “329”, ]
SSasympOrig(Lob.329$age, 100, -3.2) # response only
Asym <- 100; lrc <- -3.2
SSasympOrig(Lob.329$age, Asym, lrc) # response and gradient
getInitial(height ~ SSasympOrig(age, Asym, lrc), data = Lob.329)
## Initial values are in fact the converged values
fm1 <- nls(height ~ SSasympOrig(age, Asym, lrc), data = Lob.329)
summary(fm1)

4. SSbiexp()

Description
This selfStart model evaluates the biexponential model function and its gradient. It has an initial attribute that creates initial estimates of the parameters A1, lrc1, A2, and lrc2.

Usage
SSbiexp(input, A1, lrc1, A2, lrc2)

Arguments
input a numeric vector of values at which to evaluate the model.

A1 a numeric parameter representing the multiplier of the first exponential.

lrc1 a numeric parameter representing the natural logarithm of the rate constant of the first exponential.

A2 a numeric parameter representing the multiplier of the second exponential.

lrc2 a numeric parameter representing the natural logarithm of the rate constant of the second exponential.

Value
a numeric vector of the same length as input. It is the value of the expression A1*exp(-exp(lrc1)*input)+A2*exp(-exp(lrc2)*input). If all of the arguments A1, lrc1, A2, and lrc2 are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
Indo.1 <- Indometh[Indometh$Subject == 1, ]
SSbiexp( Indo.1$time, 3, 1, 0.6, -1.3 ) # response only
A1 <- 3; lrc1 <- 1; A2 <- 0.6; lrc2 <- -1.3
SSbiexp( Indo.1$time, A1, lrc1, A2, lrc2 ) # response and gradient
print(getInitial(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = Indo.1),
digits = 5)
## Initial values are in fact the converged values
fm1 <- nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = Indo.1)
summary(fm1)

5. SSfol()

Description
This selfStart model evaluates the first-order compartment function and its gradient. It has an initial attribute that creates initial estimates of the parameters lKe, lKa, and lCl.

Usage
SSfol(Dose, input, lKe, lKa, lCl)

Arguments
Dose a numeric value representing the initial dose.

input a numeric vector at which to evaluate the model.

lKe a numeric parameter representing the natural logarithm of the elimination rate constant.

lKa a numeric parameter representing the natural logarithm of the absorption rate constant.

lCl a numeric parameter representing the natural logarithm of the clearance.

Value
a numeric vector of the same length as input, which is the value of the expression

Dose * exp(lKe+lKa-lCl) * (exp(-exp(lKe)*input) – exp(-exp(lKa)*input))
/ (exp(lKa) – exp(lKe))

If all of the arguments lKe, lKa, and lCl are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.
Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
Theoph.1 <- Theoph[ Theoph$Subject == 1, ]
SSfol(Theoph.1$Dose, Theoph.1$Time, -2.5, 0.5, -3) # response only
lKe <- -2.5; lKa <- 0.5; lCl <- -3
SSfol(Theoph.1$Dose, Theoph.1$Time, lKe, lKa, lCl) # response and gradient
getInitial(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph.1)
## Initial values are in fact the converged values
fm1 <- nls(conc ~ SSfol(Dose, Time, lKe, lKa, lCl), data = Theoph.1)
summary(fm1)

6. SSfpl()

Description
This selfStart model evaluates the four-parameter logistic function and its gradient. It has an initial attribute that will evaluate initial estimates of the parameters A, B, xmid, and scal for a given set of data.

Usage
SSfpl(input, A, B, xmid, scal)

Arguments
input a numeric vector of values at which to evaluate the model.

A a numeric parameter representing the horizontal asymptote on the left side (very small values of input).

B a numeric parameter representing the horizontal asymptote on the right side (very large values of input).

xmid a numeric parameter representing the input value at the inflection point of the curve. The value of SSfpl will be midway between A and B at xmid.

scal a numeric scale parameter on the input axis.

Value
a numeric vector of the same length as input. It is the value of the expression A+(B-A)/(1+exp((xmid-input)/scal)). If all of the arguments A, B, xmid, and scal are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
Chick.1 <- ChickWeight[ChickWeight$Chick == 1, ]
SSfpl(Chick.1$Time, 13, 368, 14, 6) # response only
A <- 13; B <- 368; xmid <- 14; scal <- 6
SSfpl(Chick.1$Time, A, B, xmid, scal) # response and gradient
print(getInitial(weight ~ SSfpl(Time, A, B, xmid, scal), data = Chick.1),
digits = 5)
## Initial values are in fact the converged values
fm1 <- nls(weight ~ SSfpl(Time, A, B, xmid, scal), data = Chick.1)
summary(fm1)
7. SSgompertz()

Description
This selfStart model evaluates the Gompertz growth model and its gradient. It has an initial attribute that creates initial estimates of the parameters Asym, b2, and b3.

Usage
SSgompertz(x, Asym, b2, b3)

Arguments
x a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the asymptote.

b2 a numeric parameter related to the value of the function at x = 0

b3 a numeric parameter related to the scale the x axis.

Value
a numeric vector of the same length as input. It is the value of the expression Asym*exp(-b2*b3^x). If all of the arguments Asym, b2, and b3 are names of objects the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
Douglas Bates

See Also
nls, selfStart

Examples
DNase.1 <- subset(DNase, Run == 1)
SSgompertz(log(DNase.1$conc), 4.5, 2.3, 0.7) # response only
Asym <- 4.5; b2 <- 2.3; b3 <- 0.7
SSgompertz(log(DNase.1$conc), Asym, b2, b3) # response and gradient
print(getInitial(density ~ SSgompertz(log(conc), Asym, b2, b3),
data = DNase.1), digits = 5)
## Initial values are in fact the converged values
fm1 <- nls(density ~ SSgompertz(log(conc), Asym, b2, b3),
data = DNase.1)
summary(fm1)
8. SSlogis()

Description
This selfStart model evaluates the logistic function and its gradient. It has an initial attribute that creates initial estimates of the parameters Asym, xmid, and scal.

Usage
SSlogis(input, Asym, xmid, scal)

Arguments
input a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the asymptote.

xmid a numeric parameter representing the x value at the inflection point of the curve. The value of SSlogis will be Asym/2 at xmid.

scal a numeric scale parameter on the input axis.

Value
a numeric vector of the same length as input. It is the value of the expression Asym/(1+exp((xmid-input)/scal)). If all of the arguments Asym, xmid, and scal are names of objects the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
Chick.1 <- ChickWeight[ChickWeight$Chick == 1, ]
SSlogis(Chick.1$Time, 368, 14, 6) # response only
Asym <- 368; xmid <- 14; scal <- 6
SSlogis(Chick.1$Time, Asym, xmid, scal) # response and gradient
getInitial(weight ~ SSlogis(Time, Asym, xmid, scal), data = Chick.1)
## Initial values are in fact the converged values
fm1 <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = Chick.1)
summary(fm1)
9. SSmicmen()

Description
This selfStart model evaluates the Michaelis-Menten model and its gradient. It has an initial attribute that will evaluate initial estimates of the parameters Vm and K

Usage
SSmicmen(input, Vm, K)

Arguments
input a numeric vector of values at which to evaluate the model.

Vm a numeric parameter representing the maximum value of the response.

K a numeric parameter representing the input value at which half the maximum response is attained. In the field of enzyme kinetics this is called the Michaelis parameter.

Value
a numeric vector of the same length as input. It is the value of the expression Vm*input/(K+input). If both the arguments Vm and K are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
José Pinheiro and Douglas Bates

See Also
nls, selfStart

Examples
PurTrt <- Puromycin[ Puromycin$state == “treated”, ]
SSmicmen(PurTrt$conc, 200, 0.05) # response only
Vm <- 200; K <- 0.05
SSmicmen(PurTrt$conc, Vm, K) # response and gradient
print(getInitial(rate ~ SSmicmen(conc, Vm, K), data = PurTrt), digits = 3)
## Initial values are in fact the converged values
fm1 <- nls(rate ~ SSmicmen(conc, Vm, K), data = PurTrt)
summary(fm1)
## Alternative call using the subset argument
fm2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = Puromycin,
subset = state == “treated”)
summary(fm2)

10. SSweibull()

Description
This selfStart model evaluates the Weibull model for growth curve data and its gradient. It has an initial attribute that will evaluate initial estimates of the parameters Asym, Drop, lrc, and pwr for a given set of data.

Usage
SSweibull(x, Asym, Drop, lrc, pwr)

Arguments
x a numeric vector of values at which to evaluate the model.

Asym a numeric parameter representing the horizontal asymptote on the right side (very small values of x).

Drop a numeric parameter representing the change from Asym to the y intercept.

lrc a numeric parameter representing the natural logarithm of the rate constant.

pwr a numeric parameter representing the power to which x is raised.

Details
This model is a generalization of the SSasymp model in that it reduces to SSasymp when pwr is unity.

Value
a numeric vector of the same length as x. It is the value of the expression Asym-Drop*exp(-exp(lrc)*x^pwr). If all of the arguments Asym, Drop, lrc, and pwr are names of objects, the gradient matrix with respect to these names is attached as an attribute named gradient.

Author(s)
Douglas Bates

References
Ratkowsky, David A. (1983), Nonlinear Regression Modeling, Dekker. (section 4.4.5)

See Also
nls, selfStart, SSasymp

Examples
Chick.6 <- subset(ChickWeight, (Chick == 6) & (Time > 0))
SSweibull(Chick.6$Time, 160, 115, -5.5, 2.5) # response only
Asym <- 160; Drop <- 115; lrc <- -5.5; pwr <- 2.5
SSweibull(Chick.6$Time, Asym, Drop, lrc, pwr) # response and gradient
getInitial(weight ~ SSweibull(Time, Asym, Drop, lrc, pwr), data = Chick.6)
## Initial values are in fact the converged values
fm1 <- nls(weight ~ SSweibull(Time, Asym, Drop, lrc, pwr), data = Chick.6)
summary(fm1)

R – Logistic Regression

######################################################
## 13. Logistic Regression
## 종속변수가 0,1 / 죽다,살다 / 있다,없다 와 같은 Binary 인 경우 사용
## odds 의 log- 변환을 종속변수로 모형화 한다.
######################################################

## CASE1 : Respire.csv
## (1) 테스테 데이터 추출
setwd(“D:/DEVSource/CSV_DATA/”)
data1 = read.csv(“respire1.csv”)
data2 = read.csv(“respire2.csv”)
print(data1)
print(data2)

## (2) logistic test
## 귀무가설 : treat 종류에 따른 결과가 같다
## P-value : 6.7e-06 결과는 매우 유의한 차이가 있다
out = glm(outcome ~ treat, family=binomial, data=data)
summary(out)

## (3) 얼마나 더 좋은가?
## 새로운 약은 6배더 효과가 있다
exp(coef(out))

## (4) 신뢰 구간은?
## 2.8배 ~ 13.41배 효과가 있을 확률은 95% 이다
exp(confint(out))

## CASE2 : toxic.csv
data3 = read.csv(“toxic.csv”)
print(data3)

with(data3,tapply(count*response, dose, sum))
with(data3,tapply(count, dose, sum))
with(data3, tapply(count*response, dose, sum)/ tapply(count, dose, sum))
## P-value 가 0.331 로 dose 의 변화는 죽을 확률 변화와 연관이 있다
## 그 편차는(?) 1.1051 로 dose 1증가시 죽을 확률은 3.01 배 증가
out = glm(response ~ dose , weights=count, family=binomial, data=data3)
summary(out)
## 신뢰구간
## dose 가 1 증가할 때 죽을 확률은
## 1.16 ~ 9.33 배 증가할 확률이 95% 이다
exp(confint(out, parm=”dose”))

R – ANCOVA (analysis of covariance)공분산 분석

######################################################
## 13. ANCOVA (analysis of covariance)공분산 분석
## 사회과학과 같이 다른 변인의 통제가 쉽지 않은 상황에서
## 연속형 변수를 추가하여 오차를 줄이고 검정력을 높이는 방법
######################################################

## (1) 테스트 데이터 추출
setwd(“D:/DEVSource/CSV_DATA/”)
data = read.csv(“anorexia.csv”)
print(data)

## (2) 데이터에 대한 일원분산 분석
## 3개의 데이터가 같은 평균을 갖지는 않음
boxplot(Postwt – Prewt ~ Treat, data=data)
out = lm(Postwt – Prewt ~ Treat, data=data)
anova(out)
summary(out)

##전체 차이 조사 결과
##FT 와 Cont 사이에 유의한 차이가 존재함
library(multcomp)
tukey = glht(out, linfct=mcp(Treat=”Tukey”))
summary(tukey)

## Treat Level 조정
levels(data$Treat)
data$Treat = relevel(data$Treat, ref=”Cont”)
levels(data$Treat)

## (3) 공변량을 추가하여 분석
## 공변량을(현재 몸무게) 추가하여 분석을 하여도 치료 방법 별로 유의한 차이가 있다는 결론
## Treat 의 P-Value = 0.0008438
## Prewt 의 P-Value = 0.0019364
## 단 이전 몸무게 또한 유의한 차이가 있다는 결론이 있음
out2 = lm(Postwt ~ Prewt + Treat, data=data)
anova(out2)
summary(out2)

dunnett = glht(out2, linfct=mcp(Treat=”Dunnett”))
summary(dunnett )
plot(dunnett )

R – Two-way ANOVA

######################################################
## 12. Two-way ANOVA
## 두개의 그룹변수를 갖는 데이터의 분석에 사용
######################################################

## (1) 테스트 데이터 추출
setwd(“D:/DEVSource/CSV_DATA/”)
data = read.csv(“warpbreaks.csv”)
print(data)

## (2) wool 과 tension 두개의 그룹으로 구성된 데이터이다.
## wool 과 tension 의 level 지정 상태를 확인한다
levels(data$tension)
levels(data$wool)

## (3) tension 의 level 을 L, M , H 순서로 변경한다.
data$tension = factor(data$tension, levels=c(“L”,”M”,”H”))
levels(data$tension)

## (4) 각 그룹의 평균을 구한다
attach(data)
tapply(breaks, wool, mean)
tapply(breaks, tension, mean)
tapply(breaks, list(wool,tension), mean)

boxplot(breaks ~ wool+tension, col=”red”, data=data)

## 그룹 Interaction
interaction.plot(tension, wool, breaks, col=c(“red”, “blue”))

## (5) wool 과 tension 사이의 교호작용이 있는 경우
## wool*tension 을 넣어 줘야 함
##Response: breaks
## Df Sum Sq Mean Sq F value Pr(>F)
##wool 1 450.7 450.67 3.7653 0.0582130 .
##tension 2 2034.3 1017.13 8.4980 0.0006926 ***
##wool:tension 2 1002.8 501.39 4.1891 0.0210442 *
##Residuals 48 5745.1 119.69
##—
##Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘
## 분석 결과를 보면
## wool 은 그룹이 바뀌어도 유의한 차이가 없다.
## 일원 분석에서 보면 그룹간의 평균이 같다라는 식의 결론.
## 하지만 여기서는 유의한 차이가 없음이 되어버린다.
## tension 의 p-value 는 0.0006926 으로 그룹간의 평균이 전혀 달라
## 유의한 차이가 있다가 된다.
## Wool 과 tension 사이의 교호관계를 나타내는 wool:tension 은
## p-value 0.210442 로 유의한 차이가 있음으로써 교호관계가 있다는 가설 또한
## 설명이 된다.
out2 = lm(breaks ~ wool + tension + wool*tension)
anova(out2)

## 교호 작용이 없는 경우라면 아래와 같이 wool + tension 으로 관계식을 만들면 되겠다.
summary(lm(breaks ~ wool + tension ))

## (6) 회귀 진단
## 판단식에 대한 회귀 진단
par(mfrow=c(2,2))
plot(out2)
## (7) 정규성 판단
## p-value : 0.8162
## 회귀식의 잔차는 정규성을 만족함
## 식의 추정값이 실제 표본 값과 차이가 일정, 균등 함
shapiro.test(resid(out2))

R – Kruskal-Wallis Test

######################################################
## 11. Kruskal-Wallis Test
## 독립변수가 2개 이상인 경우에 각각의 그룹의 평균이 같은가를 결정함
## 일원 테스트에서 F-TEST 는 그룹간 변동량과 그룹내 변동량의 비
## – (각그룹평균 – 전체평균)^2 + (각그룹평균 – 전체평균)^2 과
## – (각각 변수 – 각 그룹 평균)^2
## 프로그램으로는 anova(lm(값 ~ 그룹)) 하면 끝나버린다.
## 여기서 귀무가설을 모든 그룹의 평균이 같다는 것이고
## 대립가설을 차이가 있는 어떤 그룹이 있다이다
## ※ 일원 분산 분석과 거의 동일한 방법으로 수행하지만
## 주로 데이터가 정규성을 띄지 않을때 Rank 방법을 이용 분석시 이용
######################################################

## (1) 테스테 데이터 추출
setwd(“D:/DEVSource/CSV_DATA/”)
data = read.csv(“PlantGrowth.csv”)
print(data)

## (2) 데이터 형태 파악
## 그룹별 평균 확인
## 그룹별 분포도 box plot 확인
## 그룹별 shapiroo test 를 통한 정규성 검사
## 3개의 그룹은 모두 정규성을 만족하는 데이터
with(data, tapply(weight, group, mean))
boxplot(weight ~ group, col=”red”, data=data)
shapiro.test(data[data$group==”ctrl” , 2])
shapiro.test(data[data$group==”trt1″ , 2])
shapiro.test(data[data$group==”trt2″ , 2])
## (3) Kruskal-Waillis Test
## 모든 데이터에 순위를 부여하고
## 그룹별 순위의 분포를 분석함
## P-Value : 0.1842
## 귀무가설 : 모든 그룹의 평균은 같다를 기각함
kruskal.test(weight ~ group , data=data)

R – Multiple Rgression

######################################################
## 8. Multiple Regression
## 독립변수가 2개 이상인 경우에 사용하며 아래와 같은 방법이 있음
## – 이론을 미리 정해 놓고 분석하는 Conframtory 분석
## – 이론을 정하지 않고 가장 적절한 모형을 고르는 Exploratory 분석
## – 다중회귀 모형의 특수한 경우인 다항회귀분석
## CASE 분석은 아래와 같은 원칙을 따른다
## – 종속변수와 독립변수간의 상관관계가 높아야 한다.
## – 독립변수들 간의 상관관계는 낮아야 한다.
## – 소수 정예의 원칙을 따른다
## 모형 도출 방법은 아래의 방법들을 사용한다.
## – all possible regression : 모든 방안을 만들어 놓고 비교하는 방법
## – forward selection : 가장 유의한 변수 부터 덜 유의한 변수 순으로 하나씩 추가
## – backward selection : 모든 변수를 넣어 놓고 가장 기여도가 낮은 변수를 하나씩 제거
## – stepwise selection : forward 와 backward 조합
######################################################

## (1) Data Set Attitude 데이터 사용
setwd(“D:/DEVSource/CSV_DATA/”)
data = read.csv(“attitude.csv”)
print(data)
## (2) 전체 변수의 기여도를 확인
## anova 를 사용
## P-Value 가 높을 수로 기여도가 낮음
out1 = lm(rating ~ . , data=data)
anova(out1)

## (3) 일부 기여도가 낮은 데이터를 제외하고 선형식 구성
out2 = lm(rating ~ complaints + learning + advance , data=data)

## (4) 일부 기여도가 낮은 데이터를 제외하고 선형식 구성
out3 = lm(rating ~ complaints + learning , data=data)

## (5) 일부 기여도가 낮은 데이터를 제외하고 선형식 구성
out4 = lm(rating ~ learning , data=data)

## (6) 모델간 비교를 통해 제거한 변수의 기여도를 확인
## – rasies , critical, privaliges 는 기여도가 거의 없음
## – advance 도 기여도가 높지 않음
## – complaints 는 기여도가 매우 높음
anova(out4, out3, out2, out1)

## (7) Summary 로 모델을 확인
##모든 항목을 사용할 경우 조정R : 0.65
##out2 : 0.693
##out3 : 0.684
##out4 : 0.0002311
##out2 의 항목이 가장 높은 조정 R 을 갖는다.
summary(out1)
summary(out2)
summary(out3)
summary(out4)

## (8) 다중선형회귀식 확인
## out2 모델을 사용하여 예측을 진행해 본다.
## 모델 분석에 사용한 데이터를 그대로 넣는 경우
print(data[1,c(3,5,8)])
plot(predict(out2, newdata=data[1,c(3,5,8)], interval=”prediction”)[,1])

## complaints , learning, advance 를 임의로 정하여 예측
## 실행 결과 예
## 100, 100, 100 => 88.39
## 10, 10, 10 => 21.0594
## 10, 100, 100 => 32.3486
input = data.frame(10,100,100)
colnames(input) <- c(“complaints”, “learning”, “advance”)
print(predict(out2, newdata=input , interval=”prediction”))

R – Simple Linear Regression

######################################################
## 7. Simple LinearRegression
## 연속 데이터에 대한 분석 방법
## 종속 변수와 독립변수의 관계가 1:1 이다.
## 예를 들면 독립변수는 차량의 속도 종속변수는 차량의 제동 거리가 될 수 있겠다.
######################################################
## (1) speed ,dist 두개의 컬럼으로 구성된 데이터 셋을 만들었다
speed = c(10,20,30,40,50,60,70)
dist = c(13,18,25,38,47,61,105)
data = data.frame(speed ,dist)
print(data)

## (2) 선형회귀식을 만든다
## – speed 1.371 0.237 5.786 0.00217 **
## ==> 기울기는 1.371 이며 이 기울기에 대한 P-Value 는 0.00217 로 회기식은 유의함
## – Adjusted R-squared: 0.8441
## ==> 조정된 설명력 : 해당 모형의 설득력을 지수화
## 아래는 회귀식의 Summary(out) 결과이다.
##
##Residuals:
## 1 2 3 4 5 6 7
## 10.286 1.571 -5.143 -5.857 -10.571 -10.286 20.000
##
##Coefficients:
## Estimate Std. Error t value Pr(>|t|)
##(Intercept) -11.000 10.600 -1.038 0.34696
##speed 1.371 0.237 5.786 0.00217 **
##—
##Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
##Residual standard error: 12.54 on 5 degrees of freedom
##Multiple R-squared: 0.87, Adjusted R-squared: 0.8441
##F-statistic: 33.48 on 1 and 5 DF, p-value: 0.002171

out = lm(dist ~ speed , data=data)
out
summary(out)
## (3) 회귀식 확인
## plot(dist ~ speed, data=data, col=”blue”) 은 파란점으로 실제 데이터를 표시
## abline(out, col=”red”) 는 선형회귀식을 붉은선으로 표시한다.
## 앞에서 설명한 것처럼 오차가 가장적은 방향으로 식이 성립된다.

plot(dist ~ speed, data=data, col=”blue”)
abline(out, col=”red”)

## (4) 회귀식 보정
## 이 케이스에서 실제로 속도가 0 이라면 제동거리도 0이 되는 것이 맞다
## -1 을 통해서 Y 절편을 제거 할 수 있다.
summary(lm(dist ~ speed -1, data=data))
out = lm(dist ~ speed – 1 , data=data)
plot(dist ~ speed -1, data=data)
abline(out, col=”red”)

## (4) sqrt 변환
## sqrt 를 적용시 데이터의 퍼지는 정도가 줄어드는 것을 볼 수 있다.
## 단, 본 케이스의 경우 log 적용시 정규분포도가 오히려 나빠짐
plot(dist ~ speed -1, data=data)
plot(log(dist) ~ speed – 1, data=data)
plot(sqrt(dist) ~ speed- 1, data=data)

##미 적용
out = lm(dist ~ speed – 1 , data=data)
plot(dist ~ speed – 1, data=data)
abline(out, col=”red”)

##log 적용
out = lm(log(dist) ~ speed – 1 , data=data)
plot(log(dist) ~ speed – 1, data=data)
abline(out, col=”red”)

##sqrt 적용
out = lm(sqrt(dist) ~ speed – 1 , data=data)
plot(sqrt(dist) ~ speed – 1, data=data)
abline(out, col=”red”)

## (5) 모델 검증
## out 의 Normal Q-Q 는 직선의 형태를 갖으며 (normal Q-q 표준화된 그래프, 정렬)
## 나머지 관계도는 별다른 추세를 보이지 않는다.
## 또한 resid(out)을 통해 구한 잔차는
## 데이터 변환 없을 경우 P-Value 0.0.004615
## log 적용시 P-Value 0.8196
## sqrt 적용시 P-Value 0.9638

par(mfrow=c(2,2))
plot(out)
##qq Normal
plot(resid(out))
qqnorm(resid(out))
qqline(resid(out))
##정규성 검정
shapiro.test(resid(out))

## (6) 결과물 out 의 내용
## names(out) 을 통해서 얻을 수 있는 값을 확인 가능
names(out)
out$rank
out$residual
##예측치
fitted(out)
##잔차
resid(out)

##Residual Standard Error
sqrt(sum(resid(out)^2)/(length(data$dist)-2))

## (7) simulation
## – 편차를 제거하고 sqrt를 적용한 회귀식을 생성
## – rnorm 을 통해 data 수와 같은 수를 갖고, 평균이 거리 예측치이며, 분산은 오차와인
## 임의수 테스트 데이터를 만들어 낸다.
## – 다시 만들어낸 dist set 과 speed 의 관계를 그래프로 표현

## 회귀식 생성
out = lm(sqrt(dist) ~ speed-1, data=data)
## 분산
print(summary(out)$sigma)

## 평균
print(mean(fitted(out)))

##위와 같은 분산, 평균을 갖는 새로운 dist 데이터 생성
dist = rnorm(n=nrow(data), mean=fitted(out), sd=summary(out)$sigma)^2
plot(dist ~ data$speed)

## 모델과 R 명령어
##Model : Y = a + bX
##R Code : Y ~ X
##
##Model : Y = bX
##R Code : Y ~ X – 1
##
##Model : sqrt(Y) = a + bX
##R Code : sqrt(Y) ~ X
##
##Model : sqrt(Y) = bX
##R Code : sqrt(Y) ~ X – 1
##
##Model : lnY = a + bX
##R Code : log(Y) ~ X

R – Correlation Analysis

######################################################
## 6. Correlation Analysis
## 상관분석이라 함은 두 변수가 같은 방향으로 움직이는지 다른 방향으로 움직이는지를 보여줌
## 보통 상관계수는 데이터 측량 단위에 영향을 받기 때문에 이를 -1 ~ 1 사이로 보정하는
## Pearson 상관계수를 많이 사용한다.
## cf) Pearson 상관계수 외에도 비모수 방법인 Kendall , Spearman 도 상관계수도 존재함
######################################################

## (1)기본 내장 데이터를 테스트 데이터로 사용한다.
## 아래는 테이블의 데어터 내용이다.
## rating complaints privileges learning raises critical advance
##1 43 51 30 39 61 92 45
##2 63 64 51 54 63 73 47
##3 71 70 68 69 76 86 48
##4 61 63 45 47 54 84 35
##5 81 78 56 66 71 83 47
##6 43 55 49 44 54 49 34
##7 58 67 42 56 66 68 35
##8 71 75 50 55 70 66 41
##9 72 82 72 67 71 83 31
##10 67 61 45 47 62 80 41
##11 64 53 53 58 58 67 34
##12 67 60 47 39 59 74 41
print(attitude)
## (2)모든 컬럼간의 공분산을 구한다.
## rating complaints privileges learning raises critical
##rating 148.17126 133.77931 63.46437 89.10460 74.68851 18.84253
##complaints 133.77931 177.28276 90.95172 93.25517 92.64138 24.73103
##privileges 63.46437 90.95172 149.70575 70.84598 56.67126 17.82529
##learning 89.10460 93.25517 70.84598 137.75747 78.13908 13.46782
##raises 74.68851 92.64138 56.67126 78.13908 108.10230 38.77356
##critical 18.84253 24.73103 17.82529 13.46782 38.77356 97.90920
##advance 19.42299 30.76552 43.21609 64.19770 61.42299 28.84598

cov(attitude)

## (3)모든 컬럼간의 상관계수를 구한다 ( -1 ~ 1)
## 기본은 pearson 상관계수이다.
## rating complaints privileges learning raises critical
##rating 1.0000000 0.8254176 0.4261169 0.6236782 0.5901390 0.1564392
##complaints 0.8254176 1.0000000 0.5582882 0.5967358 0.6691975 0.1877143
##privileges 0.4261169 0.5582882 1.0000000 0.4933310 0.4454779 0.1472331
##learning 0.6236782 0.5967358 0.4933310 1.0000000 0.6403144 0.1159652
##raises 0.5901390 0.6691975 0.4454779 0.6403144 1.0000000 0.3768830
##critical 0.1564392 0.1877143 0.1472331 0.1159652 0.3768830 1.0000000
##advance 0.1550863 0.2245796 0.3432934 0.5316198 0.5741862 0.2833432

cor(attitude)

## (4)P-Value
## rating 과 complaints 사이의 P-Value 를 구해보겠다.
## pearson, kendall, spearman 모두 매우 낮은 P-Value 를 갖음
## 귀무가설 : 둘의 상관계수는 0 이다.
## 대립가설 : 둘은 상관 관계가 있다.
## 결론 : 귀무가설을 기각하고 둘의 상관관계가 있다는 대립가설을 채택함
with(attitude, cor.test(rating , complaints , method=”pearson”, use=”pairwise.complete.obs”))

with(attitude, cor.test(rating , complaints , method=”kendall”, use=”pairwise.complete.obs”))

with(attitude, cor.test(rating , complaints , method=”spearman”, use=”pairwise.complete.obs”))