이번 포스팅에서는 R을 사용하여 예측이나 분류 모델링을 할 때 기본적으로 필요한 두가지 작업인 


(1) DataFrame을 Train set, Test set 으로 분할하기 (Split a DataFrame into Train and Test set)

   - (1-1) 무작위 샘플링 (random sampling)

   - (1-2) 순차 샘플링 (sampling in order)

   - (1-3) 층화 무작위 샘플링 (stratified random sampling)


(2) 여러개의 숫자형 변수를 가진 DataFrame을 변환하기 (Transformation of DataFrame)

  - (2-1) z-변환 (z-transformation, standardization)

  - (2-2) [0-1] 변환 ([0-1] transformation, normalization)


에 대해서 소개하겠습니다. 



예제로 사용할 Cars93 DataFrame을 MASS 패키지로 부터 불러오겠습니다. 변수가 무척 많으므로 예제를 간단하게 하기 위해 설명변수 X로 'Price', 'Horsepower', 'RPM', 'Length', 'Type', 'Origin' 만을 subset 하여 가져오고, 반응변수 y 로는 'MPG.highway' 변수를 사용하겠습니다. 



# get Cars93 DataFrame from MASS package

library(MASS)

data(Cars93)

str(Cars93)

'data.frame': 93 obs. of 27 variables: $ Manufacturer : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4... $ Model : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1... $ Type : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3... $ Min.Price : num 12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ... $ Price : num 15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ... $ Max.Price : num 18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3... $ MPG.city : int 25 18 20 19 22 22 19 16 19 16 ... $ MPG.highway : int 31 25 26 26 30 31 28 25 27 25 ... $ AirBags : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2... $ DriveTrain : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3... $ Cylinders : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4... $ EngineSize : num 1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ... $ Horsepower : int 140 200 172 172 208 110 170 180 170 200 ... $ RPM : int 6300 5500 5500 5500 5700 5200 4800 4000 4800... $ Rev.per.mile : int 2890 2335 2280 2535 2545 2565 1570 1320 1690... $ Man.trans.avail : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1... $ Fuel.tank.capacity: num 13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ... $ Passengers : int 5 5 5 6 4 6 6 6 5 6 ... $ Length : int 177 195 180 193 186 189 200 216 198 206 ... $ Wheelbase : int 102 115 102 106 109 105 111 116 108 114 ... $ Width : int 68 71 67 70 69 69 74 78 73 73 ... $ Turn.circle : int 37 38 37 37 39 41 42 45 41 43 ... $ Rear.seat.room : num 26.5 30 28 31 27 28 30.5 30.5 26.5 35 ... $ Luggage.room : int 11 15 14 17 13 16 17 21 14 18 ... $ Weight : int 2705 3560 3375 3405 3640 2880 3470 4105 3495... $ Origin : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1... 

$ Make : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...



X <- subset(Cars93, select=c('Price', 'Horsepower', 'RPM', 'Length', 'Type', 'Origin'))

head(X)

A data.frame: 6 × 6
PriceHorsepowerRPMLengthTypeOrigin
<dbl><int><int><int><fct><fct>
15.91406300177Smallnon-USA
33.92005500195Midsizenon-USA
29.11725500180Compactnon-USA
37.71725500193Midsizenon-USA
30.02085700186Midsizenon-USA
15.71105200189MidsizeUSA



table(X$Origin)

USA non-USA 48 45



y <- Cars93$MPG.highway

y

  1. 31
  2.  
  3. 25
  4.  
  5. 26
  6.  
  7. 26
  8.  
  9. 30
  10.  
  11. 31
  12.  
  13. 28
  14.  
  15. 25
  16.  
  17. 27
  18.  
  19. 25
  20.  
  21. 25
  22.  
  23. 36
  24.  
  25. 34
  26.  
  27. 28
  28.  
  29. 29
  30.  
  31. 23
  32.  
  33. 20
  34.  
  35. 26
  36.  
  37. 25
  38.  
  39. 28
  40.  
  41. 28
  42.  
  43. 26
  44.  
  45. 33
  46.  
  47. 29
  48.  
  49. 27
  50.  
  51. 21
  52.  
  53. 27
  54.  
  55. 24
  56.  
  57. 33
  58.  
  59. 28
  60.  
  61. 33
  62.  
  63. 30
  64.  
  65. 27
  66.  
  67. 29
  68.  
  69. 30
  70.  
  71. 20
  72.  
  73. 30
  74.  
  75. 26
  76.  
  77. 50
  78.  
  79. 36
  80.  
  81. 31
  82.  
  83. 46
  84.  
  85. 31
  86.  
  87. 33
  88.  
  89. 29
  90.  
  91. 34
  92.  
  93. 27
  94.  
  95. 22
  96.  
  97. 24
  98.  
  99. 23
  100.  
  101. 26
  102.  
  103. 26
  104.  
  105. 37
  106.  
  107. 36
  108.  
  109. 34
  110.  
  111. 24
  112.  
  113. 25
  114.  
  115. 29
  116.  
  117. 25
  118.  
  119. 26
  120.  
  121. 26
  122.  
  123. 33
  124.  
  125. 24
  126.  
  127. 33
  128.  
  129. 30
  130.  
  131. 23
  132.  
  133. 26
  134.  
  135. 31
  136.  
  137. 31
  138.  
  139. 23
  140.  
  141. 28
  142.  
  143. 30
  144.  
  145. 41
  146.  
  147. 31
  148.  
  149. 28
  150.  
  151. 27
  152.  
  153. 28
  154.  
  155. 26
  156.  
  157. 38
  158.  
  159. 37
  160.  
  161. 30
  162.  
  163. 30
  164.  
  165. 43
  166.  
  167. 37
  168.  
  169. 32
  170.  
  171. 29
  172.  
  173. 22
  174.  
  175. 33
  176.  
  177. 21
  178.  
  179. 30
  180.  
  181. 25
  182.  
  183. 28
  184.  
  185. 28




  (1) DataFrame을 Train set, Test set 으로 분할하기 (Split a DataFrame into Train and Test set)




(1-1) 무작위 샘플링 (random sampling)



# (1) index for splitting data into Train and Test set

set.seed(1004) # for reprodicibility

train_idx <- sample(1:nrow(X), size=0.8*nrow(X), replace=F) # train-set 0.8, test-set 0.2

test_idx <- (-train_idx)


X_train <- X[train_idx,]

y_train <- y[train_idx]

X_test <- X[test_idx,]

y_test <- y[test_idx]


print(paste0('X_train: ', nrow(X_train)))

print(paste0('y_train: ', length(y_train)))

print(paste0('X_test: ', nrow(X_test)))

print(paste0('y_test: ', length(y_test)))

[Out]:

[1] "X_train: 74" [1] "y_train: 74" [1] "X_test: 19" [1] "y_test: 19"





(1-2) 순차 샘플링 (sampling in order)






 




(1-3) 층화 무작위 샘플링 (stratified random sampling)






 




  (2) 여러개의 숫자형 변수를 가진 DataFrame을 변환하기 (Transformation of DataFrame)


(2-1) z-변환 (z-transformation, standardization)






 




(2-2) [0-1] 변환 ([0-1] transformation, normalization)







 




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)


Posted by R Friend R_Friend

댓글을 달아 주세요

지난번 포스팅에서는 시계열 자료의 구성 요인 (time series component factors)으로서 추세 요인, 순환 요인, 계절 요인, 불규칙 요인에 대해서 소개하였습니다. 


지난번 포스팅에서는 가법 모형으로 가상의 시계열 자료를 만들었다면(time series composition), ==> 이번 포스팅에서는 반대로 시계열 분해(time series decomposition)를 통해 시계열 자료를 추세(순환)(Trend), 계절성(Seasonality), 잔차(Residual)로 분해를 해보겠습니다. 


시계열 분해는 직관적으로 이해하기 쉽고 구현도 쉬워서 시계열 자료 분석의 고전적인 방법론이지만 지금까지도 꾸준히 사용되는 방법론입니다. 시계열 분해의 순서와 방법은 대략 아래와 같습니다. 


(1) 시도표 (time series plot)를 보고 시계열의 주기적 반복/계절성이 있는지, 가법 모형(additive model, y = t + s + r)과 승법 모형(multiplicative model, y = t * s * r) 중 무엇이 더 적합할지 판단을 합니다. 


(가법 모형을 가정할 시)

(2) 시계열 자료에서 추세(trend)를 뽑아내기 위해서 중심 이동 평균(centered moving average)을 이용합니다. 


(3) 원 자료에서 추세 분해값을 빼줍니다(detrend). 그러면 계절 요인과 불규칙 요인만 남게 됩니다. 


(4) 다음에 계절 주기 (seasonal period) 로 detrend 이후 남은 값의 합을 나누어주면 계절 평균(average seasonality)을 구할 수 있습니다. (예: 01월 계절 평균 = (2020-01 + 2021-01 + 2022-01 + 2023-01)/4, 02월 계절 평균 = (2020-02 + 2021-02 + 2022-02 + 2023-02)/4). 


(5) 원래의 값에서 추세와 계절성 분해값을 빼주면 불규칙 요인(random, irregular factor)이 남게 됩니다. 



시계열 분해 후에 추세와 계절성을 제외한 잔차(residual, random/irregular factor) 가 특정 패턴 없이 무작위 분포를 띠고 작은 값이면 추세와 계절성으로 모형화가 잘 되는 것이구요, 시계열 자료의 특성을 이해하고 예측하는데 활용할 수 있습니다. 만약 시계열 분해 후의 잔차에 특정 패턴 (가령, 주기적인 파동을 그린다거나, 분산이 점점 커진다거나 등..) 이 존재한다면 잔차에 대해서만 다른 모형을 추가로 적합할 수도 있겠습니다. 



예제로 사용할 시계열 자료로서 '1차 선형 추세 + 4년 주기 순환 + 1년 단위 계절성 + 불규칙 noise' 의 가법 모형 (additive model)으로 시계열 데이터를 만들어보겠습니다. 



import numpy as np

import pandas as pd

import matplotlib.pyplot as plt


dates = pd.date_range('2020-01-01', periods=48, freq='M')

dates

[Out]:

DatetimeIndex(['2020-01-31', '2020-02-29', '2020-03-31', '2020-04-30', '2020-05-31', '2020-06-30', '2020-07-31', '2020-08-31', '2020-09-30', '2020-10-31', '2020-11-30', '2020-12-31', '2021-01-31', '2021-02-28', '2021-03-31', '2021-04-30', '2021-05-31', '2021-06-30', '2021-07-31', '2021-08-31', '2021-09-30', '2021-10-31', '2021-11-30', '2021-12-31', '2022-01-31', '2022-02-28', '2022-03-31', '2022-04-30', '2022-05-31', '2022-06-30', '2022-07-31', '2022-08-31', '2022-09-30', '2022-10-31', '2022-11-30', '2022-12-31', '2023-01-31', '2023-02-28', '2023-03-31', '2023-04-30', '2023-05-31', '2023-06-30', '2023-07-31', '2023-08-31', '2023-09-30', '2023-10-31', '2023-11-30', '2023-12-31'], dtype='datetime64[ns]', freq='M')

# additive model: trend + cycle + seasonality + irregular factor


timestamp = np.arange(len(dates))

trend_factor = timestamp*1.1

cycle_factor = 10*np.sin(np.linspace(0, 3.14*2, 48))

seasonal_factor = 7*np.sin(np.linspace(0, 3.14*8, 48))

np.random.seed(2004)

irregular_factor = 2*np.random.randn(len(dates))


