ANALISIS DATASET 2
Erik Fariz N
2026-10-06
#MEMANGGIL LIBRARY YANG DIGUNAKAN#
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(MASS)
library(car)## Loading required package: carDatalibrary(pastecs)
library(pracma)##
## Attaching package: 'pracma'## The following object is masked from 'package:car':
##
## logitlibrary(Matrix)##
## Attaching package: 'Matrix'## The following objects are masked from 'package:pracma':
##
## expm, lu, tril, triu#MEMANGGIL DATA#
Catatan: pada data, kolom Y berada di kolom terakhir (X1 … X6, Y). Fungsi GCV dan uji pada sintaks ini mengambil kolom pertama (data[,1]) sebagai respon, sehingga kolom data disusun ulang agar Y berada di kolom pertama. Jika tidak disusun ulang, X1 akan ikut berperan sebagai respon.
folder="C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\TUGAS 6" #folder untuk menyimpan/membaca file csv hasil
data=read.table(file.choose(),header=TRUE)
data=data[,c("Y","X1","X2","X3","X4","X5","X6")] #Y dijadikan kolom pertama
datanrow(data)## [1] 514#ANALISIS STATISTIKA DESKRIPTIF & DETEKSI MULTIKOLINIERITAS#
stat.desc(data)#DETEKSI PENCILAN (ATURAN IQR)#
hitung_pencilan=function(v)
{
q=quantile(v,c(.25,.75))
iqr=diff(q)
sum(v<q[1]-1.5*iqr | v>q[2]+1.5*iqr)
}
pencilan=data.frame(Variabel=names(data),
Banyak_pencilan=sapply(data,hitung_pencilan),
Persen=round(100*sapply(data,hitung_pencilan)/nrow(data),2),
row.names=NULL)
pencilan#KORELASI ANTAR VARIABEL#
kor=cor(data)
round(kor,3)## Y X1 X2 X3 X4 X5 X6
## Y 1.000 -0.681 -0.667 0.457 -0.592 0.629 -0.112
## X1 -0.681 1.000 0.534 -0.518 0.437 -0.563 0.191
## X2 -0.667 0.534 1.000 -0.541 0.425 -0.523 0.118
## X3 0.457 -0.518 -0.541 1.000 -0.335 0.354 -0.019
## X4 -0.592 0.437 0.425 -0.335 1.000 -0.389 -0.080
## X5 0.629 -0.563 -0.523 0.354 -0.389 1.000 -0.273
## X6 -0.112 0.191 0.118 -0.019 -0.080 -0.273 1.000image(1:ncol(kor),1:ncol(kor),t(kor[ncol(kor):1,]),
col=colorRampPalette(c("#b2182b","white","#2166ac"))(50),
zlim=c(-1,1),axes=FALSE,xlab="",ylab="",main="Matriks Korelasi Pearson")
axis(1,at=1:ncol(kor),labels=colnames(kor))
axis(2,at=1:ncol(kor),labels=rev(colnames(kor)),las=1)
for (i in 1:ncol(kor))
{
for (j in 1:ncol(kor))
{
text(j,ncol(kor)-i+1,sprintf("%.2f",kor[i,j]),cex=0.8)
}
}
#SCATTERPLOT#
Garis merah = garis linear, garis biru putus-putus = kurva lowess. Jika kurva lowess menyimpang dari garis linear (melengkung atau patah), hal itu mendukung penggunaan regresi nonparametrik spline.
plot(data$X1,data$Y,main="Scatterplot Y dan X1"
, xlab="X1"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X1), col="darkred", lwd=3)
lines(lowess(data$X1,data$Y), col="steelblue", lwd=3, lty=2)
plot(data$X2,data$Y,main="Scatterplot Y dan X2"
, xlab="X2"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X2), col="darkred", lwd=3)
lines(lowess(data$X2,data$Y), col="steelblue", lwd=3, lty=2)
plot(data$X3,data$Y,main="Scatterplot Y dan X3"
, xlab="X3"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X3), col="darkred", lwd=3)
lines(lowess(data$X3,data$Y), col="steelblue", lwd=3, lty=2)
plot(data$X4,data$Y,main="Scatterplot Y dan X4"
, xlab="X4"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X4), col="darkred", lwd=3)
lines(lowess(data$X4,data$Y), col="steelblue", lwd=3, lty=2)
plot(data$X5,data$Y,main="Scatterplot Y dan X5"
, xlab="X5"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X5), col="darkred", lwd=3)
lines(lowess(data$X5,data$Y), col="steelblue", lwd=3, lty=2)
plot(data$X6,data$Y,main="Scatterplot Y dan X6"
, xlab="X6"
, ylab="Y"
,col="black")
abline(lm(data$Y~data$X6), col="darkred", lwd=3)
lines(lowess(data$X6,data$Y), col="steelblue", lwd=3, lty=2)
#HISTOGRAM DAN BOXPLOT TIAP PREDIKTOR#
for (j in 2:7)
{
par(mfrow=c(1,2))
hist(data[,j],breaks=30,col="steelblue",border="white",
main=paste("Histogram",names(data)[j]),xlab=names(data)[j],ylab="Frekuensi")
boxplot(data[,j],horizontal=TRUE,col="#fdae61",outcol="red",
main=paste("Boxplot",names(data)[j]),xlab=names(data)[j])
}
par(mfrow=c(1,1))#RINGKASAN SEMUA SCATTERPLOT DALAM SATU GAMBAR#
par(mfrow=c(2,3))
for (j in 2:7)
{
plot(data[,j],data$Y,main=paste("Y dan",names(data)[j]),
xlab=names(data)[j],ylab="Y",col="grey30",pch=16,cex=0.6)
abline(lm(data$Y~data[,j]), col="darkred", lwd=2)
lines(lowess(data[,j],data$Y), col="steelblue", lwd=2, lty=2)
}
par(mfrow=c(1,1))#MODEL REGRESI LINEAR BERGANDA (PEMBANDING)#
model=(lm(formula=Y~X1+X2+X3+X4+X5+X6,data=data))
model##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
##
## Coefficients:
## (Intercept) X1 X2 X3 X4 X5
## 26.47106 -0.61512 -0.21398 -0.21664 -0.42017 0.92267
## X6
## 0.03101summary(model)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.041 -3.610 2.047 5.431 21.434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.47106 11.16110 2.372 0.0181 *
## X1 -0.61512 0.07038 -8.740 < 2e-16 ***
## X2 -0.21398 0.02453 -8.721 < 2e-16 ***
## X3 -0.21664 0.29665 -0.730 0.4656
## X4 -0.42017 0.04933 -8.518 < 2e-16 ***
## X5 0.92267 0.14182 6.506 1.85e-10 ***
## X6 0.03101 0.04248 0.730 0.4657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.209 on 507 degrees of freedom
## Multiple R-squared: 0.6801, Adjusted R-squared: 0.6763
## F-statistic: 179.7 on 6 and 507 DF, p-value: < 2.2e-16vif(model)## X1 X2 X3 X4 X5 X6
## 1.930060 1.842417 1.599757 1.415399 1.754443 1.155668par(mfrow=c(2,2))
plot(model)
par(mfrow=c(1,1))shapiro.test(residuals(model))##
## Shapiro-Wilk normality test
##
## data: residuals(model)
## W = 0.91328, p-value < 2.2e-16bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 47.042, df = 6, p-value = 1.835e-08dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 1.4363, p-value = 3.087e-11
## alternative hypothesis: true autocorrelation is greater than 0#PEMILIHAN 1 TITIK KNOT#
GCV1=function(data)
{
para=0
data=as.matrix(data)
N=nrow(data)
M=ncol(data)
m=ncol(data)-para-1
dataA=data[,(para+2):M]
dataA=as.matrix(dataA)
F=diag(N)
nk=50 #nk=banyaknya alternatif titik knot yang akan dicoba
knot1=matrix(ncol=m,nrow=nk)
for (i in (1:m)) #membuat knot
{
a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk)
knot1[,i]=t(as.matrix(a))
}
a1=length(knot1[,1])
knot1=as.matrix(knot1[2:(a1-1),])
colnames(knot1)=paste0("X",1:m)
aa=rep(1,N)
data1=matrix(ncol=m,nrow=N)
data2=data[,2:M] #data x saja
nk1=nrow(knot1)
GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV"
MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE"
SSE=rep(NA,nk1)
SSR=rep(NA,nk1)
Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq"
knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke"
for (i in 1:nk1)
{
for (j in 1:m)
{
for (k in 1:N)
{
if (data[k,(j+para+1)]<knot1[i,j]) data1[k,j]=0 else
data1[k,j]=data[k,(j+para+1)]-knot1[i,j]
}
}
mx=as.matrix(cbind(aa,data2,data1))
C=pinv(t(mx)%*%mx)
B=C%*%(t(mx)%*%data[,1])
yhat=mx%*%B
res=data[,1]-yhat
SSE[i]=sum((res)^2)
SSR[i]=sum((yhat-mean(data[,1]))^2)
MSE[i]=SSE[i]/(N)
Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100
A=mx%*%C%*%t(mx)
A1=(F-A)
A2=(sum(diag(A1))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot1))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll, file=paste0(folder,"DATA ALL knot 1.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 1 knot","\n")
cat("==============================================","\n")
print(((dataG[1,1:ncol(dataG)])))
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
# --- hitung ulang B dari titik knot dengan GCV TERKECIL ---
baris_terbaik = dataG[1, "knot_ke"]
knotgcv = knot1[baris_terbaik, ]
datagcv1=matrix(ncol=m,nrow=N)
for (j in 1:m)
{
for (k in 1:N)
{
if (data[k,(j+para+1)]<knotgcv[j]) datagcv1[k,j]=0 else
datagcv1[k,j]=data[k,(j+para+1)]-knotgcv[j]
}
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
Cgcv=pinv(t(mxgcv)%*%mxgcv)
B=Cgcv%*%(t(mxgcv)%*%data[,1])
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 1","\n")
cat("==============================================","\n")
print(B)
cat("\n")
}
GCV1(data)## ==============================================
## HASIL GCV terkecil dengan 1 knot
## ==============================================
## GCV knot_ke X1 X2 X3 X4 X5 X6
## 60.891404 15.000000 13.823061 30.495918 4.788163 26.577551 62.518980 15.367347
## Nilai GCV 10 terkecil pertama
## GCV knot_ke X1 X2 X3 X4 X5 X6
## [1,] 60.89140 15 13.82306 30.49592 4.788163 26.57755 62.51898 15.36735
## [2,] 60.90658 16 14.59327 32.52898 5.022041 28.34939 62.97224 16.39184
## [3,] 60.93708 17 15.36347 34.56204 5.255918 30.12122 63.42551 17.41633
## [4,] 61.00914 18 16.13367 36.59510 5.489796 31.89306 63.87878 18.44082
## [5,] 61.03068 19 16.90388 38.62816 5.723673 33.66490 64.33204 19.46531
## [6,] 61.03475 14 13.05286 28.46286 4.554286 24.80571 62.06571 14.34286
## [7,] 61.08612 20 17.67408 40.66122 5.957551 35.43673 64.78531 20.48980
## [8,] 61.22630 21 18.44429 42.69429 6.191429 37.20857 65.23857 21.51429
## [9,] 61.33258 22 19.21449 44.72735 6.425306 38.98041 65.69184 22.53878
## [10,] 61.42555 13 12.28265 26.42980 4.320408 23.03388 61.61245 13.31837
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 1
## ==============================================
## [,1]
## [1,] 55.235597472
## [2,] -0.060206537
## [3,] -0.139780329
## [4,] 4.602967604
## [5,] -0.613837239
## [6,] 0.008921616
## [7,] 0.059789517
## [8,] -0.866308131
## [9,] -0.081699055
## [10,] -5.057603077
## [11,] 0.720183134
## [12,] 1.020919183
## [13,] -0.129832245#UJI SIGNIFIKANSI SATU TITIK KNOT#
Fungsi uji di bawah dipakai untuk 1, 2, dan 3 titik knot (argumen K = banyak titik knot). Derajat bebas galat yang dipakai adalah standar (MSE = SSE/(n-p)), sehingga konsisten dengan summary(lm).
uji=function(alpha,para,K)
{
# --- baca & siapkan titik knot terbaik (GCV terkecil) ---
knot = read.csv(paste0(folder,"DATA ALL knot ",K,".csv"), header=TRUE)
knot = as.matrix(knot)
knot = knot[,-1] # buang kolom index baris
knot = knot[order(knot[,1]), ] # urutkan berdasarkan GCV (kolom 1) terkecil
baris_knot = as.numeric(knot[1, 4:ncol(knot)]) # knot terbaik X1-X6
data=as.matrix(data)
ybar=mean(data[,1])
m=para+2
n=nrow(data)
q=ncol(data)
p=q-para-1 # jumlah variabel X (X1-X6)
dataA=as.matrix(data[,rep(m:(m+p-1),each=K)]) # tiap variabel X diulang sebanyak K knot
satu=rep(1,n)
n1=length(baris_knot) # jumlah kolom knot = p*K
data.knot=matrix(ncol=n1,nrow=n)
# --- hitung fungsi basis knot (loop i untuk knot, j untuk observasi) ---
for (i in 1:n1)
{
for (j in 1:n)
{
if(dataA[j,i] < baris_knot[i]) data.knot[j,i]=0
else data.knot[j,i] = dataA[j,i] - baris_knot[i]
}
}
mx=satu
nama="(Intercept)"
for (j in 1:p)
{
mx=cbind(mx,data[,j+1],data.knot[,((j-1)*K+1):(j*K)])
nama=c(nama,paste0("X",j),paste0("X",j,"_k",1:K))
}
mx=as.matrix(mx)
B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1]
rownames(B)=nama
n1=nrow(B)
yhat=mx%*%B
ybar=mean(data[,1])
res=data[,1]-yhat
SSE=sum((data[,1]-yhat)^2)
SSR=sum((yhat-ybar)^2)
MSE=SSE/(n-n1)
MSR=SSR/(n1-1)
SST=sum((data[,1]-ybar)^2)
Rsq=(SSR/(SSR+SSE))*100
Rsqadj=(1-(SSE/(n-n1))/(SST/(n-1)))*100
#-------------------------------------------------------#
# UJI SIMULTAN
#-------------------------------------------------------#
Fhit=MSR/MSE
pvalue=pf(Fhit,(n1-1),(n-n1),lower.tail=FALSE)
if(pvalue<=alpha)
{
cat('---------------------------------------','\n')
cat('Kesimpulan hasil uji simultan','\n')
cat('---------------------------------------','\n')
cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan','\n')
cat('','\n')
}
else
{
cat('---------------------------------------','\n')
cat('Kesimpulan hasil uji simultan','\n')
cat('---------------------------------------','\n')
cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan','\n')
cat('','\n')
}
#------------------------------------------------------#
# UJI PARSIAL
#------------------------------------------------------#
thit=rep(NA,n1)
pval=rep(NA,n1)
SE=sqrt(diag(MSE*(pinv(t(mx)%*%mx))))
cat('---------------------------------------------','\n')
cat('Kesimpulan hasil uji parsial','\n')
cat('---------------------------------------------','\n')
for (i in 1:n1)
{
thit[i]=B[i,1]/SE[i]
pval[i]=2*(pt(abs(thit[i]),(n-n1),lower.tail=FALSE))
if (pval[i]<=alpha) cat('Tolak Ho yakni parameter',nama[i],'signifikan dengan pvalue',pval[i],'\n')
else cat('Gagal tolak Ho yakni parameter',nama[i],'tidak signifikan dengan pvalue',pval[i],'\n')
}
thit=as.matrix(thit)
cat('=============================================','\n')
cat('nilai t hitung','\n')
cat('=============================================','\n')
rownames(thit)=nama
print(thit)
cat('=============================================','\n')
cat('Tabel uji parsial','\n')
cat('=============================================','\n')
tabel=data.frame(Estimasi=B[,1],Galat_baku=SE,t_hitung=thit[,1],p_value=pval,
Keputusan=ifelse(pval<=alpha,"Signifikan","Tidak signifikan"))
print(tabel,digits=4)
cat('Parameter signifikan:',sum(pval<=alpha),'dari',n1,'\n')
cat('Analysis of Variance','\n')
cat('=============================================','\n')
cat('Sumber df SS MS Fhit','\n')
cat('Regresi ',(n1-1),' ',SSR,' ',MSR,' ',Fhit,'\n')
cat('Error ',n-n1,' ',SSE,' ',MSE,'\n')
cat('Total ',n-1,' ',SST,'\n')
cat('=============================================','\n')
cat('s=',sqrt(MSE),' Rsq=',Rsq,' Rsq adj=',Rsqadj,'\n')
cat('pvalue(F)=',pvalue,'\n')
write.csv(res, file=paste0(folder,'OUTPUT UJI RESIDUAL knot',K,'.csv'))
write.csv(mx, file=paste0(folder,'OUTPUT UJI MX knot',K,'.csv'))
write.csv(yhat, file=paste0(folder,'OUTPUT UJI yhat knot',K,'.csv'))
}
uji(0.05, 0, 1)## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Gagal tolak Ho yakni parameter (Intercept) tidak signifikan dengan pvalue 0.3883954
## Gagal tolak Ho yakni parameter X1 tidak signifikan dengan pvalue 0.6170333
## Tolak Ho yakni parameter X1_k1 signifikan dengan pvalue 6.567427e-06
## Tolak Ho yakni parameter X2 signifikan dengan pvalue 0.003580801
## Gagal tolak Ho yakni parameter X2_k1 tidak signifikan dengan pvalue 0.223957
## Tolak Ho yakni parameter X3 signifikan dengan pvalue 0.005646631
## Tolak Ho yakni parameter X3_k1 signifikan dengan pvalue 0.003879527
## Tolak Ho yakni parameter X4 signifikan dengan pvalue 5.156252e-13
## Tolak Ho yakni parameter X4_k1 signifikan dengan pvalue 1.687891e-05
## Gagal tolak Ho yakni parameter X5 tidak signifikan dengan pvalue 0.9927584
## Gagal tolak Ho yakni parameter X5_k1 tidak signifikan dengan pvalue 0.3135006
## Gagal tolak Ho yakni parameter X6 tidak signifikan dengan pvalue 0.6565236
## Gagal tolak Ho yakni parameter X6_k1 tidak signifikan dengan pvalue 0.4467439
## =============================================
## nilai t hitung
## =============================================
## [,1]
## (Intercept) 0.863282194
## X1 -0.500371464
## X1_k1 -4.555637852
## X2 -2.926755428
## X2_k1 -1.217579486
## X3 2.779633244
## X3_k1 -2.901311348
## X4 -7.417012199
## X4_k1 4.344876881
## X5 0.009080621
## X5_k1 1.008918804
## X6 0.444983846
## X6_k1 -0.761454471
## =============================================
## Tabel uji parsial
## =============================================
## Estimasi Galat_baku t_hitung p_value Keputusan
## (Intercept) 55.235598 63.98325 0.863282 3.884e-01 Tidak signifikan
## X1 -0.060207 0.12032 -0.500371 6.170e-01 Tidak signifikan
## X1_k1 -0.866308 0.19016 -4.555638 6.567e-06 Signifikan
## X2 -0.139780 0.04776 -2.926755 3.581e-03 Signifikan
## X2_k1 -0.081699 0.06710 -1.217579 2.240e-01 Tidak signifikan
## X3 4.602968 1.65596 2.779633 5.647e-03 Signifikan
## X3_k1 -5.057603 1.74321 -2.901311 3.880e-03 Signifikan
## X4 -0.613837 0.08276 -7.417012 5.156e-13 Signifikan
## X4_k1 0.720183 0.16575 4.344877 1.688e-05 Signifikan
## X5 0.008922 0.98249 0.009081 9.928e-01 Tidak signifikan
## X5_k1 1.020919 1.01189 1.008919 3.135e-01 Tidak signifikan
## X6 0.059790 0.13436 0.444984 6.565e-01 Tidak signifikan
## X6_k1 -0.129832 0.17051 -0.761454 4.467e-01 Tidak signifikan
## Parameter signifikan: 6 dari 13
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 12 77085.11 6423.759 108.2327
## Error 501 29735.03 59.35135
## Total 513 106820.1
## =============================================
## s= 7.703983 Rsq= 72.16346 Rsq adj= 71.49672
## pvalue(F)= 1.293317e-130#PEMILIHAN DUA TITIK KNOT#
GCV2=function(data)
{
para=0
data=as.matrix(data)
N=nrow(data)
M=ncol(data)
m=ncol(data)-para-1
dataA=data[,(para+2):M]
dataA=as.matrix(dataA)
F=diag(N)
nk=50
knot1=matrix(ncol=m,nrow=nk)
for (i in (1:m))
{
a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk)
knot1[,i]=t(as.matrix(a))
}
a1=length(knot1[,1])
knot1=as.matrix(knot1[2:(a1-1),])
a2=nk-2
z=(a2*(a2-1)/2)
knot2=cbind(rep(NA,(z+1)))
for (i in (1:m))
{
knot=rbind(rep(NA,2))
for ( j in 1:(a2-1))
{
for (k in (j+1):a2)
{ xx=cbind(knot1[j,i],knot1[k,i])
knot=rbind(knot,xx)
}
}
knot2=cbind(knot2,knot)
}
knot2=knot2[2:(z+1),2:(2*m+1)]
colnames(knot2)=paste0(rep(paste0("X",1:m),each=2),"_k",rep(1:2,m))
a3=nrow(knot2)
aa=rep(1,N)
data1=matrix(ncol=2*m,nrow=N)
data2=data[,2:M]
nk1=nrow(knot2)
GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV"
MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE"
SSE=rep(NA,nk1)
SSR=rep(NA,nk1)
Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq"
knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke"
for (i in 1:a3)
{
for (j in 1:(2*m))
{
if (mod(j,2)==1) b=floor(j/2)+1 else b=j/2
for (k in 1:N)
{
if (data[k,(b+para+1)]<knot2[i,j]) data1[k,j]=0 else
data1[k,j]=data[k,(b+para+1)]-knot2[i,j]
}
}
mx=as.matrix(cbind(aa,data2,data1))
C=pinv(t(mx)%*%mx)
B=C%*%(t(mx)%*%data[,1])
yhat=mx%*%B
res=data[,1]-yhat
SSE[i]=sum((res)^2)
SSR[i]=sum((yhat-mean(data[,1]))^2)
MSE[i]=SSE[i]/(N)
Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100
A=mx%*%C%*%t(mx)
A1=(F-A)
A2=(sum(diag(A1))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll,file=paste0(folder,"DATA ALL knot 2.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 2 knot","\n")
cat("==============================================","\n")
print(((dataG[1,1:ncol(dataG)])))
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
# --- hitung ulang B dari titik knot dengan GCV TERKECIL ---
baris_terbaik = dataG[1, "knot_ke"]
knotgcv = knot2[baris_terbaik, ]
datagcv1=matrix(ncol=2*m,nrow=N)
for (j in 1:(2*m))
{
if (mod(j,2)==1) b=floor(j/2)+1 else b=j/2
for (k in 1:N)
{
if (data[k,(b+para+1)]<knotgcv[j]) datagcv1[k,j]=0 else
datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j]
}
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
Cgcv=pinv(t(mxgcv)%*%mxgcv)
B=Cgcv%*%(t(mxgcv)%*%data[,1])
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 2","\n")
cat("==============================================","\n")
print(B)
cat("\n")
}
GCV2(data)## ==============================================
## HASIL GCV terkecil dengan 2 knot
## ==============================================
## GCV knot_ke X1_k1 X1_k2 X2_k1 X2_k2 X3_k1
## 59.653741 245.000000 6.891224 22.295306 12.198367 52.859592 2.683265
## X3_k2 X4_k1 X4_k2 X5_k1 X5_k2 X6_k1 X6_k2
## 7.360816 10.631020 46.067755 58.439592 67.504898 6.146939 26.636735
## Nilai GCV 10 terkecil pertama
## GCV knot_ke X1_k1 X1_k2 X2_k1 X2_k2 X3_k1 X3_k2
## [1,] 59.65374 245 6.891224 22.29531 12.19837 52.85959 2.683265 7.360816
## [2,] 59.71889 239 6.891224 17.67408 12.19837 40.66122 2.683265 5.957551
## [3,] 59.71968 246 6.891224 23.06551 12.19837 54.89265 2.683265 7.594694
## [4,] 59.72621 244 6.891224 21.52510 12.19837 50.82653 2.683265 7.126939
## [5,] 59.72924 238 6.891224 16.90388 12.19837 38.62816 2.683265 5.723673
## [6,] 59.73756 905 23.065510 29.22714 54.89265 71.15714 7.594694 9.465714
## [7,] 59.76477 237 6.891224 16.13367 12.19837 36.59510 2.683265 5.489796
## [8,] 59.76563 241 6.891224 19.21449 12.19837 44.72735 2.683265 6.425306
## [9,] 59.76975 240 6.891224 18.44429 12.19837 42.69429 2.683265 6.191429
## [10,] 59.77079 242 6.891224 19.98469 12.19837 46.76041 2.683265 6.659184
## X4_k1 X4_k2 X5_k1 X5_k2 X6_k1 X6_k2
## [1,] 10.63102 46.06776 58.43959 67.50490 6.146939 26.63673
## [2,] 10.63102 35.43673 58.43959 64.78531 6.146939 20.48980
## [3,] 10.63102 47.83959 58.43959 67.95816 6.146939 27.66122
## [4,] 10.63102 44.29592 58.43959 67.05163 6.146939 25.61224
## [5,] 10.63102 33.66490 58.43959 64.33204 6.146939 19.46531
## [6,] 47.83959 62.01429 67.95816 71.58429 27.661224 35.85714
## [7,] 10.63102 31.89306 58.43959 63.87878 6.146939 18.44082
## [8,] 10.63102 38.98041 58.43959 65.69184 6.146939 22.53878
## [9,] 10.63102 37.20857 58.43959 65.23857 6.146939 21.51429
## [10,] 10.63102 40.75224 58.43959 66.14510 6.146939 23.56327
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 2
## ==============================================
## [,1]
## [1,] 232.61514814
## [2,] -0.21284620
## [3,] -0.47339109
## [4,] -3.37097474
## [5,] -0.93861870
## [6,] -2.77940373
## [7,] 0.78647431
## [8,] -0.04625130
## [9,] -1.12800600
## [10,] 0.33226449
## [11,] -0.04946566
## [12,] 6.15885749
## [13,] -4.00900376
## [14,] 0.87768543
## [15,] 0.02702941
## [16,] 4.11508131
## [17,] -0.43791153
## [18,] -0.91535404
## [19,] 0.05822887#UJI SIGNIFIKANSI DUA TITIK KNOT#
uji(0.05, 0, 2)## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Gagal tolak Ho yakni parameter (Intercept) tidak signifikan dengan pvalue 0.2033832
## Gagal tolak Ho yakni parameter X1 tidak signifikan dengan pvalue 0.5802154
## Gagal tolak Ho yakni parameter X1_k1 tidak signifikan dengan pvalue 0.9150597
## Tolak Ho yakni parameter X1_k2 signifikan dengan pvalue 0.0002143497
## Tolak Ho yakni parameter X2 signifikan dengan pvalue 0.001689131
## Tolak Ho yakni parameter X2_k1 signifikan dengan pvalue 0.04228767
## Gagal tolak Ho yakni parameter X2_k2 tidak signifikan dengan pvalue 0.580082
## Gagal tolak Ho yakni parameter X3 tidak signifikan dengan pvalue 0.53971
## Gagal tolak Ho yakni parameter X3_k1 tidak signifikan dengan pvalue 0.3171688
## Tolak Ho yakni parameter X3_k2 signifikan dengan pvalue 0.002221382
## Tolak Ho yakni parameter X4 signifikan dengan pvalue 4.68481e-11
## Tolak Ho yakni parameter X4_k1 signifikan dengan pvalue 7.791085e-06
## Gagal tolak Ho yakni parameter X4_k2 tidak signifikan dengan pvalue 0.9171649
## Gagal tolak Ho yakni parameter X5 tidak signifikan dengan pvalue 0.3726539
## Gagal tolak Ho yakni parameter X5_k1 tidak signifikan dengan pvalue 0.2028383
## Gagal tolak Ho yakni parameter X5_k2 tidak signifikan dengan pvalue 0.304175
## Gagal tolak Ho yakni parameter X6 tidak signifikan dengan pvalue 0.1014708
## Gagal tolak Ho yakni parameter X6_k1 tidak signifikan dengan pvalue 0.0736619
## Gagal tolak Ho yakni parameter X6_k2 tidak signifikan dengan pvalue 0.7026649
## =============================================
## nilai t hitung
## =============================================
## [,1]
## (Intercept) 1.2736558
## X1 -0.5534350
## X1_k1 -0.1067132
## X1_k2 -3.7291410
## X2 -3.1573696
## X2_k1 2.0359461
## X2_k2 -0.5536301
## X3 -0.6136780
## X3_k1 1.0013046
## X3_k2 -3.0749953
## X4 -6.7307292
## X4_k1 4.5186820
## X4_k2 0.1040587
## X5 -0.8923222
## X5_k1 1.2751953
## X5_k2 -1.0285894
## X6 1.6408161
## X6_k1 -1.7925095
## X6_k2 0.3819466
## =============================================
## Tabel uji parsial
## =============================================
## Estimasi Galat_baku t_hitung p_value Keputusan
## (Intercept) 232.61515 182.63580 1.2737 2.034e-01 Tidak signifikan
## X1 -0.21285 0.38459 -0.5534 5.802e-01 Tidak signifikan
## X1_k1 -0.04625 0.43342 -0.1067 9.151e-01 Tidak signifikan
## X1_k2 -1.12801 0.30248 -3.7291 2.143e-04 Signifikan
## X2 -0.47339 0.14993 -3.1574 1.689e-03 Signifikan
## X2_k1 0.33226 0.16320 2.0359 4.229e-02 Signifikan
## X2_k2 -0.04947 0.08935 -0.5536 5.801e-01 Tidak signifikan
## X3 -3.37097 5.49307 -0.6137 5.397e-01 Tidak signifikan
## X3_k1 6.15886 6.15083 1.0013 3.172e-01 Tidak signifikan
## X3_k2 -4.00900 1.30374 -3.0750 2.221e-03 Signifikan
## X4 -0.93862 0.13945 -6.7307 4.685e-11 Signifikan
## X4_k1 0.87769 0.19423 4.5187 7.791e-06 Signifikan
## X4_k2 0.02703 0.25975 0.1041 9.172e-01 Tidak signifikan
## X5 -2.77940 3.11480 -0.8923 3.727e-01 Tidak signifikan
## X5_k1 4.11508 3.22702 1.2752 2.028e-01 Tidak signifikan
## X5_k2 -0.43791 0.42574 -1.0286 3.042e-01 Tidak signifikan
## X6 0.78647 0.47932 1.6408 1.015e-01 Tidak signifikan
## X6_k1 -0.91535 0.51066 -1.7925 7.366e-02 Tidak signifikan
## X6_k2 0.05823 0.15245 0.3819 7.027e-01 Tidak signifikan
## Parameter signifikan: 6 dari 19
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 18 78383.06 4354.614 75.80013
## Error 495 28437.08 57.44864
## Total 513 106820.1
## =============================================
## s= 7.579488 Rsq= 73.37854 Rsq adj= 72.41049
## pvalue(F)= 1.905753e-129#PEMILIHAN TIGA TITIK KNOT#
Pencarian 3 knot memeriksa 17.296 kombinasi (48 pilih 3), sehingga chunk ini memerlukan waktu beberapa menit.
GCV3=function(data)
{
para=0
data=as.matrix(data)
N=nrow(data)
M=ncol(data)
m=ncol(data)-para-1
dataA=data[,(para+2):M]
dataA=as.matrix(dataA)
F=diag(N)
nk=50
knot1=matrix(ncol=m,nrow=nk)
for (i in (1:m))
{
a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk)
knot1[,i]=t(as.matrix(a))
}
a1=length(knot1[,1])
knot1=as.matrix(knot1[2:(a1-1),])
a2=nk-2
z=(a2*(a2-1)*(a2-2)/6)
knot2=cbind(rep(NA,(z+1)))
for (i in (1:m))
{
knot=rbind(rep(NA,3))
for ( j in 1:(a2-2))
{
for (k in (j+1):(a2-1))
{
for (g in (k+1):a2)
{
xx=cbind(knot1[j,i],knot1[k,i],knot1[g,i])
knot=rbind(knot,xx)
}
}
}
knot2=cbind(knot2,knot)
}
knot2=knot2[2:(z+1),2:(3*m+1)]
colnames(knot2)=paste0(rep(paste0("X",1:m),each=3),"_k",rep(1:3,m))
a3=nrow(knot2)
aa=rep(1,N)
data1=matrix(ncol=3*m,nrow=N)
data2=data[,2:M]
nk1=nrow(knot2)
GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV"
MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE"
SSE=rep(NA,nk1)
SSR=rep(NA,nk1)
Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq"
knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke"
for (i in 1:a3)
{
for (j in 1:(3*m))
{
b=ceiling(j/3)
for (k in 1:N)
{
if (data[k,(b+para+1)]<knot2[i,j]) data1[k,j]=0 else
data1[k,j]=data[k,(b+para+1)]-knot2[i,j]
}
}
mx=as.matrix(cbind(aa,data2,data1))
C=pinv(t(mx)%*%mx)
B=C%*%(t(mx)%*%data[,1])
yhat=mx%*%B
res=data[,1]-yhat
SSE[i]=sum((res)^2)
SSR[i]=sum((yhat-mean(data[,1]))^2)
MSE[i]=SSE[i]/(N)
Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100
A=mx%*%C%*%t(mx)
A1=(F-A)
A2=(sum(diag(A1))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll,file=paste0(folder,"DATA ALL knot 3.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 3 knot","\n")
cat("==============================================","\n")
print(((dataG[1,1:ncol(dataG)])))
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
# --- hitung ulang B dari titik knot dengan GCV TERKECIL ---
baris_terbaik = dataG[1, "knot_ke"]
knotgcv = knot2[baris_terbaik, ]
datagcv1=matrix(ncol=3*m,nrow=N)
for (j in 1:(3*m))
{
b=ceiling(j/3)
for (k in 1:N)
{
if (data[k,(b+para+1)]<knotgcv[j]) datagcv1[k,j]=0 else
datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j]
}
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
Cgcv=pinv(t(mxgcv)%*%mxgcv)
B=Cgcv%*%(t(mxgcv)%*%data[,1])
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 3","\n")
cat("==============================================","\n")
print(B)
cat("\n")
}
GCV3(data)## ==============================================
## HASIL GCV terkecil dengan 3 knot
## ==============================================
## GCV knot_ke X1_k1 X1_k2 X1_k3 X2_k1
## 57.871554 11127.000000 13.052857 24.605918 28.456939 28.462857
## X2_k2 X2_k3 X3_k1 X3_k2 X3_k3 X4_k1
## 58.958776 69.124082 4.554286 8.062449 9.231837 24.805714
## X4_k2 X4_k3 X5_k1 X5_k2 X5_k3 X6_k1
## 51.383265 60.242449 62.065714 68.864694 71.131020 14.342857
## X6_k2 X6_k3
## 29.710204 34.832653
## Nilai GCV 10 terkecil pertama
## GCV knot_ke X1_k1 X1_k2 X1_k3 X2_k1 X2_k2 X2_k3
## [1,] 57.87155 11127 13.05286 24.60592 28.45694 28.46286 58.95878 69.12408
## [2,] 57.90291 11144 13.05286 25.37612 27.68673 28.46286 60.99184 67.09102
## [3,] 57.91615 11126 13.05286 24.60592 27.68673 28.46286 58.95878 67.09102
## [4,] 57.97474 11655 13.82306 24.60592 28.45694 30.49592 58.95878 69.12408
## [5,] 58.00390 11672 13.82306 25.37612 27.68673 30.49592 60.99184 67.09102
## [6,] 58.01200 11654 13.82306 24.60592 27.68673 30.49592 58.95878 67.09102
## [7,] 58.01834 10566 12.28265 24.60592 28.45694 26.42980 58.95878 69.12408
## [8,] 58.05155 11145 13.05286 25.37612 28.45694 28.46286 60.99184 69.12408
## [9,] 58.07181 10583 12.28265 25.37612 27.68673 26.42980 60.99184 67.09102
## [10,] 58.08304 10565 12.28265 24.60592 27.68673 26.42980 58.95878 67.09102
## X3_k1 X3_k2 X3_k3 X4_k1 X4_k2 X4_k3 X5_k1 X5_k2
## [1,] 4.554286 8.062449 9.231837 24.80571 51.38327 60.24245 62.06571 68.86469
## [2,] 4.554286 8.296327 8.997959 24.80571 53.15510 58.47061 62.06571 69.31796
## [3,] 4.554286 8.062449 8.997959 24.80571 51.38327 58.47061 62.06571 68.86469
## [4,] 4.788163 8.062449 9.231837 26.57755 51.38327 60.24245 62.51898 68.86469
## [5,] 4.788163 8.296327 8.997959 26.57755 53.15510 58.47061 62.51898 69.31796
## [6,] 4.788163 8.062449 8.997959 26.57755 51.38327 58.47061 62.51898 68.86469
## [7,] 4.320408 8.062449 9.231837 23.03388 51.38327 60.24245 61.61245 68.86469
## [8,] 4.554286 8.296327 9.231837 24.80571 53.15510 60.24245 62.06571 69.31796
## [9,] 4.320408 8.296327 8.997959 23.03388 53.15510 58.47061 61.61245 69.31796
## [10,] 4.320408 8.062449 8.997959 23.03388 51.38327 58.47061 61.61245 68.86469
## X5_k3 X6_k1 X6_k2 X6_k3
## [1,] 71.13102 14.34286 29.71020 34.83265
## [2,] 70.67776 14.34286 30.73469 33.80816
## [3,] 70.67776 14.34286 29.71020 33.80816
## [4,] 71.13102 15.36735 29.71020 34.83265
## [5,] 70.67776 15.36735 30.73469 33.80816
## [6,] 70.67776 15.36735 29.71020 33.80816
## [7,] 71.13102 13.31837 29.71020 34.83265
## [8,] 71.13102 14.34286 30.73469 34.83265
## [9,] 70.67776 13.31837 30.73469 33.80816
## [10,] 70.67776 13.31837 29.71020 33.80816
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 3
## ==============================================
## [,1]
## [1,] 167.379165062
## [2,] 0.006721307
## [3,] -0.162048384
## [4,] 3.088515432
## [5,] -0.572073896
## [6,] -1.828642793
## [7,] 0.121548312
## [8,] -0.884363855
## [9,] 0.860349570
## [10,] -2.205508264
## [11,] -0.027020113
## [12,] -0.313063519
## [13,] 0.513304370
## [14,] -1.193279342
## [15,] -6.518365931
## [16,] 6.397219439
## [17,] 0.756445243
## [18,] 0.935416563
## [19,] -2.197128367
## [20,] 3.616134316
## [21,] -2.564442076
## [22,] 2.261039366
## [23,] -0.235273051
## [24,] -0.185090874
## [25,] 0.665478503#UJI SIGNIFIKANSI TIGA TITIK KNOT#
uji(0.05, 0, 3)## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Tolak Ho yakni parameter (Intercept) signifikan dengan pvalue 0.030933
## Gagal tolak Ho yakni parameter X1 tidak signifikan dengan pvalue 0.9589731
## Tolak Ho yakni parameter X1_k1 signifikan dengan pvalue 0.0008273402
## Gagal tolak Ho yakni parameter X1_k2 tidak signifikan dengan pvalue 0.3383506
## Gagal tolak Ho yakni parameter X1_k3 tidak signifikan dengan pvalue 0.05676414
## Tolak Ho yakni parameter X2 signifikan dengan pvalue 0.002971664
## Gagal tolak Ho yakni parameter X2_k1 tidak signifikan dengan pvalue 0.7632598
## Gagal tolak Ho yakni parameter X2_k2 tidak signifikan dengan pvalue 0.3196921
## Gagal tolak Ho yakni parameter X2_k3 tidak signifikan dengan pvalue 0.2179738
## Gagal tolak Ho yakni parameter X3 tidak signifikan dengan pvalue 0.1575939
## Gagal tolak Ho yakni parameter X3_k1 tidak signifikan dengan pvalue 0.6503021
## Tolak Ho yakni parameter X3_k2 signifikan dengan pvalue 0.0001412823
## Tolak Ho yakni parameter X3_k3 signifikan dengan pvalue 9.938082e-06
## Tolak Ho yakni parameter X4 signifikan dengan pvalue 1.323051e-09
## Tolak Ho yakni parameter X4_k1 signifikan dengan pvalue 0.002019186
## Gagal tolak Ho yakni parameter X4_k2 tidak signifikan dengan pvalue 0.4230511
## Gagal tolak Ho yakni parameter X4_k3 tidak signifikan dengan pvalue 0.1188693
## Gagal tolak Ho yakni parameter X5 tidak signifikan dengan pvalue 0.129938
## Tolak Ho yakni parameter X5_k1 signifikan dengan pvalue 0.00774821
## Tolak Ho yakni parameter X5_k2 signifikan dengan pvalue 0.0005331614
## Tolak Ho yakni parameter X5_k3 signifikan dengan pvalue 0.001533134
## Gagal tolak Ho yakni parameter X6 tidak signifikan dengan pvalue 0.4278882
## Gagal tolak Ho yakni parameter X6_k1 tidak signifikan dengan pvalue 0.2545297
## Gagal tolak Ho yakni parameter X6_k2 tidak signifikan dengan pvalue 0.6005636
## Gagal tolak Ho yakni parameter X6_k3 tidak signifikan dengan pvalue 0.1704717
## =============================================
## nilai t hitung
## =============================================
## [,1]
## (Intercept) 2.16419172
## X1 0.05146859
## X1_k1 -3.36441804
## X1_k2 0.95836847
## X1_k3 -1.90964221
## X2 -2.98564456
## X2_k1 -0.30137092
## X2_k2 -0.99610423
## X2_k3 1.23352155
## X3 1.41537923
## X3_k1 -0.45362202
## X3_k2 -3.83625662
## X3_k3 4.46530789
## X4 -6.18367326
## X4_k1 3.10410199
## X4_k2 0.80181212
## X4_k3 -1.56227849
## X5 -1.51690035
## X5_k1 2.67390671
## X5_k2 -3.48666848
## X5_k3 3.18631001
## X6 0.79347125
## X6_k1 -1.14075640
## X6_k2 -0.52393144
## X6_k3 1.37271011
## =============================================
## Tabel uji parsial
## =============================================
## Estimasi Galat_baku t_hitung p_value Keputusan
## (Intercept) 167.379165 77.34027 2.16419 3.093e-02 Signifikan
## X1 0.006721 0.13059 0.05147 9.590e-01 Tidak signifikan
## X1_k1 -0.884364 0.26286 -3.36442 8.273e-04 Signifikan
## X1_k2 0.860350 0.89772 0.95837 3.384e-01 Tidak signifikan
## X1_k3 -2.205508 1.15493 -1.90964 5.676e-02 Tidak signifikan
## X2 -0.162048 0.05428 -2.98564 2.972e-03 Signifikan
## X2_k1 -0.027020 0.08966 -0.30137 7.633e-01 Tidak signifikan
## X2_k2 -0.313064 0.31429 -0.99610 3.197e-01 Tidak signifikan
## X2_k3 0.513304 0.41613 1.23352 2.180e-01 Tidak signifikan
## X3 3.088515 2.18211 1.41538 1.576e-01 Tidak signifikan
## X3_k1 -1.193279 2.63056 -0.45362 6.503e-01 Tidak signifikan
## X3_k2 -6.518366 1.69915 -3.83626 1.413e-04 Signifikan
## X3_k3 6.397219 1.43265 4.46531 9.938e-06 Signifikan
## X4 -0.572074 0.09251 -6.18367 1.323e-09 Signifikan
## X4_k1 0.756445 0.24369 3.10410 2.019e-03 Signifikan
## X4_k2 0.935417 1.16663 0.80181 4.231e-01 Tidak signifikan
## X4_k3 -2.197128 1.40636 -1.56228 1.189e-01 Tidak signifikan
## X5 -1.828643 1.20551 -1.51690 1.299e-01 Tidak signifikan
## X5_k1 3.616134 1.35238 2.67391 7.748e-03 Signifikan
## X5_k2 -2.564442 0.73550 -3.48667 5.332e-04 Signifikan
## X5_k3 2.261039 0.70961 3.18631 1.533e-03 Signifikan
## X6 0.121548 0.15319 0.79347 4.279e-01 Tidak signifikan
## X6_k1 -0.235273 0.20624 -1.14076 2.545e-01 Tidak signifikan
## X6_k2 -0.185091 0.35327 -0.52393 6.006e-01 Tidak signifikan
## X6_k3 0.665479 0.48479 1.37271 1.705e-01 Tidak signifikan
## Parameter signifikan: 10 dari 25
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 24 79897.37 3329.057 60.46587
## Error 489 26922.77 55.05679
## Total 513 106820.1
## =============================================
## s= 7.420026 Rsq= 74.79617 Rsq adj= 73.55917
## pvalue(F)= 1.15666e-129#PERBANDINGAN GCV MINIMUM 1, 2, DAN 3 KNOT#
gcv_min=rep(NA,3)
r2_min=rep(NA,3)
for (K in 1:3)
{
d=read.csv(paste0(folder,"DATA ALL knot ",K,".csv"),header=TRUE)
d=d[order(d$GCV),]
gcv_min[K]=d$GCV[1]
r2_min[K]=d$Rsq[1]
}
tab_gcv=data.frame(Jumlah_knot=1:3,GCV_minimum=gcv_min,R2_persen=r2_min)
tab_gcvK_terbaik=tab_gcv$Jumlah_knot[which.min(tab_gcv$GCV_minimum)]
cat("Jumlah knot optimal (GCV terkecil):",K_terbaik,"\n")## Jumlah knot optimal (GCV terkecil): 3plot(tab_gcv$Jumlah_knot,tab_gcv$GCV_minimum,type="b",pch=16,col="steelblue",lwd=2,
xaxt="n",main="GCV minimum per jumlah knot",xlab="Jumlah knot",ylab="GCV")
axis(1,at=1:3)
points(K_terbaik,tab_gcv$GCV_minimum[K_terbaik],col="red",cex=2.5,lwd=2)
#MODEL SPLINE TERBAIK#
#fungsi untuk membentuk matriks basis spline truncated (intersep, lalu per variabel: Xj dan (Xj-knot)+)
buat_mx=function(Xm,knot)
{
Xm=as.matrix(Xm)
p=ncol(Xm)
K=nrow(knot)
mx=matrix(1,nrow=nrow(Xm),ncol=1)
nama="(Intercept)"
for (j in 1:p)
{
mx=cbind(mx,Xm[,j])
nama=c(nama,paste0("X",j))
for (h in 1:K)
{
mx=cbind(mx,pmax(Xm[,j]-knot[h,j],0))
nama=c(nama,paste0("X",j,"_k",h))
}
}
colnames(mx)=nama
mx
}
#fungsi untuk mengambil titik knot dengan GCV terkecil dari file csv (baris = urutan knot, kolom = variabel)
ambil_knot=function(K)
{
d=read.csv(paste0(folder,"DATA ALL knot ",K,".csv"),header=TRUE)
d=d[order(d$GCV),]
knot=matrix(as.numeric(d[1,5:ncol(d)]),nrow=K)
colnames(knot)=paste0("X",1:ncol(knot))
rownames(knot)=paste0("knot_",1:K)
knot
}
y=data[,1]
X=as.matrix(data[,2:7])
n=nrow(data)
knot_opt=ambil_knot(K_terbaik)
knot_opt## X1 X2 X3 X4 X5 X6
## knot_1 13.05286 28.46286 4.554286 24.80571 62.06571 14.34286
## knot_2 24.60592 58.95878 8.062449 51.38327 68.86469 29.71020
## knot_3 28.45694 69.12408 9.231837 60.24245 71.13102 34.83265mx_opt=buat_mx(X,knot_opt)
df_sp=as.data.frame(mx_opt[,-1])
df_sp$Y=y
model_sp=lm(Y~.,data=df_sp)
sm=summary(model_sp)
sm##
## Call:
## lm(formula = Y ~ ., data = df_sp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.742 -3.129 1.549 4.796 20.263
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.379165 77.340267 2.164 0.030933 *
## X1 0.006721 0.130590 0.051 0.958973
## X1_k1 -0.884364 0.262858 -3.364 0.000827 ***
## X1_k2 0.860350 0.897723 0.958 0.338351
## X1_k3 -2.205508 1.154933 -1.910 0.056764 .
## X2 -0.162048 0.054276 -2.986 0.002972 **
## X2_k1 -0.027020 0.089657 -0.301 0.763260
## X2_k2 -0.313064 0.314288 -0.996 0.319692
## X2_k3 0.513304 0.416129 1.234 0.217974
## X3 3.088515 2.182112 1.415 0.157594
## X3_k1 -1.193279 2.630559 -0.454 0.650302
## X3_k2 -6.518366 1.699148 -3.836 0.000141 ***
## X3_k3 6.397219 1.432649 4.465 9.94e-06 ***
## X4 -0.572074 0.092514 -6.184 1.32e-09 ***
## X4_k1 0.756445 0.243692 3.104 0.002019 **
## X4_k2 0.935417 1.166628 0.802 0.423051
## X4_k3 -2.197128 1.406362 -1.562 0.118869
## X5 -1.828643 1.205513 -1.517 0.129938
## X5_k1 3.616134 1.352379 2.674 0.007748 **
## X5_k2 -2.564442 0.735499 -3.487 0.000533 ***
## X5_k3 2.261039 0.709611 3.186 0.001533 **
## X6 0.121548 0.153186 0.793 0.427888
## X6_k1 -0.235273 0.206243 -1.141 0.254530
## X6_k2 -0.185091 0.353273 -0.524 0.600564
## X6_k3 0.665479 0.484792 1.373 0.170472
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.42 on 489 degrees of freedom
## Multiple R-squared: 0.748, Adjusted R-squared: 0.7356
## F-statistic: 60.47 on 24 and 489 DF, p-value: < 2.2e-16#tabel ANOVA
fs=sm$fstatistic
p_f=pf(fs[1],fs[2],fs[3],lower.tail=FALSE)
anova_tab=data.frame(Sumber=c("Regresi","Galat","Total"),
db=c(fs[2],fs[3],n-1),
JK=c(sum((fitted(model_sp)-mean(y))^2),sum(residuals(model_sp)^2),sum((y-mean(y))^2)))
anova_tab$KT=c(anova_tab$JK[1]/anova_tab$db[1],anova_tab$JK[2]/anova_tab$db[2],NA)
anova_tab$F_hitung=c(unname(fs[1]),NA,NA)
anova_tab$p_value=c(unname(p_f),NA,NA)
anova_tabcat(sprintf("R2 = %.2f%% | R2 adj = %.2f%% | s = %.4f\n",100*sm$r.squared,100*sm$adj.r.squared,sm$sigma))## R2 = 74.80% | R2 adj = 73.56% | s = 7.4200if (p_f<0.05) cat("Keputusan: tolak H0, minimal satu parameter signifikan.\n") else cat("Keputusan: gagal tolak H0.\n")## Keputusan: tolak H0, minimal satu parameter signifikan.#uji signifikansi parsial
koef=as.data.frame(sm$coefficients)
names(koef)=c("Estimasi","Galat_baku","t_hitung","p_value")
koef$Keputusan=ifelse(koef$p_value<0.05,"Signifikan","Tidak signifikan")
koefsig=koef[koef$p_value<0.05,]
cat("Parameter signifikan:",nrow(sig),"dari",nrow(koef),"\n")## Parameter signifikan: 10 dari 25#variabel dianggap berpengaruh jika minimal satu parameternya signifikan
nama_x=paste0("X",1:6)
var_pengaruh=nama_x[sapply(nama_x,function(v) any(grepl(paste0("^",v,"($|_k)"),rownames(sig))))]
cat("Variabel berpengaruh:",paste(var_pengaruh,collapse=", "),"\n")## Variabel berpengaruh: X1, X2, X3, X4, X5#estimasi parameter dengan selang kepercayaan 95%
kf=koef[rownames(koef)!="(Intercept)",]
kf$bawah=kf$Estimasi-1.96*kf$Galat_baku
kf$atas=kf$Estimasi+1.96*kf$Galat_baku
ypos=nrow(kf):1
warna=ifelse(kf$Keputusan=="Signifikan","#d73027","grey55")
plot(kf$Estimasi,ypos,xlim=range(c(kf$bawah,kf$atas)),pch=16,col=warna,yaxt="n",
xlab="Estimasi",ylab="",main="Estimasi parameter spline dengan selang kepercayaan 95%")
segments(kf$bawah,ypos,kf$atas,ypos,col=warna)
axis(2,at=ypos,labels=rownames(kf),las=1,cex.axis=0.7)
abline(v=0,lty=2)
legend("bottomright",legend=c("Signifikan","Tidak signifikan"),col=c("#d73027","grey55"),pch=16,bty="n")
#DIAGNOSTIK MODEL SPLINE#
res=residuals(model_sp)
hat=fitted(model_sp)
par(mfrow=c(2,2))
plot(y,hat,pch=16,col="grey30",main="Aktual vs prediksi",xlab="Y aktual",ylab="Y prediksi")
abline(0,1,col="red",lwd=2)
plot(hat,res,pch=16,col="grey30",main="Residual vs prediksi",xlab="Prediksi",ylab="Residual")
abline(h=0,col="red",lwd=2)
lines(lowess(hat,res),col="steelblue",lwd=2)
qqnorm(rstandard(model_sp),pch=16,col="grey30",main="Q-Q plot residual baku")
qqline(rstandard(model_sp),col="red",lwd=2)
hist(res,breaks=30,freq=FALSE,col="steelblue",border="white",main="Histogram residual",xlab="Residual")
curve(dnorm(x,mean=0,sd=sd(res)),add=TRUE,col="red",lwd=2)
par(mfrow=c(1,1))#uji asumsi residual model spline
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.93503, p-value = 3.613e-14bptest(model_sp)##
## studentized Breusch-Pagan test
##
## data: model_sp
## BP = 87.693, df = 24, p-value = 3.461e-09dwtest(model_sp)##
## Durbin-Watson test
##
## data: model_sp
## DW = 1.5473, p-value = 3.331e-08
## alternative hypothesis: true autocorrelation is greater than 0#PENGAMATAN BERPENGARUH (COOK’S DISTANCE)#
cd=cooks.distance(model_sp)
batas=4/n
plot(cd,type="h",col="grey50",main="Cook's distance (garis merah = 4/n)",xlab="Observasi",ylab="Cook's D")
abline(h=batas,col="red",lty=2)
besar=which(cd>4*batas)
if (length(besar)>0) text(besar,cd[besar],labels=besar,pos=3,cex=0.7)
d_cd=data.frame(obs=1:n,cd=cd)
head(d_cd[order(-d_cd$cd),],8)#VISUALISASI KURVA SPLINE PER VARIABEL#
Untuk tiap variabel, kurva menunjukkan efek parsial (variabel lain ditetapkan pada median). Titik adalah partial residual (residual + efek parsial variabel tersebut), garis putus-putus vertikal adalah letak knot.
med=apply(X,2,median)
b_sp=coef(model_sp)
b_sp[is.na(b_sp)]=0
efek_parsial=function(xj_nilai,j)
{
Xb=matrix(rep(med,each=length(xj_nilai)),ncol=ncol(X))
Xb[,j]=xj_nilai
drop(buat_mx(Xb,knot_opt)%*%b_sp)
}
par(mfrow=c(3,2))
for (j in 1:6)
{
g=seq(min(X[,j]),max(X[,j]),length.out=300)
pr=res+efek_parsial(X[,j],j)
plot(X[,j],pr,pch=16,col=rgb(0,0,0,0.3),cex=0.6,
main=paste("Efek parsial spline truncated",nama_x[j],"(",K_terbaik,"knot )"),
xlab=nama_x[j],ylab="Efek parsial pada Y")
lines(g,efek_parsial(g,j),col="darkred",lwd=3)
abline(v=knot_opt[,j],lty=2,col="steelblue")
}
par(mfrow=c(1,1))#PERSAMAAN PIECEWISE DAN PENGELOMPOKAN OBSERVASI#
for (j in 1:6)
{
nm=nama_x[j]
kj=knot_opt[,j]
b_lin=b_sp[nm]
b_h=b_sp[paste0(nm,"_k",1:K_terbaik)]
bawah=c(-Inf,kj)
atas=c(kj,Inf)
tab=data.frame(Segmen=1:(K_terbaik+1),
Batas_bawah=bawah,
Batas_atas=atas,
Kemiringan=unname(b_lin+c(0,cumsum(b_h))),
Konstanta=unname(b_sp["(Intercept)"]+c(0,-cumsum(b_h*kj))),
N_obs=sapply(1:(K_terbaik+1),function(r) sum(X[,j]>=bawah[r] & X[,j]<atas[r])))
cat("\n==============================================\n")
cat("Model per segmen untuk",nm,"(variabel lain dianggap konstan)\n")
cat("==============================================\n")
print(tab,digits=4,row.names=FALSE)
for (r in 1:nrow(tab))
{
arah=ifelse(tab$Kemiringan[r]>=0,"meningkatkan","menurunkan")
cat("- Pada segmen ",r," (",sprintf("%.3f",tab$Batas_bawah[r])," <= ",nm," < ",sprintf("%.3f",tab$Batas_atas[r]),
"), setiap kenaikan 1 satuan ",nm," ",arah," Y sebesar ",sprintf("%.3f",abs(tab$Kemiringan[r])),
" (",tab$N_obs[r]," observasi).\n",sep="")
}
}##
## ==============================================
## Model per segmen untuk X1 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 13.05 0.006721 167.4 360
## 2 13.05 24.61 -0.877643 178.9 117
## 3 24.61 28.46 -0.017293 157.8 14
## 4 28.46 Inf -2.222801 220.5 23
## - Pada segmen 1 (-Inf <= X1 < 13.053), setiap kenaikan 1 satuan X1 meningkatkan Y sebesar 0.007 (360 observasi).
## - Pada segmen 2 (13.053 <= X1 < 24.606), setiap kenaikan 1 satuan X1 menurunkan Y sebesar 0.878 (117 observasi).
## - Pada segmen 3 (24.606 <= X1 < 28.457), setiap kenaikan 1 satuan X1 menurunkan Y sebesar 0.017 (14 observasi).
## - Pada segmen 4 (28.457 <= X1 < Inf), setiap kenaikan 1 satuan X1 menurunkan Y sebesar 2.223 (23 observasi).
##
## ==============================================
## Model per segmen untuk X2 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 28.46 -0.16205 167.4 300
## 2 28.46 58.96 -0.18907 168.1 174
## 3 58.96 69.12 -0.50213 186.6 15
## 4 69.12 Inf 0.01117 151.1 25
## - Pada segmen 1 (-Inf <= X2 < 28.463), setiap kenaikan 1 satuan X2 menurunkan Y sebesar 0.162 (300 observasi).
## - Pada segmen 2 (28.463 <= X2 < 58.959), setiap kenaikan 1 satuan X2 menurunkan Y sebesar 0.189 (174 observasi).
## - Pada segmen 3 (58.959 <= X2 < 69.124), setiap kenaikan 1 satuan X2 menurunkan Y sebesar 0.502 (15 observasi).
## - Pada segmen 4 (69.124 <= X2 < Inf), setiap kenaikan 1 satuan X2 meningkatkan Y sebesar 0.011 (25 observasi).
##
## ==============================================
## Model per segmen untuk X3 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 4.554 3.089 167.4 10
## 2 4.554 8.062 1.895 172.8 129
## 3 8.062 9.232 -4.623 225.4 192
## 4 9.232 Inf 1.774 166.3 183
## - Pada segmen 1 (-Inf <= X3 < 4.554), setiap kenaikan 1 satuan X3 meningkatkan Y sebesar 3.089 (10 observasi).
## - Pada segmen 2 (4.554 <= X3 < 8.062), setiap kenaikan 1 satuan X3 meningkatkan Y sebesar 1.895 (129 observasi).
## - Pada segmen 3 (8.062 <= X3 < 9.232), setiap kenaikan 1 satuan X3 menurunkan Y sebesar 4.623 (192 observasi).
## - Pada segmen 4 (9.232 <= X3 < Inf), setiap kenaikan 1 satuan X3 meningkatkan Y sebesar 1.774 (183 observasi).
##
## ==============================================
## Model per segmen untuk X4 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 24.81 -0.5721 167.4 497
## 2 24.81 51.38 0.1844 148.6 13
## 3 51.38 60.24 1.1198 100.6 2
## 4 60.24 Inf -1.0773 232.9 2
## - Pada segmen 1 (-Inf <= X4 < 24.806), setiap kenaikan 1 satuan X4 menurunkan Y sebesar 0.572 (497 observasi).
## - Pada segmen 2 (24.806 <= X4 < 51.383), setiap kenaikan 1 satuan X4 meningkatkan Y sebesar 0.184 (13 observasi).
## - Pada segmen 3 (51.383 <= X4 < 60.242), setiap kenaikan 1 satuan X4 meningkatkan Y sebesar 1.120 (2 observasi).
## - Pada segmen 4 (60.242 <= X4 < Inf), setiap kenaikan 1 satuan X4 menurunkan Y sebesar 1.077 (2 observasi).
##
## ==============================================
## Model per segmen untuk X5 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 62.07 -1.829 167.38 10
## 2 62.07 68.86 1.787 -57.06 161
## 3 68.86 71.13 -0.777 119.54 128
## 4 71.13 Inf 1.484 -41.29 215
## - Pada segmen 1 (-Inf <= X5 < 62.066), setiap kenaikan 1 satuan X5 menurunkan Y sebesar 1.829 (10 observasi).
## - Pada segmen 2 (62.066 <= X5 < 68.865), setiap kenaikan 1 satuan X5 meningkatkan Y sebesar 1.787 (161 observasi).
## - Pada segmen 3 (68.865 <= X5 < 71.131), setiap kenaikan 1 satuan X5 menurunkan Y sebesar 0.777 (128 observasi).
## - Pada segmen 4 (71.131 <= X5 < Inf), setiap kenaikan 1 satuan X5 meningkatkan Y sebesar 1.484 (215 observasi).
##
## ==============================================
## Model per segmen untuk X6 (variabel lain dianggap konstan)
## ==============================================
## Segmen Batas_bawah Batas_atas Kemiringan Konstanta N_obs
## 1 -Inf 14.34 0.1215 167.4 83
## 2 14.34 29.71 -0.1137 170.8 326
## 3 29.71 34.83 -0.2988 176.3 65
## 4 34.83 Inf 0.3667 153.1 40
## - Pada segmen 1 (-Inf <= X6 < 14.343), setiap kenaikan 1 satuan X6 meningkatkan Y sebesar 0.122 (83 observasi).
## - Pada segmen 2 (14.343 <= X6 < 29.710), setiap kenaikan 1 satuan X6 menurunkan Y sebesar 0.114 (326 observasi).
## - Pada segmen 3 (29.710 <= X6 < 34.833), setiap kenaikan 1 satuan X6 menurunkan Y sebesar 0.299 (65 observasi).
## - Pada segmen 4 (34.833 <= X6 < Inf), setiap kenaikan 1 satuan X6 meningkatkan Y sebesar 0.367 (40 observasi).#daftar nomor observasi per segmen disimpan ke CSV (agar laporan tidak penuh angka)
id_segmen=do.call(rbind,lapply(1:6,function(j)
{
seg=findInterval(X[,j],knot_opt[,j])+1
data.frame(Variabel=nama_x[j],Segmen=seg,Observasi=1:n)
}))
write.csv(id_segmen,paste0(folder,"pengelompokan_observasi_per_segmen.csv"),row.names=FALSE)
#banyak observasi per segmen (baris = variabel, kolom = nomor segmen)
table(id_segmen$Variabel,id_segmen$Segmen)##
## 1 2 3 4
## X1 360 117 14 23
## X2 300 174 15 25
## X3 10 129 192 183
## X4 497 13 2 2
## X5 10 161 128 215
## X6 83 326 65 40#PERBANDINGAN MODEL#
hitung_gcv=function(model)
{
sse=sum(residuals(model)^2)
p=sum(!is.na(coef(model)))
(sse/n)/((n-p)/n)^2
}
fit_spline=function(K)
{
mxK=buat_mx(X,ambil_knot(K))
dfK=as.data.frame(mxK[,-1])
dfK$Y=y
lm(Y~.,data=dfK)
}
semua_model=list(model,fit_spline(1),fit_spline(2),fit_spline(3))
bandingkan=data.frame(Model=c("Regresi linear",paste0("Spline ",1:3," knot")),
Parameter=sapply(semua_model,function(m) length(coef(m))),
R2=NA,RMSE=NA,GCV=NA)
for (i in 1:4)
{
m=semua_model[[i]]
bandingkan$R2[i]=100*summary(m)$r.squared
bandingkan$RMSE[i]=sqrt(mean(residuals(m)^2))
bandingkan$GCV[i]=hitung_gcv(m)
}
bandingkanpar(mfrow=c(1,2))
barplot(bandingkan$R2,names.arg=bandingkan$Model,col="steelblue",main="R2 (%)",las=2,cex.names=0.7)
barplot(bandingkan$RMSE,names.arg=bandingkan$Model,col="steelblue",main="RMSE",las=2,cex.names=0.7)
par(mfrow=c(1,1))#KESIMPULAN#
Isi bagian ini setelah file di-knit. Poin yang bisa diisi:
- Jumlah knot optimal dan nilai GCV minimumnya (tabel perbandingan GCV).
- Variabel yang berpengaruh signifikan (bagian uji parsial) beserta R2 model.
- Perbandingan model spline dengan regresi linear (bagian perbandingan model).
- Hasil uji asumsi residual dan pengamatan berpengaruh.
#INFORMASI SESI#
sessionInfo()## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Asia/Singapore
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Matrix_1.7-6 pracma_2.4.6 pastecs_1.4.2 car_3.1-5 carData_3.0-6
## [6] MASS_7.3-66 lmtest_0.9-40 zoo_1.9-0
##
## loaded via a namespace (and not attached):
## [1] vctrs_0.7.3 cli_3.6.6 knitr_1.52 rlang_1.3.0
## [5] xfun_0.60 Formula_1.2-6 jsonlite_2.0.0 htmltools_0.5.9
## [9] sass_0.4.10 rmarkdown_2.32 grid_4.6.1 evaluate_1.0.5
## [13] jquerylib_0.1.4 abind_1.4-8 fastmap_1.2.0 yaml_2.3.12
## [17] lifecycle_1.0.5 compiler_4.6.1 rstudioapi_0.19.0 lattice_0.22-9
## [21] digest_0.6.39 R6_2.6.1 bslib_0.12.0 tools_4.6.1
## [25] boot_1.3-32 cachem_1.1.0