Package loading and data preparation
library(ggplot2)
library(dplyr)
library(chngpt)
data(diamonds,package="ggplot2")#use diamond data from ggplot2 package
set.seed(888)#set seed for random number generator
dia=sample_n(diamonds, 100)#sample 100 data
names(dia)
## [1] "carat" "cut" "color" "clarity" "depth" "table" "price"
## [8] "x" "y" "z"
Model 1: Step model
fit1=chngptm(formula.1=price~color, formula.2=~carat, family="gaussian", data=dia,type="step", var.type="bootstrap", weights=NULL,test.statistic="lr")
summary(fit1)
## Change point model threshold.type: step
##
## Coefficients:
## est Std. Error* (lower upper) p.value*
## (Intercept) 2388.78139 357.8547 1614.7555 3017.546 2.467557e-11
## color.L -117.72317 1262.3233 -1766.5451 3181.762 9.256977e-01
## color.Q 22.17552 1327.0908 -2086.7443 3115.452 9.866681e-01
## color.C -203.13155 1288.2854 -1964.8426 3085.236 8.747122e-01
## color^4 542.00841 680.0187 -561.0645 2104.609 4.254224e-01
## color^5 876.32034 609.5447 -486.3500 1903.065 1.505295e-01
## color^6 450.28602 536.3013 -320.4673 1781.834 4.011249e-01
## carat>chngpt 7796.50278 1594.0398 4374.4380 10623.074 1.003077e-06
##
## Threshold:
## est Std. Error (lower upper)
## 1.170000 0.247449 0.320000 1.290000
plot(fit1)
Model 2: Hinge model
fit2=chngptm(formula.1=price~color, formula.2=~carat, family="gaussian", data=dia,type="hinge", var.type="bootstrap", weights=NULL,test.statistic="lr")
summary(fit2)
## Change point model threshold.type: hinge
##
## Coefficients:
## est Std. Error* (lower upper) p.value*
## (Intercept) 775.1345 181.4693 367.7466 1079.1063 1.942189e-05
## color.L -1049.0613 526.4154 -2257.5414 -193.9929 4.627906e-02
## color.Q -539.4041 477.8035 -1573.3707 299.6190 2.589296e-01
## color.C -306.0323 404.5805 -1153.5712 432.3844 4.493981e-01
## color^4 -160.4980 352.6540 -897.5327 484.8711 6.490268e-01
## color^5 -272.2834 247.5320 -765.4933 204.8320 2.713352e-01
## color^6 571.2956 229.2629 103.6430 1002.3537 1.270691e-02
## (carat-chngpt)+ 9932.2172 957.1257 7826.3223 11578.2552 3.151096e-25
##
## Threshold:
## est Std. Error (lower upper)
## 0.5600000 0.2040816 0.4000000 1.2000000
plot(fit2)
Model 3: Segmented model
fit3=chngptm(formula.1=price~color, formula.2=~carat, family="gaussian", data=dia,type="segmented", var.type="bootstrap", weights=NULL)
summary(fit3)
## Change point model threshold.type: segmented
##
## Coefficients:
## est Std. Error* (lower upper) p.value*
## (Intercept) -1384.1387 352.4071 -2113.23556 -731.7996 8.577298e-05
## color.L -1240.0324 819.9517 -3330.77934 -116.5688 1.304515e-01
## color.Q -748.8264 760.0964 -2510.21942 469.3586 3.245392e-01
## color.C -526.1910 567.5413 -1826.81567 397.9464 3.538532e-01
## color^4 -217.6639 448.6655 -1288.93820 469.8304 6.275797e-01
## color^5 -377.7491 245.7066 -852.57939 110.5904 1.241955e-01
## color^6 435.4351 232.2382 -40.21034 870.1635 6.079959e-02
## carat 5589.7985 585.7028 4226.79977 6522.7548 1.377650e-21
## (carat-chngpt)+ 6161.6928 1463.0538 3787.65213 9522.8228 2.536484e-05
##
## Threshold:
## est Std. Error (lower upper)
## 0.9600000 0.2168367 0.3200000 1.1700000
plot(fit3)
Model 4: Three-segmented model
fit4=chngptm(formula.1=price~color, formula.2=~carat, family="gaussian", data=dia,type="M111", var.type="bootstrap", weights=NULL)
summary(fit4)
## Change point model threshold.type: M111
##
## Coefficients:
## est Std. Error* (lower upper)
## (Intercept) -6072.8588 5164.9726 -12762.06927 7484.62337
## color.L -1201.8119 783.8722 -3021.33378 51.44517
## color.Q -692.8943 759.4676 -2393.47326 583.63989
## color.C -539.1010 576.5238 -1791.77970 468.19352
## color^4 -225.4975 469.2834 -1227.16060 612.43032
## color^5 -305.2114 274.6834 -866.11058 210.64819
## color^6 410.2322 239.0187 -72.99389 863.95950
## carat 11024.4067 3188.5911 2646.31688 15145.59392
## (carat-chngpt1)- -50841.8269 55747.1575 -107801.80661 110727.05085
## (carat-chngpt2)- 45805.2499 55269.7962 -113126.63231 103530.96876
## p.value*
## (Intercept) 0.2396837811
## color.L 0.1252331692
## color.Q 0.3615886595
## color.C 0.3497424613
## color^4 0.6308615711
## color^5 0.2665085842
## color^6 0.0861038272
## carat 0.0005453055
## (carat-chngpt1)- 0.3617647584
## (carat-chngpt2)- 0.4072417258
##
## Threshold:
## est Std. Error (lower upper)
## chngpt.1 1.17 0.2448980 0.31 1.27
## chngpt.2 1.20 0.2015306 0.71 1.50
plot(fit4)
cat("Smaller AIC value is better","\n","Step ","Hinge "," Segmented","3-Segmented","\n",AIC(fit1),AIC(fit2),AIC(fit3),AIC(fit4))
## Smaller AIC value is better
## Step Hinge Segmented 3-Segmented
## 1824.967 1688.038 1668.149 1667.915