df = pd.DataFrame({'timeseries': trend_factor + cycle_factor + seasonal_factor + irregular_factor, 

                   'trend': trend_factor, 

                   'cycle': cycle_factor, 

                   'trend_cycle': trend_factor + cycle_factor,

                   'seasonal': seasonal_factor, 

                   'irregular': irregular_factor},

                   index=dates)


df

[Out]:

timeseriestrendcycletrend_cycleseasonalirregular
2020-01-312.5961190.00.0000000.0000000.0000002.596119
2020-02-296.7461601.11.3321982.4321983.5656840.748278
2020-03-318.1121002.22.6406474.8406476.136825-2.865371
2020-04-308.2559413.33.9020217.2020216.996279-5.942358
2020-05-3116.8896554.45.0938349.4938345.9043271.491495
2020-06-3016.1823575.56.19483911.6948393.1655361.321981
2020-07-3114.1280876.67.18540913.785409-0.4561870.798865
2020-08-3111.9433137.78.04788615.747886-3.9506710.146099
2020-09-309.7280958.88.76689217.566892-6.343231-1.495567
2020-10-3112.4834899.99.32961219.229612-6.9665330.220411
2020-11-3012.14180811.09.72601320.726013-5.646726-2.937480
2020-12-3115.14333412.19.94902922.049029-2.751930-4.153764
2021-01-3121.77451613.29.99468423.1946840.910435-2.330604
2021-02-2828.43289214.39.86216424.1621644.318862-0.048134
2021-03-3132.35058315.49.55383224.9538326.5226690.874082
2021-04-3030.59655616.59.07518425.5751846.907169-1.885797
2021-05-3132.51052317.68.43475326.0347535.3651181.110653
2021-06-3030.42551918.77.64395526.3439552.3266241.754939
2021-07-3124.30095819.86.71689026.516890-1.360813-0.855119
2021-08-3120.45091720.95.67008226.570082-4.668691-1.450475
2021-09-3018.87088122.04.52219526.522195-6.674375-0.976939
2021-10-3121.32631023.13.29369026.393690-6.8184381.751059
2021-11-3022.90244824.22.00646926.206469-5.0606991.756678
2021-12-3126.62057825.30.68347825.983478-1.8914262.528526
2022-01-3127.62649926.4-0.65169625.7483041.8054040.072791
2022-02-2831.85892327.5-1.97525325.5247474.9986701.335506
2022-03-3135.93046928.6-3.26359825.3364026.7977043.796363
2022-04-3030.17787029.7-4.49376225.2062386.700718-1.729087
2022-05-3130.01616530.8-5.64381625.1561844.7347640.125217
2022-06-3026.59172931.9-6.69325825.2067421.448187-0.063200
2022-07-3121.11848133.0-7.62337925.376621-2.242320-2.015820
2022-08-3116.63603134.1-8.41759925.682401-5.307397-3.738973
2022-09-3017.68261335.2-9.06175926.138241-6.892132-1.563496
2022-10-3121.16329836.3-9.54437526.755625-6.5545090.962182
2022-11-3022.45567237.4-9.85684427.543156-4.388699-0.698786
2022-12-3126.91952938.5-9.99359528.506405-0.998790-0.588086
2023-01-3133.96462339.6-9.95219129.6478092.6697021.647112
2023-02-2837.45977640.7-9.73336930.9666315.5935590.899586
2023-03-3140.79376641.8-9.34103132.4589696.9572571.377540
2023-04-3043.83841542.9-8.78217134.1178296.3804333.340153
2023-05-3141.30178044.0-8.06675135.9332494.0239751.344556
2023-06-3039.21786645.1-7.20752637.8924740.5451470.780245
2023-07-3135.12550246.2-6.21981339.980187-3.085734-1.768952
2023-08-3133.84192647.3-5.12121942.178781-5.855940-2.480916
2023-09-3038.77051148.4-3.93132944.468671-6.9928031.294643
2023-10-3137.37121649.5-2.67135646.828644-6.179230-3.278198
2023-11-3046.58763350.6-1.36376049.236240-3.6421420.993536
2023-12-3146.40332651.7-0.03185351.668147-0.089186-5.175634





  (1) Python을 이용한 시계열 분해 (Time series decomposition using Python)


Python의 statsmodels 라이브러리를 사용해서 가법 모형(additive model) 가정 하에 시계열 분해를 해보겠습니다. 



from statsmodels.tsa.seasonal import seasonal_decompose


ts = df.timeseries

result = seasonal_decompose(ts, model='additive')


plt.rcParams['figure.figsize'] = [12, 8]

result.plot()

plt.show()





원래의 시계열 구성요소(추세+순환, 계절성, 불규칙 요인)와 시계열 분해(time series decomposition)를 통해 분리한 추세(&순환), 계절성, 잔차(불규칙 요인)를 겹쳐서 그려보았습니다. (즉, 원래 데이터의 추세요인과 시계열 분해를 통해 분리한 추세를 겹쳐서 그려보고, 원래 데이터의 계절요인과 시계열 분해를 통해 분리한 계절을 겹쳐서 그려보고, 원래 데이터의 불규칙 요인과 시계열 분해를 통해 분리한 잔차를 겹쳐서 그려봄) 


원래의 데이터와 얼추 비슷하게, 그럴싸하게 시계열 분해를 한 것처럼 보이지요?   



# ground truth & timeseries decompostion all together

# -- observed data

plt.figure(figsize=(12, 12))

plt.subplot(4,1, 1)

result.observed.plot()

plt.grid(True)

plt.ylabel('Observed', fontsize=14)


# -- trend & cycle factor

plt.subplot(4, 1, 2)

result.trend.plot()        # from timeseries decomposition

df.trend_cycle.plot()     # groud truth

plt.grid(True)

plt.ylabel('Trend', fontsize=14)


# -- seasonal factor

plt.subplot(4, 1, 3)

result.seasonal.plot()  # from timeseries decomposition

df.seasonal.plot()        # groud truth

plt.grid(True)

plt.ylabel('Seasonality', fontsize=14)


# -- irregular factor (noise)

plt.subplot(4, 1, 4)

result.resid.plot()    # from timeseries decomposition

df.irregular.plot()    # groud truth

plt.grid(True)

plt.ylabel('Residual', fontsize=14)


plt.show()




원래의 관측치(observed), 추세(trend), 계절성(seasonal), 잔차(residual) 데이터 아래처럼 시계열 분해한 객체에서 obsered, trend, seasonal, resid 라는 attributes 를 통해서 조회할 수 있습니다. 


Observed

 Trend ( & Cycle)

 print(result.observed)

[Out]:
2020-01-31     2.596119
2020-02-29     6.746160
2020-03-31     8.112100
2020-04-30     8.255941
2020-05-31    16.889655
2020-06-30    16.182357
2020-07-31    14.128087
2020-08-31    11.943313
2020-09-30     9.728095
2020-10-31    12.483489
2020-11-30    12.141808
2020-12-31    15.143334
2021-01-31    21.774516
2021-02-28    28.432892
2021-03-31    32.350583
2021-04-30    30.596556
2021-05-31    32.510523
2021-06-30    30.425519
2021-07-31    24.300958
2021-08-31    20.450917
2021-09-30    18.870881
2021-10-31    21.326310
2021-11-30    22.902448
2021-12-31    26.620578
2022-01-31    27.626499
2022-02-28    31.858923
2022-03-31    35.930469
2022-04-30    30.177870
2022-05-31    30.016165
2022-06-30    26.591729
2022-07-31    21.118481
2022-08-31    16.636031
2022-09-30    17.682613
2022-10-31    21.163298
2022-11-30    22.455672
2022-12-31    26.919529
2023-01-31    33.964623
2023-02-28    37.459776
2023-03-31    40.793766
2023-04-30    43.838415
2023-05-31    41.301780
2023-06-30    39.217866
2023-07-31    35.125502
2023-08-31    33.841926
2023-09-30    38.770511
2023-10-31    37.371216
2023-11-30    46.587633
2023-12-31    46.403326
Freq: M, Name: timeseries, 
dtype: float64

print(result.trend)

[Out]

2020-01-31          NaN
2020-02-29          NaN
2020-03-31          NaN
2020-04-30          NaN
2020-05-31          NaN
2020-06-30          NaN
2020-07-31    11.994971
2020-08-31    13.697685
2020-09-30    15.611236
2020-10-31    17.552031
2020-11-30    19.133760
2020-12-31    20.378094
2021-01-31    21.395429
2021-02-28    22.173782
2021-03-31    22.909215
2021-04-30    23.658616
2021-05-31    24.475426
2021-06-30    25.402005
2021-07-31    26.124056
2021-08-31    26.510640
2021-09-30    26.802553
2021-10-31    26.934270
2021-11-30    26.812893
2021-12-31    26.549220
2022-01-31    26.256876
2022-02-28    25.965319
2022-03-31    25.756854
2022-04-30    25.700551
2022-05-31    25.675143
2022-06-30    25.668984
2022-07-31    25.945528
2022-08-31    26.442986
2022-09-30    26.878992
2022-10-31    27.650819
2022-11-30    28.690242
2022-12-31    29.686565
2023-01-31    30.796280
2023-02-28    32.096818
2023-03-31    33.692393
2023-04-30    35.246385
2023-05-31    36.927214
2023-06-30    38.744537
2023-07-31          NaN
2023-08-31          NaN
2023-09-30          NaN
2023-10-31          NaN
2023-11-30          NaN
2023-12-31          NaN 

Freq: M, Name: timeseries,

dtype: float64 



 Seasonality 

 Residual (Noise)

print(result.seasonal)

[Out]:
2020-01-31    1.501630
2020-02-29    5.701170
2020-03-31    8.768065
2020-04-30    6.531709
2020-05-31    5.446174
2020-06-30    2.002476
2020-07-31   -1.643064
2020-08-31   -6.011071
2020-09-30   -7.807785
2020-10-31   -5.858728
2020-11-30   -5.849710
2020-12-31   -2.780867
2021-01-31    1.501630
2021-02-28    5.701170
2021-03-31    8.768065
2021-04-30    6.531709
2021-05-31    5.446174
2021-06-30    2.002476
2021-07-31   -1.643064
2021-08-31   -6.011071
2021-09-30   -7.807785
2021-10-31   -5.858728
2021-11-30   -5.849710
2021-12-31   -2.780867
2022-01-31    1.501630
2022-02-28    5.701170
2022-03-31    8.768065
2022-04-30    6.531709
2022-05-31    5.446174
2022-06-30    2.002476
2022-07-31   -1.643064
2022-08-31   -6.011071
2022-09-30   -7.807785
2022-10-31   -5.858728
2022-11-30   -5.849710
2022-12-31   -2.780867
2023-01-31    1.501630
2023-02-28    5.701170
2023-03-31    8.768065
2023-04-30    6.531709
2023-05-31    5.446174
2023-06-30    2.002476
2023-07-31   -1.643064
2023-08-31   -6.011071
2023-09-30   -7.807785
2023-10-31   -5.858728
2023-11-30   -5.849710
2023-12-31   -2.780867
Freq: M, Name: timeseries, 
dtype: float64

print(result.resid)

[Out]:
2020-01-31         NaN
2020-02-29         NaN
2020-03-31         NaN
2020-04-30         NaN
2020-05-31         NaN
2020-06-30         NaN
2020-07-31    3.776179
2020-08-31    4.256699
2020-09-30    1.924644
2020-10-31    0.790186
2020-11-30   -1.142242
2020-12-31   -2.453893
2021-01-31   -1.122544
2021-02-28    0.557940
2021-03-31    0.673303
2021-04-30    0.406231
2021-05-31    2.588922
2021-06-30    3.021039
2021-07-31   -0.180034
2021-08-31   -0.048653
2021-09-30   -0.123887
2021-10-31    0.250769
2021-11-30    1.939265
2021-12-31    2.852225
2022-01-31   -0.132007
2022-02-28    0.192434
2022-03-31    1.405550
2022-04-30   -2.054390
2022-05-31   -1.105152
2022-06-30   -1.079730
2022-07-31   -3.183983
2022-08-31   -3.795884
2022-09-30   -1.388594
2022-10-31   -0.628793
2022-11-30   -0.384861
2022-12-31    0.013830
2023-01-31    1.666713
2023-02-28   -0.338212
2023-03-31   -1.666692
2023-04-30    2.060321
2023-05-31   -1.071608
2023-06-30   -1.529146
2023-07-31         NaN
2023-08-31         NaN
2023-09-30         NaN
2023-10-31         NaN
2023-11-30         NaN
2023-12-31         NaN 

Freq: M, Name: timeseries,

dtype: float64 




# export to csv file

df.to_csv('ts_components.txt', sep=',', index=False)

 





  (2) R을 이용한 시계열 분해 (Time series Decomposition using R)


위에서 가법 모형을 적용해서 Python으로 만든 시계열 자료를 text 파일로 내보낸 후, 이를 R에서 읽어서 시계열 분해 (time series decomposition)를 해보겠습니다. 


ts_components.txt



# read time series text file

df <- read.table('ts_components.txt', sep=',', header=T)

head(df)

A data.frame: 6 × 6
timeseriestrendcycletrend_cycleseasonalirregular
<dbl><dbl><dbl><dbl><dbl><dbl>
2.5961190.00.0000000.0000000.0000002.5961193
6.7461601.11.3321982.4321983.5656840.7482781
8.1121002.22.6406474.8406476.136825-2.8653712
8.2559413.33.9020217.2020216.996279-5.9423582
16.8896554.45.0938349.4938345.9043271.4914946
16.1823575.56.19483911.6948393.1655361.3219812

 



이렇게 불러와서 만든 df DataFrame의 칼럼 중에서 시계열 분해를 할 'timeseries' 칼럼만을 가져와서 ts() 함수를 사용하여 1년 12개월 이므로 frequency =  12로 설정해 R의 시계열 자료 형태로 변환합니다. 


그 다음에 decompose() 함수를 사용하여 시계열 분해를 하는데요, 이때 가법 모형 (additive model)을 적용할 것이므로 decompose(ts, "additive") 라고 설정해줍니다. 


시계열 분해를 한 결과를 모아놓은 리스트 ts_decompose 객체를 프린트해보면 원래의 값 $x, 계절 요인 $seasonal, 추세(&순환) 요인 $trend, 불규칙 요인 $random 분해값이 순서대로 저장되어 있음을 알 수 있습니다. 



# transforming data to time series with 12 months frequency

ts <- ts(df$timeseries, frequency = 12) # 12 months


# time series decomposition

ts_decompose <- decompose(ts, "additive") # additive model


# decomposition results

ts_decompose

[Out]:

$x Jan Feb Mar Apr May Jun Jul 1 2.596119 6.746160 8.112100 8.255941 16.889655 16.182357 14.128087 2 21.774516 28.432892 32.350583 30.596556 32.510523 30.425519 24.300958 3 27.626499 31.858923 35.930469 30.177870 30.016165 26.591729 21.118481 4 33.964623 37.459776 40.793766 43.838415 41.301780 39.217866 35.125502 Aug Sep Oct Nov Dec 1 11.943313 9.728095 12.483489 12.141808 15.143334 2 20.450917 18.870881 21.326310 22.902448 26.620578 3 16.636031 17.682613 21.163298 22.455672 26.919529 4 33.841926 38.770511 37.371216 46.587633 46.403326 $seasonal Jan Feb Mar Apr May Jun Jul 1 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064 2 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064 3 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064 4 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064 Aug Sep Oct Nov Dec 1 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867 2 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867 3 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867 4 -6.011071 -7.807785 -5.858728 -5.849710 -2.780867 $trend Jan Feb Mar Apr May Jun Jul Aug 1 NA NA NA NA NA NA 11.99497 13.69769 2 21.39543 22.17378 22.90922 23.65862 24.47543 25.40200 26.12406 26.51064 3 26.25688 25.96532 25.75685 25.70055 25.67514 25.66898 25.94553 26.44299 4 30.79628 32.09682 33.69239 35.24639 36.92721 38.74454 NA NA Sep Oct Nov Dec 1 15.61124 17.55203 19.13376 20.37809 2 26.80255 26.93427 26.81289 26.54922 3 26.87899 27.65082 28.69024 29.68657 4 NA NA NA NA $random Jan Feb Mar Apr May Jun 1 NA NA NA NA NA NA 2 -1.12254398 0.55793990 0.67330316 0.40623129 2.58892205 3.02103851 3 -0.13200705 0.19243403 1.40555039 -2.05439035 -1.10515213 -1.07973023 4 1.66671294 -0.33821202 -1.66669164 2.06032097 -1.07160802 -1.52914637 Jul Aug Sep Oct Nov Dec 1 3.77617915 4.25669908 1.92464365 0.79018590 -1.14224232 -2.45389325 2 -0.18003389 -0.04865275 -0.12388741 0.25076875 1.93926482 2.85222467 3 -3.18398336 -3.79588443 -1.38859434 -0.62879275 -0.38486059 0.01383049 4 NA NA NA NA NA NA $figure [1] 1.501630 5.701170 8.768065 6.531709 5.446174 2.002476 -1.643064 [8] -6.011071 -7.807785 -5.858728 -5.849710 -2.780867 $type [1] "additive" attr(,"class") [1] "decomposed.ts"

 




위의 분해 결과가 숫자만 잔뜩 들어있으니 뭐가 뭔지 잘 눈에 안들어오지요?  그러면 이제 원래의 값 (observed)과 시계열 분해된 결과인 trend, seasonal, random 을 plot() 함수를 사용하여 다같이 시각화해보겠습니다. 



# change plot in jupyter

library(repr)


# Change plot size to 12 x 10

options(repr.plot.width=12, repr.plot.height=10)


plot(ts_decompose)

 



위의 분해 결과를 Trend (ts_decompose$trend), Seasonal (ts_decompose$seasonal), Random (ts_decompose$random) 의 각 요소별로 나누어서 시각화해볼 수도 있습니다. 



# change the plot size

options(repr.plot.width=12, repr.plot.height=5)


# Trend

plot(as.ts(ts_decompose$trend))

# Seasonality

plot(as.ts(ts_decompose$seasonal))

# Random (Irregular factor)

plot(as.ts(ts_decompose$random))




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



Posted by R Friend R_Friend

댓글을 달아 주세요

Lag, Lead window function은 시계열 데이터를 처리할 때 많이 사용하는 매우 유용한 함수입니다. 


이번 포스팅에서는 PostgreSQL, Python (pandas), R (dplyr) 을 이용해서 그룹별로 행을 하나씩 내리기, 올리기 (lag or lead a row by group using PostgreSQL, Python, R) 하는 방법을 소개하겠습니다. 





  1. PostgreSQL로 그룹별로 특정 칼럼의 행을 하나씩 내리기, 올리기 

    (lag, lead a row by group using PostgreSQL lag(), lead() window function)


연월일(dt), 그룹ID(id), 측정값(val) 의 세 개 칼럼을 가진 시계열 데이터의 테이블을 PostgreSQL DB에 만들어보겠습니다. 



DROP TABLE IF EXISTS ts;

CREATE TABLE ts (

    dt date not null

    , id text not null  

    , val numeric not null

);


INSERT INTO ts VALUES 

  ('2019-12-01', 'a', 5)

, ('2019-12-02', 'a', 6)

, ('2019-12-03', 'a', 7)

, ('2019-12-04', 'a', 8)

, ('2019-12-01', 'b', 13)

, ('2019-12-02', 'b', 14)

, ('2019-12-03', 'b', 15)

, ('2019-12-04', 'b', 16);


SELECT * FROM ts ORDER BY id, dt;




PostgreSQL 의 LAG(value, offset, default), LEAD(value, offset, default) Window function을 이용해서 그룹ID('id') 별로 측정값('val')의 행을 하나씩 내리기(lag), 올리기(lead) 해보겠습니다. 행을 내리거나 올린 후에 빈 셀의 값은 'NULL'로 지정해주었습니다. 


LAG(), LEAD() 함수를 사용할 때 그룹ID('id')별로 년월일('dt') 을 기준으로 내림차순 정렬(OVER(PARTITIO BY id ORDER BY dt)) 을 해줍니다. 



-- lead() windows function

SELECT 

    *

    , LAG(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lag_1

    , LEAD(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lead_2

FROM ts;

 



lag(), lead() 함수를 사용해서 lag_1, lead_2 라는 새로운 칼럼을 추가한 'ts_lag_lead' 라는 이름의 테이블을 만들어보겠습니다. 



DROP TABLE IF EXISTS ts_lag_lead;

CREATE TABLE ts_lag_lead AS (

SELECT 

    *

    , LAG(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lag_1

    , LEAD(val, 1, NULL) OVER (PARTITION BY id ORDER BY dt) AS val_lead_2

FROM ts

);


SELECT * FROM ts_lag_lead ORDER BY id, dt;

 





  2. Python pandas 로 DataFrame 내 그룹별 특정 칼럼의 행을 하나씩 내리기, 올리기 

     (shift a row by group using Python pandas library)


위에서 PostgreSQL의 lag(), lead() window function과 똑같은 작업을 Python pandas 를 가지고 수행해보겠습니다. 


먼저 dt, id, val의 칼럼을 가진 pandas DataFrame 시계열 데이터를 만들어보겠습니다. 



import pandas as pd


ts = pd.DataFrame({'dt': ['2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04', 

                          '2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04'], 

                  'id': ['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b'], 

                  'val': [5, 6, 7, 8, 13, 14, 15, 16]})


ts

dtidval
02019-12-01a5
12019-12-02a6
22019-12-03a7
32019-12-04a8
42019-12-01b13
52019-12-02b14
62019-12-03b15
72019-12-04b16

 



shift() 함수를 쓰기 전에 sort_values() 함수로 정렬을 해주는데요, lag 는 내림차순 정렬, lead는 오름차순 정렬임에 주의해야 합니다. (PostgreSQL, R 대비 Python이 좀 불편하긴 하네요 -,-;)


(a) lagsort_values() 함수를 이용해서 년월일('dt')를 기준으로 내림차순 정렬 (ascending=True) 한 후, 'id' 그룹별로 'val' 값을 하나씩 내려기 groupby('id')['val'].shift(1)


(b) lead: sort_values() 함수를 이용해서 년월일('dt')를 기준으로 오름차순 정렬 (ascending=False) 한 후, 'id' 그룹별로 'val' 값을 하나씩 올리기 groupby('id')['val].shift(1)



# lag a row by group 'id'

ts['val_lag_1'] =  ts.sort_values(by='dt', ascending=True).groupby('id')['val'].shift(1)


# lead a row by group 'id'

ts['val_lead_1'] = ts.sort_values(by='dt', ascending=False).groupby('id')['val'].shift(1)


ts.sort_values(by=['id', 'dt'])

dtidvalval_lag_1val_lead_1
02019-12-01a5NaN6.0
12019-12-02a65.07.0
22019-12-03a76.08.0
32019-12-04a87.0NaN
42019-12-01b13NaN14.0
52019-12-02b1413.015.0
62019-12-03b1514.016.0
72019-12-04b1615.0NaN

 





  3. R dplyr 로 dataframe 내 그룹별 특정 칼럼의 행을 하나씩 내리기, 올리기 

     (lag, lead a row by group using R dplyr library)



위에서 PostgreSQL의 lag(), lead() window function과 똑같은 작업을 R dplyr library를 가지고 수행해보겠습니다. 


먼저 dt, id, val의 칼럼을 가진 R DataFrame 시계열 데이터를 만들어보겠습니다. 



#install.packages("dplyr")

library(dplyr)


dt <- c(rep(c('2019-12-01', '2019-12-02', '2019-12-03', '2019-12-04'), 2))

id <- c(rep('a', 4), rep('b', 4)) 

val <- c(5, 6, 7, 8, 13, 14, 15, 16)


ts <- data.frame(dt, id, val)

ts

A data.frame: 8 × 3
dtidval
<fct><fct><dbl>
2019-12-01a5
2019-12-02a6
2019-12-03a7
2019-12-04a8
2019-12-01b13
2019-12-02b14
2019-12-03b15
2019-12-04b16

 



R은 Postgresql 처럼 lag(), lead() window function을 가지고 있고 dplyr library의 chain operator를 써서 arrange() 함수로 'dt' 기준 내림차순 정렬하고, group_by(id)를 써서 그룹ID('id')별로 lag(), lead()를 무척 편리하게 적용해서 새로운 변수를 생성(mutate)할 수 있습니다. 



ts <- ts %>% 

    arrange(dt) %>%

    group_by(id)  %>% 

    mutate(val_lag_1 = lag(val, 1), 

          val_lead_1 = lead(val, 1))

 


arrange(ts, id, dt)

A grouped_df: 8 × 5
dtidvalval_lag_1val_lead_1
<fct><fct><dbl><dbl><dbl>
2019-12-01a5NA6
2019-12-02a657
2019-12-03a768
2019-12-04a87NA
2019-12-01b13NA14
2019-12-02b141315
2019-12-03b151416
2019-12-04b1615NA




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요.



Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 2019.12.09 13:32  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

  2. 2019.12.09 19:36  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend R_Friend 2019.12.12 11:06 신고  댓글주소  수정/삭제

      안녕하세요.

      제가 지금 해외출장 나와있는데요, 어찌된 일인지 모바일 티스토리 앱으로 댓글 업데이트 알람을 받지 못해서요, 이제서야 (12월12일, 오전10시) 웹으로 블로그 로그인 한 후에야 추가로 댓글 달아주신것을 확인했습니다.

      제가 업무시간에는 시간을 낼 수 없어서요, 퇴근 후에 저녁에 답글 달도록 하겠습니다. 너무나 많이 늦어져서 죄송합니다.

    • 2019.12.12 20:35  댓글주소  수정/삭제

      비밀댓글입니다

    • R Friend R_Friend 2019.12.15 03:01 신고  댓글주소  수정/삭제

      안녕하세요 차차님,

      댓글에 남겨주신 말씀처럼, pandas 시계열 데이터 DataFrame 을 처리할 때 특정 시간 단위 구간별 집계/요약은 resample() 메소드를 사용하는 것이 가장 간편한 방법입니다.

      resample() 메소드를 사용해서 여러 통계량을 좀더 쓰기 쉽도록 사용자 정의 함수(UDF)를 만들어보았는데요, 아래의 포스팅을 참고하시기 바랍니다.

      https://rfriend.tistory.com/494


      만약, 대용량 시계열 데이터를 처리해야 하는 경우라면 관계형 Database를 사용하는 것이 빠르고 편리한 방법입니다.
      PostgreSQL DB를 사용해서 python pandas resample() 메소드와 동일한 결과를 얻기 위한 방법에 대한 소개는 아래 포스팅을 참고하세요.

      https://rfriend.tistory.com/495

      포스팅 2개 하다보니 주말이 훌쩍 가버리네요. ^^;

이번 포스팅에서는 STEFANY COXE, STEPHEN G. WEST, AND LEONA S. AIKEN, 2009, "The Analysis of Count Data: A Gentle Introduction to Poisson Regression and Its Alternatives" 논문을 요약해서 한글로 번역한 내용을 소개하겠습니다. 


빈도(count data) 데이터에 대한 분석을 할 때 어떤 분포를 가정하고, 어떤 회귀모형을 적합하는 것이 좋을지 분석 방향을 잡는데 도움이 될 것입니다. 






  왜 포아송 회귀모형이 필요한가(Why Poisson Regression)?


특정 시간 동안 발생한 사건의 건수에 대한 도수 자료(count data, 음수가 아닌 정수)를 목표변수, 특히 평균이 10 미만일 경우, 최소제곱법(Ordinary Least Squares) 회귀모형을 적합하면 표준 오차와 유의수준이 편향되는 문제가 발생한다. OLS 회귀모형은 오차의 조건부 정규성(conditional normality), 등분산성(homoscedasiticity), 독립성(independence)을 가정한다. 하지만 평균이 작은 도수 데이터는 0 이상의 정수만 있고, 작은 계급에 많은 관측치가 몰려있으며, 우측으로 길게 치우친 분포를 띠어 회귀모형의 가정을 위배하고 음수값도 결과값으로 반환하는 문제가 있다. 종속변수가 도수 데이터이면서 오차가 정규분포, 등분산이 아닌 경우 포아송 회귀모형을 사용한다.

 


  • 포아송 회귀모형은 무엇인가? (Overview of Poisson Regression)

포아송 회귀모형은 일반화선형모형(Generalized Linear Model, GLM)의 한 종류로서,

-       Random component: Poisson Distribution, 


-       Systematic component: Mixed linear relationship

-       Link function: Log Link (ln)

 

표현 예

해석  이며, X1을 한 단위 증가시켰을 때 효과를 보면   이므로, 다른 변수들을 통제했을 때 X1이 한단위 증가하면 도수의 추정치 만큼 승법적으로 변화(multiplicative changes)한다.


선형회귀모형이 최소제곱법(OLS)로 모수를 추정한다면, 포아송 회귀모형은 최대가능도추정(MLE, Maximum Likelihood Estimation)을 통해 모수를 추정한다.


선형회귀모형은 모델에 의해 설명되는 잔차 제곱합 SS (Sum of Squares)를 총 잔차제곱합(TSS, Total Sum of Squares)로 나눈 값인  R^2로 모델 적합도(Goodness of Fit of the model) 평가하지만, 포아송 회귀모형은 설명변수의 추가에 따른 이탈도의 감소 비율(proportional reduction in deviance)로 표현되는 pseudo-R^2를 사용하여 완전모델(perfect model)에 얼마나 가깝게 적합이 되었는지를 평가한다



선형회귀모형은 F-test를 통해 모델의 유의성을 검정하지만, 포아송 회귀모형은 base modelcomplete model 간의 이탈도 차이에 대한 Chi-square test를 통해 모델 유의성을 검정한다 

선형회귀모형은 모형 타당성 평가 (Assessing model adequacy)를 위해 잔차도(residual plot)를 그려서 확인해보는 반면, 포아송 회귀모형은 이탈도 잔차도(deviance residuals vs. predicted values)를 그려보거나 실측치와 예측치를 비교해보는 것이다.

 

* [참고] 포아송 분포 : https://rfriend.tistory.com/101




  과대산포(Over-dispersion) 문제는 무엇이며어떻게 해결할 수 있나?


포아송분포는 평균과 분산이 동일(equidispersion)하다고 가정하므로, 분산이 평균보다 큰 과대산포 (Overdispersion) 시에는 포아송 회귀모형을 사용할 수 없다. 과대산포 시에 적용할 수 있는 대안 모델로는 과대산포 포아송 회귀(Overdispersed Poisson Regression), 음이항회귀(Negative Binomial Regression)가 있다.


-   Overdispersed Poisson Regression: 과대산포를 모델링하기 위해 두번째 모수인 조건부 분산 overdispersion scaling parameter 를 추정하며, 과대산포 포아송 회귀모형의 오차 분포는  를 가진다.


-   Negative Binomial Regression Models: 음이항회귀모형은 평균에는 영향이 없으면서 과대산포를 유발하는, 설명이 되지 않는 추가적인 가변성이 있다고 가정한다. 음이항 회귀모형은 똑같은 설명변수 값을 가지는 관측치가 다른 평균 모수를 가지고 포아송 회귀모형에 적합될 수 있도록 해준다. 오차 함수는 포아송 분포와 감마분포의 혼합 분포이다.


과대산포 검정은 가능도비 검정(Likelihood ratio test) 또는 score 검정을 사용한다.

모델 선택을 위해서는 AIC(Akaike Information Criterion), BIC(Bayesian Information Criterion) 을 사용하며, AIC, BIC의 값이 작은 모델을 선택하는데, 샘플 크기가 작을 때는 간단한 모델을 선택하고, 샘플 크기가 클 때는 좀더 복잡한 모델을 선택한다.

 



  과대영(Excess Zeros) 문제는 무엇이며어떻게 해결할 수 있나?


특정 평균을 모수로 가지는 포아송 분포에서 나타나는 ‘0’보다 많은 ‘0’을 가지는 샘플의 경우 과대영(excess zeros)이 발생했다고 말한다. 과대영은 (a) 특정 행동을 나타내지 않은 그룹의 포함 (: 비흡연자 그룹), (b) 특정 행동을 나타내지만 샘플링 계획에서 제외 (: 병원 입원한 사람만 연구) 하는 경우에 발생할 수 있다. 과대영일 경우 ZIP 회귀모형을 대안으로 사용할 수 있다.


-   Zero-Inflated Poisson(ZIP) Regression: 두 부분으로 나뉘는데, (1) 로지스틱 회귀모형을 사용하여 항상 ‘0’인 집단(structural zeros) 또는 항상 ‘0’은 아니지만 조사 시점에 ‘0’이라고 응답한 집단에 속할 확률을 구하고, (2) 포아송 회귀 또는 음이항 회귀모형으로 structural zeros를 포함하지 않은 나머지 관측치를 추정하는 모델을 적합한다.





  음이항분포(Negative Binomial Distribution) 


* Reference : https://en.wikipedia.org/wiki/Negative_binomial_distribution


 
확률 이론과 통계에서, 음 이항 분포는 특정 (무작위) 실패 횟수 (r로 표시)가 발생하기 전에 독립적이고 동일하게 분산 된 베르누이 시행 순서에서 성공 횟수의 이산 확률 분포이다. 예를 들어 1을 실패로, 1이 아닌 모든 것을 성공으로 정의하고 1이 세 번째로 나타날 때까지 반복적으로 던지면 (r = 3 번의 실패), 출현한 1이 아닌 수(non-1)의 확률 분포는 음 이항 분포가 된다.


 파스칼 분포 (Blaise Pascal 이후) Polya 분포 (George Pólya의 경우)는 음 이항 분포의 특수한 경우이다. 엔지니어, 기후 학자 및 다른 사람들 사이의 협약은 정수 값 정지 시간 매개 변수 r의 경우 "음 이항" 또는 "파스칼"을 사용하고, 실수 값의 경우에는 "Polya"를 사용한다. 토네이도 발생과 같은 "전염성 있는" 이산 이벤트가 발생하는 경우 Polya 분포를 사용하여 Poisson과 달리 평균 및 분산을 다르게 하여 Poisson 분포보다 더 정확한 모델을 제공 할 수 있다. "전염성 있는" 사건은 양의 공분산으로 인해 사건이 독립적이었던 경우보다 양의 상관관계가 있어 Poisson 분포에 의한 분산보다 더 큰 분산을 일으킨다. Poisson 분포는 평균과 분산이 동일하다고 가정하기 때문에 평균보다 분산이 큰 분포(Overdispersed)의 경우 적절하지 않으며, 2개의 모수를 가져서 과대산포분포를 적합할 수 있는 음 이항 분포를 사용할 수 있다.


 일련의 독립적 인 베르누이 재판이 있다고 가정하자. 따라서 각 임상 시험은 "성공" "실패"라는 두 가지 결과를 나타낸다. 각 시험에서 성공 확률은 p이고 실패 확률은 (1 - p)이다. 미리 정의 된 실패 횟수 r이 발생할 때까지의 서열을 관찰한다. 그러면 우리가 본 성공의 무작위 수 X는 음 이항 (또는 파스칼) 분포를 가질 것이다





 R을 활용한 Poisson model, Negative binomial model, Zero-Inflated model


Poisson model, Negative binomial model, Zero-Inflated model 의 유형, 분포, 모수 추정 방법, R패키지와 함수를 표로 정리하면 아래 [1]을 참고하세요. 


[1] Overview of count regression models


유형

분포

추정

방법

Description

R package & function

GLM

Poisson

ML

Poisson Regression: classical GLM, estimated by maximum likelihood(ML)

{stats} 패키지의 glm()

NB

ML

NB Regression: exteded GLM, estimated by ML including additional shape parameter

{stats} 패키지의 glm(), {MASS} 패키지의 glm.nb() 함수

Zero-augmented

Poisson

ML

Zero-Inflated Poission(ZIP)

{pscl} 패키지의 zeroinfl() 함수

 



각 회귀모형별로 R 패키지와 함수의 활용 방법은 아래와 같다.


-       (1) R을 활용한 Poisson model : glm() in the {stats} package



glm(formula, data, subset, na.action, weights, offset,

    family = poisson, start = NULL, control = glm.control(...),

    model = TRUE, y = TRUE, x = FALSE, ...)


 



-       (2) R을 활용한 Negative binomial model
   :
glm() in the {stats} package



glm(formula, data, subset, na.action, weights, offset,

    family = negative.binomial, start = NULL, control = glm.control(...),

    model = TRUE, y = TRUE, x = FALSE, ...)


 


   : glm.nb() in the {MASS} package



glm.nb(formula, data, weights, subset, na.action,

       start = NULL, etastart, mustart,

       control = glm.control(...), method = "glm.fit",

       model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, ...,

       init.theta, link = log)


 



-       (3) R을 활용한 Zero-inflated regression model: zeroinfl() in the {pscl} package



zeroinfl(formula, data, subset, na.action, weights, offset,

dist = "poisson", link = "logit", control = zeroinfl.control(...),

model = TRUE, y = TRUE, x = FALSE, ...)




많은 도움이 되었기를 바랍니다. 



Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 2019.10.27 15:26  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

    • R Friend R_Friend 2019.10.27 16:45 신고  댓글주소  수정/삭제

      안녕하세요 Mini님,

      GAM은 GLM을 포괄하는, smoother 조정을 통해 좀더 flexible하고 nonlinear한 적합이 가능한 반면, 말씀하신대로 non-parametric approach로서 RR, CI 등은 알 수 없습니다.

      만약 RR, CI 등을 구해야 하는 상황이라면 parametric approach (모델 기반)인 GLM을 사용하시는게 맞을거 같습니다.

  2. mini 2019.10.27 21:09  댓글주소  수정/삭제  댓글쓰기

    답변 너무 감사합니다.
    exp()로 RR은 구했는데요. 95% 신뢰구간을 구하는 방법을 모르겠습니다.
    아래와 같은 방법으로 구하는 것이 맞을까요?? 구글링을 많이 해봤는데, 딱 맞는 것은 못찾아 문의 드리게 되었습니다.

    > exp(-0.0041)
    [1] 0.9959084
    > 0.9959084-(1.96*0.0007)
    [1] 0.9945364
    > 0.9959084+(1.96*0.0007)
    [1] 0.9972804

  3. 2019.10.27 21:53  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

  4. Mini 2019.11.01 09:03  댓글주소  수정/삭제  댓글쓰기

    Rr과 신뢰구간은 모두 구했습니다~
    Exp(coef(모델))
    EXP(confint(모델))로요
    도움 감사합니다.
    그려진 플롯에서 변곡점을 찾아야되는데요.
    아무리 인터넷을 찾아봐도
    gam에서 적용할만한 내용이 안나오네요.
    실례가 안된다면 여쭤봐도 될까요?

  5. mini 2019.11.02 18:21  댓글주소  수정/삭제  댓글쓰기

    감사합니다.
    GAm에서는 안되는 것인지, 제가 못하는 것인지 확실치 않은데요. 안되네요.
    그래도, 도움 주셔서 너무 감사합니다. !

R의 다양한 패키지들을 사용하다 보면 함수를 사용한 후의 산출물 결과가 자신이 생각했던 결과와 다를 경우 R 패키지 함수의 소스 코드가 어떻게 되어있는 것인지 궁금할 때가 있습니다. 


이번 포스팅에서는 R 패키지에 있는 함수, 메소드의 소스 코드를 들여다볼 수 있는 2가지 방법을 소개하겠습니다. 


(1) getAnywhere(function_name)

(2) package_name:::function_name





가령, 두 벡터의 차집합(difference of sets)을 구할 때 사용하는 함수가 base 패키지의 setdiff() 함수입니다. 


> x <- c(1, 1, 2, 3, 4, 5)

> y <- c(4, 5, 6, 7)

> setdiff(x, y)

[1] 1 2 3




위에서 예로 든 base 패키지의 setdiff() 함수에 대한 소스 코드를 위의 2가지 방법을 사용해서 살펴보겠습니다. 



 (1) getAnywhere(function_name)


R package 의 이름을 몰라도 함수, 메소드의 이름만 알고 있으면 R package 이름과 namespace, 그리고 소스 코드를 친절하게 다 볼 수 있으므로 매우 편리합니다. 



> getAnywhere(setdiff)

A single object matching ‘setdiff’ was found

It was found in the following places

  package:base

  namespace:base

with value


function (x, y) 

{

    x <- as.vector(x)

    y <- as.vector(y)

    unique(if (length(x) || length(y)) 

        x[match(x, y, 0L) == 0L]

    else x)

}

<bytecode: 0x000000000ff1e8d0>

<environment: namespace:base>

 




 (2) R_package_name:::function_name


R package 이름을 알고 있다면 package_name:::function_name 의 형식으로 ':::' 를 사용해서 함수의 소스 코드를 들여다볼 수 있습니다. 위 (1)번의 getAnywhere() 함수가 기억이 잘 안날 때 쉽게 ':::' 를 사용해서 소스 코드 볼 때 편리합니다. 



> base:::setdiff

function (x, y) 

{

    x <- as.vector(x)

    y <- as.vector(y)

    unique(if (length(x) || length(y)) 

        x[match(x, y, 0L) == 0L]

    else x)

}

<bytecode: 0x000000000ff1e8d0>

<environment: namespace:base>




많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



Posted by R Friend R_Friend

댓글을 달아 주세요

이번 포스팅에서는 (1) 시계열 정수의 순차 3개 묶음 패턴 별 개수를 구하고, (2) 각 패턴별로 발생 빈도 기준으로 내림차순 정렬하는 방법을 소개하겠습니다. 



먼저 S = {1, 2, 3, 4, 5} 의 정수 집합에서 복원추출(replacement) 하여 무작위 샘플 1,000개를 생성해보겠습니다. 



> #----------------------------------

> # counting and sorting by patterns

> #----------------------------------

> set.seed(123) # for reproducibility

> sample <- sample(x=1:5, size=1000, replace=TRUE)

> sample

   [1] 2 4 3 5 5 1 3 5 3 3 5 3 4 3 1 5 2 1 2 5 5 4 4 5 4 4 3 3 2 1 5 5 4 4 1 3 4 2 2 2 1 3 3 2 1 1 2

  [48] 3 2 5 1 3 4 1 3 2 1 4 5 2 4 1 2 2 5 3 5 5 4 3 4 4 4 1 3 2 2 4 2 1 2 4 3 4 1 3 5 5 5 1 1 4 2 4

  [95] 2 1 4 1 3 3 3 2 3 5 3 5 5 4 3 1 5 2 1 5 4 1 3 5 3 3 4 2 2 2 2 5 1 1 1 4 4 5 4 4 3 4 5 4 5 3 2

 [142] 3 1 1 5 2 2 1 2 4 5 3 2 2 1 2 3 2 3 2 3 2 4 2 2 3 4 2 3 2 4 1 5 4 4 4 2 3 5 3 5 2 4 2 3 3 2 3

 [189] 5 5 2 2 5 4 5 3 3 4 1 3 2 5 4 3 3 5 2 2 1 1 3 2 2 4 1 4 2 3 5 5 2 5 4 4 1 2 3 3 4 5 4 3 3 1 2

 [236] 2 1 5 1 5 3 4 1 4 2 4 2 5 5 4 2 2 3 2 3 4 1 3 3 5 5 5 4 5 3 3 2 2 1 3 5 1 1 1 4 4 5 3 1 4 4 1

 [283] 2 2 1 2 1 2 1 4 2 1 1 5 4 5 5 1 1 4 4 1 4 4 4 3 1 1 3 3 2 3 4 1 2 5 5 2 2 5 5 2 1 4 1 1 5 1 2

 [330] 5 4 2 4 5 2 3 4 4 4 5 3 1 3 2 3 2 5 2 2 4 1 5 3 1 4 5 4 3 2 5 1 1 5 2 2 4 1 1 3 3 3 2 5 3 3 5

 [377] 4 2 2 5 3 4 2 4 3 5 1 2 1 2 5 2 4 4 1 2 2 1 1 5 5 1 5 3 2 3 4 1 2 4 2 5 2 3 2 1 5 2 3 1 3 4 3

 [424] 2 5 4 2 2 1 5 5 3 2 5 2 4 1 3 1 1 2 3 4 1 5 4 4 1 2 4 3 2 2 1 2 5 5 4 1 1 5 3 3 1 3 2 2 4 2 5

 [471] 1 5 3 4 4 3 4 1 2 4 2 5 2 4 5 1 1 4 3 5 5 1 3 4 1 5 2 1 4 4 2 2 2 1 2 1 3 3 5 2 3 1 5 2 4 1 3

 [518] 4 1 5 2 2 2 1 3 4 5 1 5 2 5 4 3 5 1 4 2 5 3 2 5 1 5 2 4 5 1 5 5 2 2 2 3 5 3 3 4 3 3 1 1 5 4 2

 [565] 4 4 2 5 3 5 3 5 3 3 5 4 3 1 3 3 5 5 5 2 5 1 2 4 5 3 3 1 5 1 1 5 3 2 5 1 2 4 2 2 1 5 1 5 2 2 4

 [612] 1 1 5 3 4 5 4 5 1 5 5 4 2 1 2 2 5 2 2 4 5 3 4 1 2 5 4 5 3 3 4 5 2 1 4 4 1 2 4 4 5 3 3 4 5 5 1

 [659] 3 2 5 5 1 4 2 1 5 1 4 1 1 3 1 2 3 3 5 2 3 3 4 4 5 5 5 3 3 1 3 2 1 4 3 1 2 3 1 1 3 2 5 2 1 5 2

 [706] 4 3 2 1 3 4 3 4 4 3 2 1 1 5 1 3 4 4 5 4 4 5 4 4 1 3 5 4 5 1 2 5 5 2 1 5 1 1 3 2 3 4 1 4 5 2 2

 [753] 1 1 3 3 3 5 1 2 3 2 3 3 3 1 1 3 4 4 5 5 4 5 1 2 2 3 4 5 2 5 1 4 5 3 2 1 2 5 5 3 2 3 5 3 1 2 2

 [800] 2 3 2 1 1 2 5 3 3 4 4 1 4 2 5 4 1 2 1 3 3 4 2 4 5 4 1 3 1 2 4 1 4 3 4 4 3 4 4 2 1 5 4 4 1 3 2

 [847] 5 1 2 1 2 4 5 4 3 3 5 1 1 5 3 1 4 3 4 2 1 1 3 3 1 1 4 3 5 4 1 4 5 4 1 5 4 1 4 5 3 2 2 2 3 4 4

 [894] 3 5 4 1 2 1 4 5 3 5 3 4 3 4 5 5 1 1 1 5 3 2 5 3 3 2 2 5 1 2 5 4 4 2 3 1 4 5 2 4 3 5 1 3 5 5 1

 [941] 5 3 2 3 1 1 3 2 3 2 3 1 5 1 2 5 4 1 3 3 3 3 4 1 5 3 3 1 4 3 4 1 5 2 2 4 5 1 2 3 3 4 5 4 5 4 1

 [988] 1 1 3 1 2 3 3 3 5 4 2 4 1

 




다음으로 순차적으로 3개의 정수를 묶음으로 하여 1개씩 이동하면서 패턴을 정의하고(패턴1 {X1, X2, X3}, 패턴2 {X2, X3, X4}, ..., 패턴n-2 {Xn-2, Xn-1, X}, 위의 정수 난수 샘플을 가지고 예를 들면, 첫번째 3개 묶음의 패턴은 '2_4_3', 두번째 3개 묶음의 패턴은 '4_3_5', 세번째 3개 묶음의 패턴은 '3_5_5', ... ), 각 패턴별로 발생 빈도를 세어보겠습니다


먼저 비어있는 pattern_list 라는 이름의 리스트(list) 를 만들어놓구요, 


for loop 반복문을 사용하여 '패턴1' {X1, X2, X3}, '패턴2' {X2, X3, X4}, ..., '패턴n-2' {Xn-2, Xn-1, X} 을 생성합니다. 


if else 조건문을 사용하여, 만약 3개 정수 묶음의 패턴이 pattern_list의 키(key)에 이미 들어있는 것이면 +1 을 추가해주구요, 그렇지 않다면 (처음 발생하는 패턴이라면) pattern_list의 키(key)에 새로 발생한 패턴을 키로 추가해주고 값(value)으로 1을 부여해줍니다. 



> # blank list to store the patterns' count

> pattern_list <- {}

> # count per patterns

> for (i in 1:(length(sample)-2)){

+   

+   pattern <- paste(sample[i], sample[i+1], sample[i+2], sep="_")

+   

+   if (pattern %in% names(pattern_list)) {

+     pattern_list[[pattern]] <- pattern_list[[pattern]] + 1

+   } else {

+     pattern_list[[pattern]] <- 1

+   }

+ }

> pattern_list

2_4_3 4_3_5 3_5_5 5_5_1 5_1_3 1_3_5 3_5_3 5_3_3 3_3_5 5_3_4 3_4_3 4_3_1 3_1_5 1_5_2 5_2_1 2_1_2 

    6     7    10     9     6     6    10    15    11     7     5     5     5    13     7    12 

1_2_5 2_5_5 5_5_4 5_4_4 4_4_5 4_5_4 4_4_3 4_3_3 3_3_2 3_2_1 2_1_5 1_5_5 4_4_1 4_1_3 1_3_4 3_4_2 

   11     8     9    11     9    13     7     5     8     9    10     5    11    14     9     6 

4_2_2 2_2_2 2_2_1 2_1_3 1_3_3 2_1_1 1_1_2 1_2_3 2_3_2 3_2_5 2_5_1 3_4_1 1_3_2 2_1_4 1_4_5 4_5_2 

    8     8    14     6    11     8     3     9    11    13    10    15    12     8     8     6 

5_2_4 2_4_1 4_1_2 1_2_2 2_2_5 2_5_3 5_3_5 5_4_3 4_3_4 3_4_4 4_4_4 3_2_2 2_2_4 2_4_2 4_2_1 1_2_4 

    9    10    12     7     7     7     6     8    10    11     4     8     8     9     7    10 

5_5_5 5_1_1 1_1_4 1_4_2 4_2_4 1_4_1 3_3_3 3_2_3 2_3_5 1_5_4 5_4_1 3_3_4 1_1_1 1_4_4 3_4_5 5_4_5 

    4    10     6     7     7     3     7    15     6     7    10    11     4     7     9     9 

4_5_3 5_3_2 2_3_1 3_1_1 1_1_5 5_2_2 2_4_5 3_2_4 2_2_3 2_3_4 4_2_3 4_1_5 4_4_2 3_5_2 2_3_3 5_5_2 

   13    11     7     8    12    12     9     2     6     9     5     8     5     4     7     7 

2_5_4 1_1_3 4_1_4 5_2_5 3_3_1 3_1_2 1_5_1 5_1_5 1_5_3 4_2_5 5_4_2 3_5_1 5_3_1 3_1_4 1_2_1 4_5_5 

   10    11     8     4     8     6     8     9    11     7     7     6     5     5     7     5 

4_1_1 5_1_2 5_2_3 3_1_3 2_5_2 4_3_2 3_5_4 2_4_4 5_5_3 1_3_1 4_5_1 1_4_3 5_1_4 

    6    11     5     5     7     5     6     3     3     4     7     6     4 

 




이제 발생 빈도를 기준으로 패턴 리스트를 내림차순 정렬(sort in descending order)을 해보겠습니다. 

'5_3_3', '3_4_1', '3_2-3' 패턴이 총 15번 발생해서 공동 1등을 하였네요. 



> # sorting pattern_list in descending order

> sort(pattern_list, decreasing = TRUE)

5_3_3  3_4_1  3_2_3  4_1_3 2_2_1 1_5_2 4_5_4 3_2_5 4_5_3 2_1_2 1_3_2 4_1_2 1_1_5 5_2_2 3_3_5 1_2_5 

   15     15      15     14    14    13    13    13    13    12    12    12    12    12    11    11 

5_4_4 4_4_1 1_3_3 2_3_2 3_4_4 3_3_4 5_3_2 1_1_3 1_5_3 5_1_2 3_5_5 3_5_3 2_1_5 2_5_1 2_4_1 4_3_4 

   11    11    11    11    11    11    11    11    11    11    10    10    10    10    10    10 

1_2_4 5_1_1 5_4_1 2_5_4 5_5_1 5_5_4 4_4_5 3_2_1 1_3_4 1_2_3 5_2_4 2_4_2 3_4_5 5_4_5 2_4_5 2_3_4 

   10    10    10    10     9     9     9     9     9     9     9     9     9     9     9     9 

5_1_5 2_5_5 3_3_2 4_2_2 2_2_2 2_1_1 2_1_4 1_4_5 5_4_3 3_2_2 2_2_4 3_1_1 4_1_5 4_1_4 3_3_1 1_5_1 

    9     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8 

4_3_5 5_3_4 5_2_1 4_4_3 1_2_2 2_2_5 2_5_3 4_2_1 1_4_2 4_2_4 3_3_3 1_5_4 1_4_4 2_3_1 2_3_3 5_5_2 

    7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7 

4_2_5 5_4_2 1_2_1 2_5_2 4_5_1 2_4_3 5_1_3 1_3_5 3_4_2 2_1_3 4_5_2 5_3_5 1_1_4 2_3_5 2_2_3 3_1_2 

    7     7     7     7     7     6     6     6     6     6     6     6     6     6     6     6 

3_5_1 4_1_1 3_5_4 1_4_3 3_4_3 4_3_1 3_1_5 4_3_3 1_5_5 4_2_3 4_4_2 5_3_1 3_1_4 4_5_5 5_2_3 3_1_3 

    6     6     6     6     5     5     5     5     5     5     5     5     5     5     5     5 

4_3_2 4_4_4 5_5_5 1_1_1 3_5_2 5_2_5 1_3_1 5_1_4 1_1_2 1_4_1 2_4_4 5_5_3 3_2_4 

    5     4     4     4     4     4     4     4     3     3     3     3     2

 




위의 패턴별 발생 빈도수 기준으로 정렬된 결과를 DataFrame으로 변환해서 상위 10개 패턴을 프린트해보겠습니다. 



> # convert a list to DataFrame

> cnt_per_pattern_sorted <- sort(pattern_list, decreasing = TRUE)

> pattern <- names(cnt_per_pattern_sorted)

> df <- data.frame(pattern, cnt_per_pattern_sorted)

> rownames(df) <- NULL

> # display top 10 patterns

> df[1:10,]

   pattern     cnt_per_pattern_sorted

1    5_3_3                     15

2    3_4_1                     15

3    3_2_3                     15

4    4_1_3                     14

5    2_2_1                     14

6    1_5_2                     13

7    4_5_4                     13

8    3_2_5                     13

9    4_5_3                     13

10   2_1_2                     12

 



많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. ^^



Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 최성원 2019.10.16 14:15  댓글주소  수정/삭제  댓글쓰기

    항상 감사드립니다. 질문이 있습니다!
    수집된 패턴들을 가지고 다른 분석을 하고 싶은데요, 시계열 군집분석이나 다변량 군집분석이요...
    그래서 데이터 타입을 "list" 형태로 변환해봤는데 에러가 나더라구요.
    제가 데이터 타입에 대한 깊은 이해가... 부족해서 그런것 같은데
    질문은...
    위의 방법으로 패턴 수집했을때 데이터가 어떤 형태로 저장되는 건가요?
    그리고 리스트나 데이터 프레임의 형태로 저장은 어떻게 하나요?

    • R Friend R_Friend 2019.10.16 14:22 신고  댓글주소  수정/삭제

      안녕하세요.

      현재 분석 결과는 리스트로 저장을 해둔 상태입니다.

      만약 패턴별 발생빈도 리스트를 데이터프레임으로 변경을 하고 싶으시면,

      df <- data.frame(matrix(unlist(pattern_list), nrow=length(pattern_list), byrow=T), stringAsFactors=FALSE)

      로 하시면 됩니다. (만약 character 를 factor type으로 해서 분석할 경우라면 stringAsFactors 옵션을 빼주면 됩니다)

      ps. 제가 이번주는 너무 바빠서 더이상은 못도와드리겠스니다. 죄송합니다. 이해해주세요. ㅜ.ㅜ

    • 최성원 2019.10.16 16:28  댓글주소  수정/삭제

      지금도 너무 감사드려요 ㅠㅠ

    • R Friend R_Friend 2019.10.16 16:29 신고  댓글주소  수정/삭제

      네, 이해해주셔서 고맙습니다. 아무쪼록 순탄하게 잘 해결하시기를 바래요.

    • R Friend R_Friend 2019.10.17 23:50 신고  댓글주소  수정/삭제

      본문 하단에 DataFrame으로 변환하는 코드도 update 했습니다. ^^

  2. 질문드립니다. 2019.10.24 10:16  댓글주소  수정/삭제  댓글쓰기

    매번 친절하게 알려주셔서 감사드립니다^^

    데이터 프레임에 제품 On/off 시간에 대해서 수집하고 있는데
    제품을 동작하고 끄기까지의 사용시간을 계산하고 싶은데 너무 초보라
    이해도가 부족해서 코딩을 하기가 힘들어서 문의드립니다..
    아래처럼 시작 시간 정지 시간 별로 시간 데이터가 입력이 되는 구조입니다.

    DEVICE | time | state | sum |
    ㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡㅡ
    A |2018-07-01 10:20:23| ON | |
    A |2018-07-01 10:30:55| ON | |
    A |2018-07-01 11:11:01| ON | |
    A |2018-07-01 11:51:41| OFF | |
    B |2018-07-02 07:11:02| ON | |
    B |2018-07-02 09:00:33| ON | |
    B |2018-07-02 09:20:24| OFF | |
    C |2018-07-02 12:12:21| ON | |
    C |2018-07-02 14:01:09| OFF | |
    C |2018-07-02 18:11:41| ON | |
    C |2018-07-02 19:21:51| OFF | |

    DEVICE 별로 시작과 정지 사이의 시간을 계산해서 합한 후에 열을 추가해서 기입하고
    A 제품처럼 ON과 OFF 사이에 ON 이 여러개 들어있으면 이것들을 제거하고 처음ON과
    마지막 OFF 사이만 계산하게끔 하고 싶습니다..ㅠ
    최종적인 아웃풋이 제품 별 평균 사용시간을 정의하는 것입니다.
    도움 부탁드립니다..ㅠ

  3. JHK 2019.12.02 01:40  댓글주소  수정/삭제  댓글쓰기

    이전의 질문에 대한 답변 감사했습니다!

    본문내용과는 좀 다르지만 전처리에 대한 질문을 하고싶어서요!

    데이터에서
    한id변수에 대해 기록된 정보가 여러개라서 예를 들어 id 변수에서 001이라는 아이디를 가진 사람의 정보가 약 200여개씩 존재합니다

    id
    001
    001
    001
    .
    .
    .
    이때 다른 것은 각 행마다 기록되는 시간이 다른데요, 그 기록된 시간마다 work라는 변수의 내용이 다릅니다.

    id work time
    001 작동준비 2018-10-23 03:22
    001 작동시작 2018-10-23 03:25
    001 ab투여 2018-10-23 03:43
    .
    .
    .

    이런식으로 구성되어 있는데요
    저는 work변수의 작동시작이라는 값이 해당하는 행에서 가지는 time변수의 시각을 기준으로 이전의 시각을 가지는 행들을 불러오고 싶습니다.

    어떻게 해야할지 잘 모르겠네요 ㅜㅜ
    time변수는 다루기 편하게 문자열과 공백을 지우고 대소비교가 가능하도록
    201810230322 이런식으로 바꾸어 놓은 상태입니다!

    • R Friend R_Friend 2019.12.02 10:23 신고  댓글주소  수정/삭제

      제가 저녁에 답글 남길께요

    • R Friend R_Friend 2019.12.02 23:44 신고  댓글주소  수정/삭제

      안녕하세요 JHK님,

      만약 id 별로 work 중에 '작동시작' 이 "단 한번만 존재"한다는 가정 하에

      (1) id별로 work = '작동시작' 인 행만을 가져와서 dataframe을 만들고 (df2)

      (2) 이를 id 키를 기준으로 원래의 dataframe(df)에 df2를 merge 한 후에

      (3) df2의 작동시작 시간(work_prepair_time) 보다 작은 time 인 경우만 가져와서 새로운 dataframe(df3)을 만드는 순서로 코드 작성해 보았습니다.

      (아래 코드에서 'w2' = '작동시작'으로 간주하시면 됩니다)

      > id <- c(rep('001', 3), rep('002', 4))
      > work <- c('w1', 'w2', 'w3', 'w0', 'w1', 'w2', 'w3')
      > time <- c(201810230322, 201810230325, 201810230343, 201810230222, 201810230322, 201810230325, 201810230343)
      >
      > df <- data.frame(id, work, time)
      > df
      id work time
      1 001 w1 201810230322
      2 001 w2 201810230325
      3 001 w3 201810230343
      4 002 w0 201810230222
      5 002 w1 201810230322
      6 002 w2 201810230325
      7 002 w3 201810230343
      >
      > work_prepair_time <- subset(df,
      + select = c('id', 'time'),
      + subset = (work == 'w2'))
      > names(work_prepair_time) <- c('id', 'work_prepair_time')
      > work_prepair_time
      id work_prepair_time
      2 001 201810230325
      6 002 201810230325
      >
      > df2 <- merge(df, work_prepair_time, by = 'id')
      > df2
      id work time work_prepair_time
      1 001 w1 201810230322 201810230325
      2 001 w2 201810230325 201810230325
      3 001 w3 201810230343 201810230325
      4 002 w0 201810230222 201810230325
      5 002 w1 201810230322 201810230325
      6 002 w2 201810230325 201810230325
      7 002 w3 201810230343 201810230325
      >
      > df3 <- subset(df2,
      + select = c('id', 'work', 'time'),
      + subset = (time < work_prepair_time))
      >
      > df3
      id work time
      1 001 w1 201810230322
      4 002 w0 201810230222
      5 002 w1 201810230322


      만약, 각 id 별로 work에서 '작동시작' 이 두번 이상 발생하는 경우도 있다면, 이럴 경우 처리하는 로직이 무엇이냐에 따라서 코드를 다르게 짜야 할텐데요, 혹시 이런 경우라면 댓글 남겨주세요.

    • JHK 2019.12.03 07:06  댓글주소  수정/삭제

      정말 친절한 설명 감사해요!ㅜㅜ!!
      막상 데이터분석을 해보니 전처리가 정말
      힘든 것 같아요! 감사합니다 :)

    • R Friend R_Friend 2019.12.03 09:28 신고  댓글주소  수정/삭제

      문제가 해결되었다니 기쁘네요.

이번 포스팅에서는 


(1) R ggplot2를 이용하여 커널 밀도 곡선(Kernel Density Curve)을 그리기

(2) 커밀 밀도 곡선의 최대 피크값의 좌표 구하기 
   (X, Y coordinates of the peak value in kernel density curve)

(3) 커널 밀도 곡선의 최대 피크값 위치에 수직선을 추가하는 방법 (adding a vertical line at peak value in kernel density plot) 을 소개하겠습니다. 




  (1) R ggplot2로 커널 밀도 곡선 그리기 (Kernel Density Curve by R ggplot2)


예제로 사용할 데이터는 MASS 패키지에 내장되어 있는 Cars93 데이터프레임입니다. 자동차의 가격(Price 변수)에 대해서 R ggplot2로 커널 밀도 곡선 (not frequency, but density) 을 그려보겠습니다. 


> library(MASS)

> library(ggplot2)

> ggplot(Cars93, aes(x=Price)) + 

+   geom_density(fill = 'yellow') + 

+   geom_line(stat = "density") + 

+   expand_limits(y = 0) + 

+   ggtitle("Kernel Density Curve") 






(2) 커널 밀도 곡선의 최대 피크값의 좌표 구하기 

     (X, Y coordinates of the peak value in kernel density curve)


density(df$var) 함수를 사용하면 x 와 y 값의 요약통계량 (최소, Q1, 중앙값, Q3, 최대값)을 확인할 수 있습니다. y 값의 최대값(Max)은 5.025e-02 값이네요. 


 

> density(Cars93$Price)


Call:

density.default(x = Cars93$Price)


Data: Cars93$Price (93 obs.); Bandwidth 'bw' = 3.011


       x                y            

 Min.   :-1.634   Min.   :1.608e-05  

 1st Qu.:16.508   1st Qu.:9.773e-04  

 Median :34.650   Median :5.329e-03  

 Mean   :34.650   Mean   :1.377e-02  

 3rd Qu.:52.792   3rd Qu.:2.078e-02  

 Max.   :70.934   Max.   :5.025e-02




이중에서 y 값 peak 값의 x 좌표를 알고 싶으므로 (a) which.max(density(Cars93$Price)$y) 로 y값의 최대 peak 값을 가지는 데이터의 index를 찾고 (127번째 데이터), 이 index를 이용해서 x 값의 좌표(16.25942)를 indexing 해올 수 있습니다. 



> which.max(density(Cars93$Price)$y)

[1] 127

> x_coord_of_y_max = density(Cars93$Price)$x[127]

> x_coord_of_y_max

[1] 16.25942

 


즉, y (Price) 최대 피크값의 좌표는 (x, y) 는 (16.25922, 0.05025) 가 되겠습니다. 




 (3) 커널 밀도 곡선의 최대 피크값 위치에 수직선을 추가하기

     (adding a vertical line at peak value in kernel density plot)


geom_vline() 으로 (2)번에서 구한 y 피크값의 x 좌표를 입력해주면 커널밀도곡선의 최대 피크값 위치에 수직선(vertical line)을 추가할 수 있습니다. 



> ggplot(Cars93, aes(x=Price)) + 

+   geom_density(fill = 'yellow') + 

+   geom_line(stat = "density") + 

+   expand_limits(y = 0) + 

+   ggtitle("Kernel Density Curve w/ Vertical Line at Peak Point") +

+   geom_vline(xintercept = x_coord_of_y_max, color='blue', size=2, linetype="dotted")


 





  (4) 여러개 그룹의 커널 밀도 곡선을 겹쳐 그린 경우 최대 peak 값 찾고 수직선 그리기


위의 예는 1개의 데이터셋에 대해 1개의 커널 밀도 곡선을 그린 후 최대 peak 값을 찾는 것이었습니다. 이제부터는 여러개의 하위 그룹으로 나뉘어진 데이터셋으로 여러개의 커널 밀도 곡선을 겹쳐서 그린 경우에 최대 peak 값을 찾고, 그 위치에 빨간색으로 수직선을 추가해보겠습니다. 


예제로 Cars93 데이터프레임에서 차종(Type) 별 가격(Price)의 커널 밀도 곡선을 겹쳐서 그려보겠습니다. 


> ggplot(Cars93, aes(x=Price, colour = Type)) + 

+   geom_density(fill = NA) + 

+   geom_line(stat = "density") + 

+   expand_limits(y = 0) + 

+   ggtitle("Kernel Density Curve by Car Types")





차종(Type) 중에서 Van 의 커널 밀도 곡선이 최대 Peak 값에 해당하므로 --> 차종 중에서 Van 의 데이터만 가져와서 y 최대값과 이에 해당하는 관측치의 index 위치를 찾아보겠습니다. 



> density(Cars93[Cars93$Type=='Van', ]$Price)


Call:

density.default(x = Cars93[Cars93$Type == "Van", ]$Price)


Data: Cars93[Cars93$Type == "Van", ]$Price (9 obs.); Bandwidth 'bw' = 0.303


       x               y            

 Min.   :15.39   Min.   :0.0000072  

 1st Qu.:17.45   1st Qu.:0.0040424  

 Median :19.50   Median :0.0495649  

 Mean   :19.50   Mean   :0.1215343  

 3rd Qu.:21.55   3rd Qu.:0.1807290  

 Max.   :23.61   Max.   :0.5321582  

> which.max(density(Cars93[Cars93$Type=='Van', 'Price'])$y)

[1] 238

> x_coord_of_y_max = density(Cars93[Cars93$Type=='Van','Price'])$x[238]

> x_coord_of_y_max

[1] 19.20249

 


밀도 y의 최대값(max)은 0.532 이며, 이때의 관측치의 index는 238번째 값 이고, 관측치 238번째 값의 x값의 좌표는 19.20 입니다. 


따라서, 밀도 y 의 최대 peak 값의 좌표는 (x, y) = (19.20, 0.53) 이네요. 



마지막으로, geom_vline() 을 사용해서 밀도 최대 peak 값에 해당하는 x=19.20 을 기준으로 빨간색 점선 수직선(vertical line)을 추가해보겠습니다. 



> ggplot(Cars93, aes(x=Price, colour = Type)) + 

+   geom_density(fill = NA) + 

+   geom_line(stat = "density") + 

+   expand_limits(y = 0) + 

+   ggtitle("Kernel Density Curve w/ Vertical Line at Peak Point") +

+   theme_bw() +

+   geom_vline(xintercept = x_coord_of_y_max, color='red', size=1, linetype="dashed")






R ggplot2로 히스토그램, 커널밀도곡선을 그리는 다양한 방법은 https://rfriend.tistory.com/67 를 참고하세요. 


많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. ^^


Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 스물스물 2019.10.04 14:36  댓글주소  수정/삭제  댓글쓰기

    올려주신 답변 보고 작성해주신 글 확인하였습니다. 감사합니다.

    추가적으로 질문을 드리면, 올려주신 글에서는 Cars93$Price에 대한 density plot을 그리고 그에 해당하는 x,y 좌표값을 구하는 방식 & 그래프에 최대 peak 를 나타내는 방법인것 같습니다.

    만약, density plot을 여려개를 겹쳐 그려서, 예를 들면, 내장되어 있는(Cars93$price)에는 num[1:93]까지의 데이터를 이루고 있지만, 만약 price1, pirce2, price3...이렇게 여러개의 데이터를 하나로 이루어서 density plot을 그렸다면, 그림자체도 여러개의 그래프가 겹쳐서 나타나게 될 것이고, 그에 따른 최대 peak값을 따로 구해야할것 같은데,
    개인적으로는, 전체 합쳐진 데이터내에서 최대 peak값을 구하고 싶은데 이럴 때에는 어떻게 해야 할까요?

이번 포스팅에서는 R의 DataFrame에서 특정 조건에 해당하는 값의 행과, 해당 행의 앞, 뒤 2개 행(above and below 2 rows) 을 동시에 제거하는 방법을 소개하겠습니다. 




예를 들어서, 아래와 같이 var1, var2의 두 개의 변수를 가지는 df라는 이름의 DataFrame이 있다고 했을 때, var2의 값 중 음수(-)인 값을 가지는 행과, 해당 행의 위, 아래 2개 행을 같이 제거(remove, filter) 해서 df2 라는 이름의 새로운 DataFrame을 만들어보겠습니다. 



> var1 <- c(1:12)

> var2 <- c(100, -200, 101, 1102, 50, 300, 100, 400, -100, 82, 90, 80)

> df <- data.frame(var1, var2)

> df

   var1 var2

1     1  100

2     2 -200

3     3  101

4     4 1102

5     5   50

6     6  300

7     7  100

8     8  400

9     9 -100

10   10   82

11   11   90

12   12   80

 



먼저, 칼럼 var2 에서 음수(-)인 값의 조건을 만족하는 값의 행의 위치를 indexing 한 negative_idx 벡터를 만들어보겠습니다. 



> # condition: if var2 value is negative, then TRUE

> negative_idx <- df$var2 < 0

> negative_idx

 [1] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE




다음으로 칼럼 var2의 값이 음수(-)인 값의 위치를 기준으로 해당 값의 행과 앞, 뒤 2개행까지는 제거(remove, filter)하고, 그 외의 값은 유지(keep_idx = TRUE) 하는 keep_idx 라는 벡터를 for loop 반복문과 if 조건문을 사용해서 만들어보겠습니다. 



> keep_idx <- c(rep(TRUE, length(df$var2)))

> keep_idx

 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

> for (i in 1:length(negative_idx)) {

 if (negative_idx[i] == TRUE) {

+     keep_idx[i-2] = FALSE

+     keep_idx[i-1] = FALSE

+     keep_idx[i] = FALSE

+     keep_idx[i+1] = FALSE

+     keep_idx[i+2] = FALSE

 }

+ }

> keep_idx

 [1] FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE

 



이제 df라는 DataFrame에서 위에서 구한 keep_idx = TRUE 인 행만 indexing 해서 (즉, keep_idx = FALSE 인 행은 제거) 새로운 df2 라는 DataFrame을 만들어보겠습니다. 



> # subset only rows with keep_idx = TRUE

> df_filstered <- df[keep_idx, ]

> df_filstered

   var1 var2

5     5   50

6     6  300

12   12   80




위는 예는 조건문과 반복문을 사용해서 indexing 해오는 방법이었구요, windows 함수인 lag(), lead()를 사용해서도 동일한 기능을 수행하는 프로그램을 짤 수도 있습니다. 다만, 이번 포스팅에서 소개한 코드가 좀더 범용적이고 코드도 짧기 때문에 lag(), lead() 함수를 사용한 방법은 추가 설명은 하지 않겠습니다. 


많은 도움이 되었기를 바랍니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요.



Posted by R Friend R_Friend

댓글을 달아 주세요

이번 포스팅에서는 R Shiny를 사용하여

 

- (1) 두 개의 독립 표본 집단 간 평균 차이 t-검정 (independent two-sample t-test) 

- (2) 신뢰수준 별 신뢰구간 (confidence interval by confidence level)

을 구하는 웹 애플리케이션을 만들어보겠습니다. 

 

예제로 사용한 데이터는 MASS 패키지의 Cars93 데이터프레임이며, 'Origin'  변수의 두 개 집단(USA, Non-USA) 별로 Price, RPM, Length, Width 변수 간의 t-test 를 해보겠습니다. 

 

왼쪽 사이드바 패널에는 t-검정을 할 변수를 선택할 수 있는 콤보 박스와, 신뢰수준(confidence level)을 지정할 수 있는 슬라이드 입력(default = 0.95) 화면을 만들었습니다. 

그리고 메인 패널에는 변수 이름, 신뢰수준, t-검정 결과를 텍스트로 보여줍니다. 

 

R Shiny - independent two-sample t-test

 

library(shiny)

# Define UI for application that does t-test

ui <- fluidPage(

   

   # Application title

   titlePanel("Independent two sample t-test"),

   

   # Select a variable to do t-test by Origin (USA vs. non-USA groups)

   sidebarPanel(

     selectInput("var_x", 

                 label = "Select a Variable to do t-test", 

                 choices = list("Price"="Price", 

                                "RPM"="RPM", 

                                "Length"="Length", 

                                "Width"="Width"), 

                 selected = "Price"), 

     

     sliderInput("conf_level", "Confidence Level:", 

                  min = 0.05, max = 0.99, value = 0.95)

     ),

   

   # Show a t-test result

   mainPanel(

     h5("Variable is:"),

     verbatimTextOutput("var_x_name"), 

     h5("Confidence Level is:"),

     verbatimTextOutput("conf_level"), 

     h4("t-test Result by Origin (USA vs. non-USA independent two samples)"), 

     verbatimTextOutput("t_test")

   )

)

# Define server logic required to do t-test

server <- function(input, output) {

  

  library(MASS)

  cars93_sub <- Cars93[, c('Origin', 'Price', 'RPM', 'Length', 'Width')]

  

  output$var_x_name <- renderText({

    as.character(input$var_x)

  })

  

  output$conf_level <- renderText({

    as.character(input$conf_level)

  })

  

  output$t_test <- renderPrint({

    # independent two-sample t-test

    t.test(cars93_sub[,input$var_x] ~ Origin, 

           data = cars93_sub, 

           alternative = c("two.sided"), 

           var.equal = FALSE, 

           conf.level = input$conf_level)

  })

}

 

# Run the application 

shinyApp(ui = ui, server = server)

 

많은 도움이 되었기를 바랍니다. 

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 꾸리꾸리 2019.09.10 18:17  댓글주소  수정/삭제  댓글쓰기

    최근 shiny 를 공부하고 있는데 많은 도움이 되었습니다.
    통계 관련해서는 자세히 모르지만, 통계분석도 시각화가 가능한가요?

이번 포스팅에서는 R Shiny를 사용하여 interactive하게 편하고 쉽게 두 연속형 변수 간 관계를 알아볼 수 있도록 

(1) 상관계수 (Correlation Coefficient)를 구하고, 

(2) 산점도 (Scatter Plot)을 그린 후, 

(3) 선형 회귀선 (Linear Regression Line) 을 겹쳐서 그려주는 

웹 애플리케이션을 만들어보겠습니다. 

 

 

R Shiny - Correlation Coefficient, Scatter Plot, Linear Regression Line

 

예제로 사용할 데이터는 MASS 패키지에 내장되어 있는 Cars93 데이터프레임으로서, Price, MPG.city, MPG.highway, EngineSize, Horsepower, RPM 의 6개 연속형 변수를 선별하였습니다. 

 

왼쪽 사이드바 패널에는 산점도의 X축, Y축에 해당하는 연속형 변수를 콤보박스로 선택할 수 있도록 하였으며, 상관계수와 선형회귀모형을 적합할 지 여부를 선택할 수 있는 체크박스도 추가하였습니다. 

 

오른쪽 본문의 메인 패널에는 두 연속형 변수에 대한 산점도 그래프에 선형 회귀선을 겹쳐서 그려서 보여주고, 상관계수와 선형회귀 적합 결과를 텍스트로 볼 수 있도록 화면을 구성하였습니다. 

 

library(shiny)

# Define UI for application that analyze correlation coefficients and regression

ui <- fluidPage(

   

   # Application title

   titlePanel("Correlation Coefficients and Regression"),

   

   # Select 2 Variables, X and Y 

   sidebarPanel(

     selectInput("var_x", 

                 label = "Select a Variable of X", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "RPM"),

     

     selectInput("var_y", 

                 label = "Select a Variable of Y", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "MPG.highway"), 

     

     checkboxInput(inputId = "corr_checked", 

                   label = strong("Correlation Coefficients"), 

                   value = TRUE), 

     

     checkboxInput(inputId = "reg_checked", 

                   label = strong("Simple Regression"), 

                   value = TRUE)

      ),

      # Show a plot of the generated distribution

      mainPanel(

        h4("Scatter Plot"), 

        plotOutput("scatterPlot"), 

        

        h4("Correlation Coefficent"), 

        verbatimTextOutput("corr_coef"), 

        

        h4("Simple Regression"), 

        verbatimTextOutput("reg_fit")

      )

   )

# Define server logic required to analyze correlation coefficients and regression

server <- function(input, output) {

  library(MASS)

  scatter <- Cars93[,c("Price", "MPG.city", "MPG.highway", "EngineSize", "Horsepower", "RPM")]

  

  # scatter plot

  output$scatterPlot <- renderPlot({

    var_name_x <- as.character(input$var_x)

    var_name_y <- as.character(input$var_y)

    

    plot(scatter[, input$var_x],

         scatter[, input$var_y],

         xlab = var_name_x,

         ylab = var_name_y,

         main = "Scatter Plot")

    

    # add linear regression line

    fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

    abline(fit)

    })

  

  # correlation coefficient

  output$corr_coef <- renderText({

    if(input$corr_checked){

      cor(scatter[, input$var_x], 

          scatter[, input$var_y])

      }

    })

  

  # simple regression

  output$reg_fit <- renderPrint({

    if(input$reg_checked){

      fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

      names(fit$coefficients) <- c("Intercept", input$var_x)

      summary(fit)$coefficients

    }

  })

}

# Run the application 

shinyApp(ui = ui, server = server)

 

많은 도움이 되었기를 바랍니다. 

 

Posted by R Friend R_Friend

댓글을 달아 주세요

  1. 2019.12.05 01:44  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다