Praktikum 13-14 APG : Analisis Diskriminan

Terkadang Ambisi Tinggi dapat Menjatuhkan ke Jurang Terdalam, Lakukan versi terbaik dan Tetaplah Bersyukur

Pengantar

Summary AD

Analisis diskriminan: untuk mendeskripsikan atau menjelaskan perbedaan antara dua atau lebih kelompok. Tujuannya adalah untuk menemukan jumlah dimensi terendah yang diperlukan untuk menggambarkan perbedaan kelompok.
Analisis klasifikasi: untuk mengurutkan objek (observasi) menjadi dua atau lebih kelas berlabel. Penekanannya adalah mendapatkan aturan yang dapat digunakan untuk menetapkan objek baru secara optimal ke kelas berlabel.
Analisis klasifikasi merupakan aspek prediksi dari analisis diskriminan yang merupakan aspek deskriptif.

KONSEP DASAR ANALISIS DISKRIMINAN

Analisis diskriminan berusaha menemukan sekumpulan fungsi yang dapat digunakan untuk memisahkan pengamatan ke dalam kelompok yang diketahui
Untuk menggunakan prosedur ini beberapa elemen harus diketahui sebelum analisis:
Jumlah grup harus diketahui.
Harus ada kumpulan data “pelatihan” dengan indikator keanggotaan kelompok untuk setiap mata pelajaran.
Setelah analisis tersebut diselesaikan, seseorang kemudian akan dapat mengklasifikasikan pengamatan baru tanpa mengetahui keanggotaan kelompok secara apriori.

ASUMSI ANALISIS DISKRIMINAN

normalitas multivariat,
independensi kasus,
homogenitas kovarian kelompok

Penerapan 1 Suhu Baja (Buku Alvin)

$y_1$ = titik luluh
$y_2$ = ultimate strength

Load data

y1<-c(33,36,35,38,40,35,36,38,39,41,43,41)
y2<-c(60,61,64,63,65,57,59,59,61,63,65,59)
grup<-c(1,1,1,1,1,2,2,2,2,2,2,2)
data1<-data.frame(y1,y2,grup)
names(data1)<-c("y1","y2","group")
data1$group<-as.factor(data1$group)
data1

##    y1 y2 group
## 1  33 60     1
## 2  36 61     1
## 3  35 64     1
## 4  38 63     1
## 5  40 65     1
## 6  35 57     2
## 7  36 59     2
## 8  38 59     2
## 9  39 61     2
## 10 41 63     2
## 11 43 65     2
## 12 41 59     2

Menentukan nilai a

untuk memaksimalkan fungsi diskriminan $z=a^Ty$ maka nilai dapat dicari dengan
$a=(S_{pl})^-1\times (\bar y_1 -\bar y_2)$
Sedangkan S pooled dapat dihitung sebagai berikut $S_{pl}=\frac {1}{n_1+n_2-1}\times(W_1+W_2)$ dengan $W_1=(n_1-1)S_1$ dan $W_2=(n_2-1)S_2$

library(dplyr)
datat1<-data1 %>%filter(group==1) %>% select(y1,y2)
datat1

##   y1 y2
## 1 33 60
## 2 36 61
## 3 35 64
## 4 38 63
## 5 40 65

datat2<-data1 %>%filter(group==2) %>% select(y1,y2)
datat2

##   y1 y2
## 1 35 57
## 2 36 59
## 3 38 59
## 4 39 61
## 5 41 63
## 6 43 65
## 7 41 59

S1<-cov(datat1)
S1

##     y1  y2
## y1 7.3 4.2
## y2 4.2 4.3

S2<-cov(datat2)
S2

##          y1       y2
## y1 8.333333 6.666667
## y2 6.666667 7.619048

# Matriks rata-rata untuk setiap grup
ybar1<-apply(datat1[,1:2],2,mean)
ybar2<-apply(datat2[,1:2],2,mean)
ybar1

##   y1   y2 
## 36.4 62.6

ybar2

##       y1       y2 
## 39.00000 60.42857

# Menentukan Spooled n1=5,n2=7
Spooled<-function(S1,S2,n1,n2){
  W1<-(n1-1)*S1
  W2<-(n2-1)*S2
  Sp<-(1/(n1+n2-2))*(W1+W2)
  return(Sp)
}
Spl<-Spooled(S1,S2,5,7)
Spl

##      y1       y2
## y1 7.92 5.680000
## y2 5.68 6.291429

# Menentukan nilai a
a<-solve(Spl)%*%(ybar1-ybar2)
a

##         [,1]
## y1 -1.633377
## y2  1.819779

Interpretasi Hasil menunjukkan untuk nilai a adalah -1.633 dan 1.820, maka dari itu fungsi diskriminan nya adalah sebagai berikut $z=a^T y=-1.633 y_1 +1.819 y_2$

Penerapan 2

soal pertama

Empat pengukuran dilakukan pada dua spesies kutu kumbang (Lubischew 1962). jenis kumbang adalah

jenis kumbangnya 1: Haltica Oleracea 2: Haltica Carduorum

x1 = jarak alur melintang dari batas posterior prothorax (μm), x2 = panjang elytra (dalam 0,01 mm), x3 = panjang sambungan antena kedua (μm), x4 = panjang sambungan antena ketiga (μm).

Pertanyaan: (a) Tentukan vektor koefisien fungsi diskriminan. (b) Hitung uji-t untuk masing-masing variabel.

#Anggapannya kita tidak membagi datanya untuk testing dan training

library(readxl)
library(dplyr)

(datameasur <- read_excel("Pertemuan 13_2.xlsx"))#masukkan datanya

## # A tibble: 39 x 5
##       x1    x2    x3    x4 Haltica
##    <dbl> <dbl> <dbl> <dbl>   <dbl>
##  1   189   245   137   163       1
##  2   192   260   132   217       1
##  3   217   276   141   192       1
##  4   221   299   142   213       1
##  5   171   239   128   158       1
##  6   192   262   147   173       1
##  7   213   278   136   201       1
##  8   192   255   128   185       1
##  9   170   244   146   192       1
## 10   201   276   128   186       1
## # ... with 29 more rows

(oleracea<-filter(datameasur[1:4],datameasur$Haltica=="1"))#Partisi hanya jenis haltica oleracea

## # A tibble: 19 x 4
##       x1    x2    x3    x4
##    <dbl> <dbl> <dbl> <dbl>
##  1   189   245   137   163
##  2   192   260   132   217
##  3   217   276   141   192
##  4   221   299   142   213
##  5   171   239   128   158
##  6   192   262   147   173
##  7   213   278   136   201
##  8   192   255   128   185
##  9   170   244   146   192
## 10   201   276   128   186
## 11   195   242   147   192
## 12   205   263   121   192
## 13   180   252   121   167
## 14   192   283   138   183
## 15   200   294   138   188
## 16   192   277   150   177
## 17   200   287   136   173
## 18   181   255   146   183
## 19   192   287   141   198

(calduarum<-filter(datameasur[1:4],datameasur$Haltica=="2"))#Partisi hanya jenis haltica calduarum

## # A tibble: 20 x 4
##       x1    x2    x3    x4
##    <dbl> <dbl> <dbl> <dbl>
##  1   181   305   184   209
##  2   158   237   133   188
##  3   184   300   166   231
##  4   171   273   162   213
##  5   181   297   163   224
##  6   181   308   160   223
##  7   177   301   166   221
##  8   198   308   141   197
##  9   180   286   146   214
## 10   177   299   171   192
## 11   176   317   166   213
## 12   192   312   166   209
## 13   176   285   141   200
## 14   169   287   162   214
## 15   164   265   147   192
## 16   181   308   157   204
## 17   192   276   154   209
## 18   181   278   149   235
## 19   175   271   140   192
## 20   197   303   170   205

(n1 <- nrow(oleracea))#Jumlah Obs dari haltica olearacea

## [1] 19

(n2 <-nrow(calduarum))#Jumlah Obs dari haltica calduarum

## [1] 20

Uji Asumsi Multivariate Normal

library(MVN) # Uji multivariate normal jika skala datanya numerik 
mvn(data = datameasur[,c(1:4)], multivariatePlot = "qq")

## $multivariateNormality
##            Test        HZ   p value MVN
## 1 Henze-Zirkler 0.8838133 0.1016201 YES
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling    x1        0.4975    0.1997    YES   
## 2 Anderson-Darling    x2        0.4363    0.2832    YES   
## 3 Anderson-Darling    x3        0.5137    0.1817    YES   
## 4 Anderson-Darling    x4        0.2497    0.7287    YES   
## 
## $Descriptives
##     n     Mean  Std.Dev Median Min Max  25th  75th       Skew   Kurtosis
## x1 39 186.8205 14.03168    184 158 221 177.0 193.5  0.3790264 -0.1856496
## x2 39 279.2308 22.42116    278 237 317 262.5 299.0 -0.2462722 -1.0623505
## x3 39 147.3590 14.98330    146 121 184 137.5 161.0  0.3126021 -0.6577109
## x4 39 197.8974 18.48868    197 158 235 187.0 213.0 -0.1040536 -0.6351910

Dari grafik Chi-Square QQ Plot diatas, dapat dilihat bahwasanya secara umum terbentuk garis linear X = Y, maka dapat dikatakan bahwa data berdistribusi multivariate normal.

Matriks Ragam-peragam antar kategori sama

library(biotools)
#melakukan uji statistik Box’s M di R dapat menggunakan fungsi boxM.
boxM(data = datameasur[, c(1:4)], grouping = datameasur$Haltica)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datameasur[, c(1:4)]
## Chi-Sq (approx.) = 7.3838, df = 10, p-value = 0.6888

Output diatas menunjukkan bahwa dengan tingkat signifikansi 5%, didapat keputusan untuk gagal tolak hipotesis nol sehingga tidak terdapat perbedaan matriks ragam peragamnya dan dapat digunakan model diskriminan linear

Terdapat perbedaan rata-rata antar kategori H0 : Tidak terdapat perbedaan rata-rata skor diskriminan antar kelompok secara multivariat H1 : terdapat perbedaan rata-rata skor diskriminan antar kelompok secara multivariat

(wilk <- manova(formula = cbind(datameasur$x1, datameasur$x2, datameasur$x3,datameasur$x4)~datameasur$Haltica))

## Call:
##    manova(formula = cbind(datameasur$x1, datameasur$x2, datameasur$x3, 
##     datameasur$x4) ~ datameasur$Haltica)
## 
## Terms:
##                 datameasur$Haltica Residuals
## resp 1                    2170.057  5311.687
## resp 2                    5494.776 13608.147
## resp 3                    3975.774  4555.200
## resp 4                    5290.892  7698.697
## Deg. of Freedom                  1        37
## 
## Residual standard errors: 11.98162 19.17779 11.09565 14.42473
## Estimated effects may be unbalanced

summary(object = wilk, test = "Wilks")

##                    Df   Wilks approx F num Df den Df    Pr(>F)    
## datameasur$Haltica  1 0.22512   29.258      4     34 1.388e-10 ***
## Residuals          37                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

otput diatas menunjukkan bahwa dengan tingkat signifikansi 5%, didapat keputusan untuk menolak hipotesis nol atau dengan kata lain terdapat perbedaan rata-rata (nilai variabel prediktor) antar kategori

Selain itu, untuk memeriksa korelasi dan multikolinieritas antar variabel independen dapat menjalankan perintah berikut

#Correlation dan multikolinearity check
kor<-cor(datameasur[,1:4])
kor

##             x1        x2         x3          x4
## x1  1.00000000 0.1809792 -0.2735568 -0.07351397
## x2  0.18097923 1.0000000  0.6452198  0.59184021
## x3 -0.27355681 0.6452198  1.0000000  0.58977338
## x4 -0.07351397 0.5918402  0.5897734  1.00000000

Kemudian dilakukan canonical discriminant analysis (analisis diskriminan kanonik) untuk melihat hubungan antara variabel diskriminan dengan variabel independen secara multivariat

library(candisc)
#Canonical Correlation (CC)
cc<-candisc(wilk)
cc

## 
## Canonical Discriminant Analysis for datameasur$Haltica:
## 
##    CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.77488     3.4421                100        100
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F numDF denDF   Pr(> F)    
## 1      0.22512   29.258     4    34 1.388e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Eigenvalue menunjukan ada tidaknya multikolinearitas antar variabel. Didapatkan nilai Eigenvalue = 3.4421. Dikarenakan nilai eigenvalue tidak mendekati 0, maka tidak terjadi multikolinearitas antar variabel independen. #Canonical Corellation Canonical Correlation (CC)=0.77488, artinya bahwa hubungan antara variabel diskriminan dengan variabel independen secara multivariat sebesar 077488 atau besarnya (CC)²= (0.77488)²= 0,600.

#Jadi dapat disimpulkan bahwa 60% variasi antara Haltica Oleracea dan Haltica Carduorum dapat dijelaskan oleh variabel X1, x2,x3 dan x4).

Memulai Analisis Diskriminan

secara manual Hitung rata-rata

(xbaroleracea <-colMeans(oleracea))

##       x1       x2       x3       x4 
## 194.4737 267.0526 137.0000 185.9474

(xbarcalduarum <-colMeans(calduarum))

##     x1     x2     x3     x4 
## 179.55 290.80 157.20 209.25

Hitung Spool

(s.pool<-((n1-1)*cov(datameasur[,1:4])+(n2-
1)*cov(datameasur[,1:4]))/(n1+n2-2))

##           x1        x2        x3        x4
## x1 196.88799  56.93725 -57.51282 -19.07152
## x2  56.93725 502.70850 216.75709 245.34008
## x3 -57.51282 216.75709 224.49933 163.37989
## x4 -19.07152 245.34008 163.37989 341.83131

Hitung Koefisien a

(koef.a<-solve(s.pool)%*%(xbaroleracea-xbarcalduarum))

##           [,1]
## x1  0.07925472
## x2 -0.03501070
## x3 -0.01190657
## x4 -0.03292940

###Menghitung cut off/ cuting score###

(m<-0.5*((xbaroleracea-xbarcalduarum)%*%solve(s.pool)%*%(xbaroleracea+xbarcalduarum)))

##           [,1]
## [1,] -3.202098

sehingga di peroleh persamaan pembeda formulasimya ada di buku jhonson ya

y1 <- t(koef.a)%*%xbaroleracea
y1

##           [,1]
## [1,] -1.691078

y2 <-t(koef.a)%*%xbarcalduarum
y2

##           [,1]
## [1,] -4.713117

Analisis Diskriminan Dengan Menggunakan Package membentuk fungsi diskriminan

Linear Discriminant An. with Jacknifed Prediction

library(MASS)
fit <- lda(Haltica~., data =  datameasur)
fit

## Call:
## lda(Haltica ~ ., data = datameasur)
## 
## Prior probabilities of groups:
##         1         2 
## 0.4871795 0.5128205 
## 
## Group means:
##         x1       x2    x3       x4
## 1 194.4737 267.0526 137.0 185.9474
## 2 179.5500 290.8000 157.2 209.2500
## 
## Coefficients of linear discriminants:
##            LD1
## x1 -0.09481540
## x2  0.04188462
## x3  0.01424428
## x4  0.03939467

berdasarkan hasil output, maka persamaan yang terbentuk adalah terbentuk satu persamaan : \[ y cap = -0.09481540x1+0.04188462x2+0.01424428x3+0.03939467x4 \] Kemudian akan dilakukan prediksi untuk menentukan kelas peluang anggota setiap kelas dan diskriminan linear, serta membandingkan data asli dengan data prediksi. dapat dijalankan dengan syntax sebagai berikut :

#membuat prediksi 
predik <-predict(fit,datameasur)
predik

## $class
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2
## [39] 2
## Levels: 1 2
## 
## $posterior
##               1            2
## 1  9.999865e-01 1.346284e-05
## 2  9.268287e-01 7.317130e-02
## 3  9.999924e-01 7.632217e-06
## 4  9.986802e-01 1.319767e-03
## 5  9.991995e-01 8.005012e-04
## 6  9.995609e-01 4.390736e-04
## 7  9.998866e-01 1.133694e-04
## 8  9.996841e-01 3.159040e-04
## 9  5.646499e-01 4.353501e-01
## 10 9.995995e-01 4.004713e-04
## 11 9.998862e-01 1.137660e-04
## 12 9.999767e-01 2.327413e-05
## 13 9.993414e-01 6.586304e-04
## 14 9.731329e-01 2.686705e-02
## 15 9.811852e-01 1.881477e-02
## 16 9.912920e-01 8.707966e-03
## 17 9.992928e-01 7.071807e-04
## 18 9.745908e-01 2.540915e-02
## 19 6.666361e-01 3.333639e-01
## 20 6.877758e-05 9.999312e-01
## 21 1.744691e-01 8.255309e-01
## 22 4.515096e-05 9.999548e-01
## 23 4.984444e-04 9.995016e-01
## 24 8.042812e-05 9.999196e-01
## 25 2.046323e-05 9.999795e-01
## 26 1.463367e-05 9.999854e-01
## 27 4.285410e-01 5.714590e-01
## 28 3.002200e-03 9.969978e-01
## 29 9.515822e-04 9.990484e-01
## 30 2.878003e-06 9.999971e-01
## 31 2.607210e-03 9.973928e-01
## 32 8.378414e-03 9.916216e-01
## 33 2.615259e-05 9.999738e-01
## 34 6.507909e-03 9.934921e-01
## 35 3.574139e-04 9.996426e-01
## 36 5.305999e-01 4.694001e-01
## 37 6.129373e-04 9.993871e-01
## 38 1.411753e-01 8.588247e-01
## 39 7.541599e-02 9.245840e-01
## 
## $x
##            LD1
## 1  -3.16272092
## 2  -0.76280678
## 3  -3.31970615
## 4  -1.89408906
## 5  -2.03252338
## 6  -2.19873905
## 7  -2.57334465
## 8  -2.28983658
## 9  -0.13246892
## 10 -2.22420349
## 11 -2.57237878
## 12 -3.01130702
## 13 -2.08651979
## 14 -1.05341383
## 15 -1.15423282
## 16 -1.37015824
## 17 -2.06683383
## 18 -1.06925955
## 19 -0.25222239
## 20  2.59051553
## 21  0.36936918
## 22  2.70693207
## 23  2.04256626
## 24  2.54722885
## 25  2.92583214
## 26  3.01857772
## 27  0.01906754
## 28  1.54521394
## 29  1.86358430
## 30  3.46838961
## 31  1.58434148
## 32  1.25984406
## 33  2.85797639
## 34  1.33024472
## 35  2.13460047
## 36 -0.09443614
## 37  1.98534261
## 38  0.43887312
## 39  0.63270135

Selanjutnya, kita bisa melihat seberapa banyak kesalahan klasifikasi yang dilakukan oleh model yang terbentuk.

#confussion table 
conf <-table(datameasur$Haltica,predict(fit)$class)
conf

##    
##      1  2
##   1 19  0
##   2  1 19

#hit Ratio

#cek akurasinya dengan Hit Ratio   
hit <- sum(diag(prop.table(conf)))
hit

## [1] 0.974359

#APER

aper= 1-hit
aper

## [1] 0.02564103

Penerapan 3

Pertanyaannya

Data yang digunakan pada tutorial ini menunjukkan komposisi kimia dari 48 spesimen tembikar Romawi-Inggris, ditentukan oleh spektrofotometri serapan atom, untuk sembilan oksida (Tubb et all., 1980). Selain komposisi kimia pot, lokasi tungku pembakaran di mana tembikar ditemukan juga dicatat. Untuk data ini, yang menjadi perhatian adalah apakah tembikar memiliki komposisi kimia yang berbeda terkait dengan lokasi tungku pembakaran. Dan komposisi apa yang menyebabkan perbedaan tersebut?

jawab:

langkah-langkah: 1.Menyiapkan data dengan input data di excel atau juga bisa langsung dari R

# library yang digunakan dalam tutorial ini
# jika belum terinstall, silakan diinstall terlebih dahulu kemudian diaktifkan librarynya
library(MASS)
library(stringi)
library(tidyverse)
library(caret)
library(MVN)
library(biotools)

Import Data

pot <- read.csv("pertemuan 13.csv")
knitr::kable(head(pot, 10))

Al2O3	Fe2O3	MgO	CaO	Na2O	K2O	TiO2	MnO	BaO	kiln
18.8	9.52	2.00	0.79	0.40	3.20	1.01	0.077	0.015	1
16.9	7.33	1.65	0.84	0.40	3.05	0.99	0.067	0.018	1
18.2	7.64	1.82	0.77	0.40	3.07	0.98	0.087	0.014	1
16.9	7.29	1.56	0.76	0.40	3.05	1.00	0.063	0.019	1
17.8	7.24	1.83	0.92	0.43	3.12	0.93	0.061	0.019	1
18.8	7.45	2.06	0.87	0.25	3.26	0.98	0.072	0.017	1
16.5	7.05	1.81	1.73	0.33	3.20	0.95	0.066	0.019	1
18.0	7.42	2.06	1.00	0.28	3.37	0.96	0.072	0.017	1
15.8	7.15	1.62	0.71	0.38	3.25	0.93	0.062	0.017	1
14.6	6.87	1.67	0.76	0.33	3.06	0.91	0.055	0.012	1

str(pot)

## 'data.frame':    45 obs. of  10 variables:
##  $ Al2O3: num  18.8 16.9 18.2 16.9 17.8 18.8 16.5 18 15.8 14.6 ...
##  $ Fe2O3: num  9.52 7.33 7.64 7.29 7.24 7.45 7.05 7.42 7.15 6.87 ...
##  $ MgO  : num  2 1.65 1.82 1.56 1.83 2.06 1.81 2.06 1.62 1.67 ...
##  $ CaO  : num  0.79 0.84 0.77 0.76 0.92 0.87 1.73 1 0.71 0.76 ...
##  $ Na2O : num  0.4 0.4 0.4 0.4 0.43 0.25 0.33 0.28 0.38 0.33 ...
##  $ K2O  : num  3.2 3.05 3.07 3.05 3.12 3.26 3.2 3.37 3.25 3.06 ...
##  $ TiO2 : num  1.01 0.99 0.98 1 0.93 0.98 0.95 0.96 0.93 0.91 ...
##  $ MnO  : num  0.077 0.067 0.087 0.063 0.061 0.072 0.066 0.072 0.062 0.055 ...
##  $ BaO  : num  0.015 0.018 0.014 0.019 0.019 0.017 0.019 0.017 0.017 0.012 ...
##  $ kiln : int  1 1 1 1 1 1 1 1 1 1 ...

Ubah Target Data (kiln) menjadi tipe Factor

pot$kiln <- as.factor(pot$kiln)
summary(pot)

##      Al2O3           Fe2O3            MgO             CaO        
##  Min.   :10.10   Min.   :0.920   Min.   :0.530   Min.   :0.0100  
##  1st Qu.:13.80   1st Qu.:5.390   1st Qu.:1.560   1st Qu.:0.1300  
##  Median :16.50   Median :6.920   Median :1.920   Median :0.3000  
##  Mean   :15.71   Mean   :5.756   Mean   :2.488   Mean   :0.5136  
##  3rd Qu.:18.00   3rd Qu.:7.330   3rd Qu.:3.770   3rd Qu.:0.8300  
##  Max.   :20.80   Max.   :9.520   Max.   :7.230   Max.   :1.7300  
##       Na2O             K2O            TiO2             MnO         
##  Min.   :0.0300   Min.   :1.75   Min.   :0.5600   Min.   :0.00100  
##  1st Qu.:0.1000   1st Qu.:2.93   1st Qu.:0.7200   1st Qu.:0.03500  
##  Median :0.2000   Median :3.15   Median :0.9100   Median :0.07200  
##  Mean   :0.2429   Mean   :3.20   Mean   :0.8767   Mean   :0.07051  
##  3rd Qu.:0.3800   3rd Qu.:3.70   3rd Qu.:0.9800   3rd Qu.:0.09400  
##  Max.   :0.8300   Max.   :4.89   Max.   :1.3400   Max.   :0.16300  
##       BaO          kiln  
##  Min.   :0.00900   1:20  
##  1st Qu.:0.01500   2:13  
##  Median :0.01700   3:12  
##  Mean   :0.01651         
##  3rd Qu.:0.01900         
##  Max.   :0.02300

Uji Manova Untuk Melihat Adanya Perbeda

mnv <- manova(cbind(Al2O3,Fe2O3,MgO,CaO,Na2O,K2O,TiO2,MnO,BaO)~kiln, 
data=pot)
summary(mnv)

##           Df Pillai approx F num Df den Df    Pr(>F)    
## kiln       2 1.6836   20.695     18     70 < 2.2e-16 ***
## Residuals 42                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Terlihat bahwa tolak H0 sehingga ada perbedaan yang signifikan Uji Asumsi Analisis Diskriminan

Uji Homogenitas

$$ Hipotesis statistik:\ H_0 : Sigma_1 = sigma_2 = sigma_3\ H_1 : Minimal ada satu yang berbeda sigma_l tidak sama dengan sigma_m\

boxM(pot[, -length(pot)], pot[, length(pot)])

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  pot[, -length(pot)]
## Chi-Sq (approx.) = 346.69, df = 90, p-value < 2.2e-16

Terlihat bahwa tidak memenuhi asumsi homogenitas (tolak H0). Namun pada tutorial kali ini akan dianggap memenuhi asumsi homogenitas

Uji Normalitas

mvn(pot[, -length(pot)])

## $multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 1.382664       0  NO
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling   Al2O3      0.9318  0.0165      NO    
## 2 Anderson-Darling   Fe2O3      3.8481  <0.001      NO    
## 3 Anderson-Darling    MgO       2.3480  <0.001      NO    
## 4 Anderson-Darling    CaO       1.9517  <0.001      NO    
## 5 Anderson-Darling   Na2O       0.9428  0.0155      NO    
## 6 Anderson-Darling    K2O       1.1034  0.0061      NO    
## 7 Anderson-Darling   TiO2       0.8442  0.0274      NO    
## 8 Anderson-Darling    MnO       0.9409  0.0157      NO    
## 9 Anderson-Darling    BaO       0.3557  0.4438      YES   
## 
## $Descriptives
##        n        Mean     Std.Dev Median    Min    Max   25th   75th        Skew
## Al2O3 45 15.70888889 2.703013657 16.500 10.100 20.800 13.800 18.000 -0.37467668
## Fe2O3 45  5.75622222 2.405811419  6.920  0.920  9.520  5.390  7.330 -0.94504825
## MgO   45  2.48844444 1.742110731  1.920  0.530  7.230  1.560  3.770  0.92662287
## CaO   45  0.51355556 0.454278427  0.300  0.010  1.730  0.130  0.830  0.71967637
## Na2O  45  0.24288889 0.178244243  0.200  0.030  0.830  0.100  0.380  0.88614431
## K2O   45  3.20000000 0.852725577  3.150  1.750  4.890  2.930  3.700  0.11120514
## TiO2  45  0.87666667 0.179810506  0.910  0.560  1.340  0.720  0.980  0.37222550
## MnO   45  0.07051111 0.046801137  0.072  0.001  0.163  0.035  0.094  0.07573141
## BaO   45  0.01651111 0.002981932  0.017  0.009  0.023  0.015  0.019 -0.25722580
##          Kurtosis
## Al2O3 -0.99975057
## Fe2O3 -0.58982936
## MgO   -0.17937201
## CaO   -0.36077408
## Na2O   0.64002428
## K2O   -0.82108483
## TiO2   0.01789675
## MnO   -0.71297643
## BaO   -0.12472728

Terlihat bahwa data tidak semuanya berdistribusi normal. Tetapi pada tutorial kali ini akan dianggap memenuhi asumsi normalitas.

Buat model

set.seed(123)
k <- length(levels(pot$kiln))
n <- nrow(pot)

# data training sebanyak
idx <- sample(1:n, 36)
pot.tr <- pot[idx,] #data training
pot.ts <- pot[-idx,] #data testing

# Analisis Diskriminan
fit<-lda(kiln~.,data=pot.tr)
# Output dari analisis diskriminan
fit

## Call:
## lda(kiln ~ ., data = pot.tr)
## 
## Prior probabilities of groups:
##         1         2         3 
## 0.4444444 0.2777778 0.2777778 
## 
## Group means:
##   Al2O3    Fe2O3   MgO   CaO  Na2O      K2O     TiO2       MnO      BaO
## 1 16.70 7.298125 1.820 0.975 0.325 3.073125 0.928125 0.0686875 0.017125
## 2 12.50 6.350000 5.078 0.199 0.256 4.109000 0.697000 0.1248000 0.016000
## 3 16.38 2.343000 1.283 0.087 0.054 2.561000 0.920000 0.0214000 0.014800
## 
## Coefficients of linear discriminants:
##               LD1           LD2
## Al2O3  -0.2401947  1.097991e-01
## Fe2O3  -1.8502033 -4.307396e-01
## MgO     0.1185946 -1.277381e+00
## CaO    -3.1388123  1.087416e+00
## Na2O   -2.9499199 -1.585130e+00
## K2O     2.5061212  1.971959e+00
## TiO2    2.9728578  1.699298e-04
## MnO    11.6804843 -1.398383e+01
## BaO   105.1017528 -4.571878e+01
## 
## Proportion of trace:
##    LD1    LD2 
## 0.8327 0.1673

Scores

zscore <- predict(fit,pot[,-10])
knitr::kable(head(zscore$x, 10))

LD1	LD2
-9.048968	-0.5918104
-4.975976	0.3512298
-5.788318	0.0098015
-4.573428	0.4066486
-5.311813	0.4745376
-4.807566	0.6447346
-6.579490	1.5663482
-4.840231	0.8805399
-3.755891	0.7466634
-4.096007	0.7570749

Korelasi

Korelasi antara var original dengan discriminant score

kor1 <- t(cor(zscore$x[,1],pot[,-10]))
kor2 <- t(cor(zscore$x[,2],pot[,-10]))
tabkor<-data.frame(kor1,kor2)
tabkor

##              kor1        kor2
## Al2O3 -0.33154167  0.66685671
## Fe2O3 -0.78055014 -0.56159595
## MgO    0.20503398 -0.90844076
## CaO   -0.89317269  0.05035898
## Na2O  -0.62309881 -0.40896664
## K2O   -0.02256282 -0.77098113
## TiO2  -0.23688581  0.56947412
## MnO   -0.16076466 -0.87160328
## BaO   -0.20341415 -0.04259347

Lakukan Prediksi Data Testing

pred <- predict(fit,pot.ts[,-10])
pred

## $class
## [1] 1 1 1 1 1 2 2 3 3
## Levels: 1 2 3
## 
## $posterior
##               1            2            3
## 1  1.000000e+00 1.814248e-30 5.279499e-45
## 2  1.000000e+00 2.160024e-17 2.552795e-25
## 6  1.000000e+00 2.970516e-17 2.771709e-24
## 16 1.000000e+00 7.062105e-16 1.228850e-27
## 21 1.000000e+00 5.003768e-16 1.414140e-22
## 24 6.095450e-14 9.999227e-01 7.732419e-05
## 30 9.614902e-18 9.999283e-01 7.169748e-05
## 39 6.815275e-28 4.869188e-09 1.000000e+00
## 45 7.469821e-22 6.335509e-08 9.999999e-01
## 
## $x
##          LD1        LD2
## 1  -9.048968 -0.5918104
## 2  -4.975976  0.3512298
## 6  -4.807566  0.6447346
## 16 -5.200484 -1.1149023
## 21 -4.446521  0.6856987
## 24  2.622321 -1.3610901
## 30  3.510645 -1.7565857
## 39  5.828779  2.3262918
## 45  4.502776  2.4422471

Penilaian akurasi hasil prediksi

# percent correct for each category of k
ct <- table(pot.ts$kiln, pred$class)
ct

##    
##     1 2 3
##   1 4 0 0
##   2 1 2 0
##   3 0 0 2

Terlihat bahwa model berhasil mengklasifikasikan semua data dengan benar

# total percent correct
sum(diag(prop.table(ct)))

## [1] 0.8888889

Plot LDA

lda.data <- cbind(pot.tr, predict(fit)$x)
ggplot(lda.data, aes(LD1, LD2)) +
geom_point(aes(color = as.factor(kiln))) +
theme_minimal()

Penerapan 4

SOAL 3 SOAL UAS APG 2020/2021 NOMOR 3 Ada di File drive ya teman2

#Yang diketahui soal
(xbar.HM<-matrix(c(204.4,406.6),2))

##       [,1]
## [1,] 204.4
## [2,] 406.6

(xbar.M<-matrix(c(130,355),2))

##      [,1]
## [1,]  130
## [2,]  355

s11.HM<-270.4
s22.HM<-117.6
s11.M<-236.6
s22.M<-92.9
r12.HM<-0.5492
r12.M<-0.4925

Pengerjaan:

#Matrik ragam peragaman HM
(s.HM<-matrix(c(s11.HM^2,r12.HM*s11.HM*s22.HM,r12.HM*s11.HM*s22.HM,s22.HM^2),2))

##          [,1]     [,2]
## [1,] 73116.16 17464.03
## [2,] 17464.03 13829.76

#Matrik ragam peragam  M
(s.M<-matrix(c(s11.M^2,r12.M*s11.M*s22.M,r12.M*s11.M*s22.M,s22.M^2),2))

##          [,1]     [,2]
## [1,] 55979.56 10825.22
## [2,] 10825.22  8630.41

#masing- masing jumlah observasinya

n1<-60
n2<-45
(sp<-((n1-1)/((n1-1)+(n2-1)))*s.HM+((n2-1)/((n1-1)+(n2-1)))*s.M)

##          [,1]     [,2]
## [1,] 65795.67 14628.03
## [2,] 14628.03 11608.68

##PERSAMAAN PEMBEDA##
(acap<-t(xbar.HM-xbar.M)%*%solve(sp))

##              [,1]        [,2]
## [1,] 0.0001980263 0.004195419

###Menghitung cut off###
(m<-0.5*(t(xbar.HM-xbar.M)%*%solve(sp)%*%(xbar.HM+xbar.M)))

##          [,1]
## [1,] 1.630725

x0<-matrix(c(150,400),2)
(ycap0<-(acap)%*%x0)

##          [,1]
## [1,] 1.707871

Penerapan 5

EXAMPLE 11.7 JOHNSON HAL 627 (X~normal dan matriks covariance diasumsikan sama)

library(rrcov)
data("salmon")
head(salmon)

##   Gender Freshwater Marine  Origin
## 1      2        108    368 Alaskan
## 2      1        131    355 Alaskan
## 3      1        105    469 Alaskan
## 4      2         86    506 Alaskan
## 5      1         99    402 Alaskan
## 6      2         87    423 Alaskan

salmon<-salmon[ ,-1]
summary(salmon)

##    Freshwater        Marine           Origin  
##  Min.   : 53.0   Min.   :301.0   Alaskan :50  
##  1st Qu.: 99.0   1st Qu.:367.0   Canadian:50  
##  Median :117.5   Median :396.5                
##  Mean   :117.9   Mean   :398.1                
##  3rd Qu.:140.0   3rd Qu.:428.2                
##  Max.   :179.0   Max.   :511.0

##Buat Subset alaska###
(alaska<-subset(salmon, salmon$Origin=="Alaskan"))

##    Freshwater Marine  Origin
## 1         108    368 Alaskan
## 2         131    355 Alaskan
## 3         105    469 Alaskan
## 4          86    506 Alaskan
## 5          99    402 Alaskan
## 6          87    423 Alaskan
## 7          94    440 Alaskan
## 8         117    489 Alaskan
## 9          79    432 Alaskan
## 10         99    403 Alaskan
## 11        114    428 Alaskan
## 12        123    372 Alaskan
## 13        123    372 Alaskan
## 14        109    420 Alaskan
## 15        112    394 Alaskan
## 16        104    407 Alaskan
## 17        111    422 Alaskan
## 18        126    423 Alaskan
## 19        105    434 Alaskan
## 20        119    474 Alaskan
## 21        114    396 Alaskan
## 22        100    470 Alaskan
## 23         84    399 Alaskan
## 24        102    429 Alaskan
## 25        101    469 Alaskan
## 26         85    444 Alaskan
## 27        109    397 Alaskan
## 28        106    442 Alaskan
## 29         82    431 Alaskan
## 30        118    381 Alaskan
## 31        105    388 Alaskan
## 32        121    403 Alaskan
## 33         85    451 Alaskan
## 34         83    453 Alaskan
## 35         53    427 Alaskan
## 36         95    411 Alaskan
## 37         76    442 Alaskan
## 38         95    426 Alaskan
## 39         87    402 Alaskan
## 40         70    397 Alaskan
## 41         84    511 Alaskan
## 42         91    469 Alaskan
## 43         74    451 Alaskan
## 44        101    474 Alaskan
## 45         80    398 Alaskan
## 46         95    433 Alaskan
## 47         92    404 Alaskan
## 48         99    481 Alaskan
## 49         94    491 Alaskan
## 50         87    480 Alaskan

##Buat Subset Canadian###
(canada<-subset(salmon, salmon$Origin=="Canadian"))

##     Freshwater Marine   Origin
## 51         129    420 Canadian
## 52         148    371 Canadian
## 53         179    407 Canadian
## 54         152    381 Canadian
## 55         166    377 Canadian
## 56         124    389 Canadian
## 57         156    419 Canadian
## 58         131    345 Canadian
## 59         140    362 Canadian
## 60         144    345 Canadian
## 61         149    393 Canadian
## 62         108    330 Canadian
## 63         135    355 Canadian
## 64         170    386 Canadian
## 65         152    301 Canadian
## 66         153    397 Canadian
## 67         152    301 Canadian
## 68         136    438 Canadian
## 69         122    306 Canadian
## 70         148    383 Canadian
## 71          90    385 Canadian
## 72         145    337 Canadian
## 73         123    364 Canadian
## 74         145    376 Canadian
## 75         115    354 Canadian
## 76         134    383 Canadian
## 77         117    355 Canadian
## 78         126    345 Canadian
## 79         118    379 Canadian
## 80         120    369 Canadian
## 81         153    403 Canadian
## 82         150    354 Canadian
## 83         154    390 Canadian
## 84         155    349 Canadian
## 85         109    325 Canadian
## 86         117    344 Canadian
## 87         128    400 Canadian
## 88         144    403 Canadian
## 89         163    370 Canadian
## 90         145    355 Canadian
## 91         133    375 Canadian
## 92         128    383 Canadian
## 93         123    349 Canadian
## 94         144    373 Canadian
## 95         140    388 Canadian
## 96         150    339 Canadian
## 97         124    341 Canadian
## 98         125    346 Canadian
## 99         153    352 Canadian
## 100        108    339 Canadian

###Running Plot###
plot(alaska$Freshwater,alaska$Marine,pch=20,col=2,xlim=c(50,200),ylim=c(300,550),main="Plot of Scale Size of Salmon",xlab="Freshwater Scale diameter",ylab="Marine Scale Diameter")
points(canada$Freshwater,canada$Marine,col=3,pch=15)
legend("topright",legend=c("Alaskan Salmon","Canadian Salmon"),pch=c(20,15), col=c(2:3))

###HItung group Means dan S pooled###
n1<-nrow(alaska)
n2<-nrow(canada)

(alaska.means<-apply(alaska[,1:2],2,mean))

## Freshwater     Marine 
##      98.38     429.66

(canada.means<-apply(canada[,1:2],2,mean))

## Freshwater     Marine 
##     137.46     366.62

(alaska_cov<-cov(alaska[,1:2]))

##            Freshwater    Marine
## Freshwater   260.6078 -188.0927
## Marine      -188.0927 1399.0861

(canada_cov<-cov(canada[,1:2]))

##            Freshwater   Marine
## Freshwater   326.0902 133.5049
## Marine       133.5049 893.2608

(sp<-((n1-1)/((n1-1)+(n2-1)))*alaska_cov+((n2-1)/((n1-1)+(n2-1)))*canada_cov)

##            Freshwater     Marine
## Freshwater  293.34898  -27.29388
## Marine      -27.29388 1146.17347

###Menghitung cut off###
(m<-0.5*((alaska.means-canada.means)%*%solve(sp)%*%(alaska.means+canada.means)))

##          [,1]
## [1,] 5.541204

(func<-(alaska.means-canada.means)%*%solve(sp))

##      Freshwater     Marine
## [1,] -0.1283873 0.05194311

species.prediction<-apply(salmon[,1:2],1,function(y) {
  func<-func%*%y
  ifelse(func>m,"Alaska","Canadian")
})

species.prediction

##   [1] "Canadian" "Canadian" "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##   [7] "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Canadian"
##  [13] "Canadian" "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##  [19] "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##  [25] "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Canadian"
##  [31] "Alaska"   "Canadian" "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##  [37] "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##  [43] "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"   "Alaska"  
##  [49] "Alaska"   "Alaska"   "Canadian" "Canadian" "Canadian" "Canadian"
##  [55] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [61] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [67] "Canadian" "Canadian" "Canadian" "Canadian" "Alaska"   "Canadian"
##  [73] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [79] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [85] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [91] "Canadian" "Canadian" "Canadian" "Canadian" "Canadian" "Canadian"
##  [97] "Canadian" "Canadian" "Canadian" "Canadian"

table(salmon$Origin,species.prediction,dnn = c("Actual Group","Predicted Group"))

##             Predicted Group
## Actual Group Alaska Canadian
##     Alaskan      44        6
##     Canadian      1       49

Penerapan 6

The admission officer of a business school has used an “index” of undergraduate grade point average (GPA) and graduate management aptitude test (GMAT) scores to help decide which applicants should be admitted to the school’s graduate programs. Keterangan pada kolom Admit : • Kode 1 = Admit • Kode 2 = Do not Admit • Kode 3 = Borderline a. Tentukan fungsi diskriminan yang terbentuk! b. Hitung t-test untuk setiap variabel! #Anggapannya kita tidak membagi datanya untuk testing dan training

library(readxl)
library(dplyr)

(tugas <- read_excel("PENUGASAN.xlsx"))#masukkan datanya

## # A tibble: 85 x 3
##      GPA  GMAT Admit
##    <dbl> <dbl> <dbl>
##  1  2.96   596     1
##  2  3.14   473     1
##  3  3.22   482     1
##  4  3.29   527     1
##  5  3.69   505     1
##  6  3.46   693     1
##  7  3.03   626     1
##  8  3.19   663     1
##  9  3.63   447     1
## 10  3.59   588     1
## # ... with 75 more rows

str(tugas)

## tibble [85 x 3] (S3: tbl_df/tbl/data.frame)
##  $ GPA  : num [1:85] 2.96 3.14 3.22 3.29 3.69 3.46 3.03 3.19 3.63 3.59 ...
##  $ GMAT : num [1:85] 596 473 482 527 505 693 626 663 447 588 ...
##  $ Admit: num [1:85] 1 1 1 1 1 1 1 1 1 1 ...

tugas$Admit <- as.factor(tugas$Admit)
summary(tugas)

##       GPA             GMAT       Admit 
##  Min.   :2.130   Min.   :313.0   1:31  
##  1st Qu.:2.600   1st Qu.:425.0   2:28  
##  Median :3.010   Median :482.0   3:26  
##  Mean   :2.975   Mean   :488.4         
##  3rd Qu.:3.300   3rd Qu.:538.0         
##  Max.   :3.800   Max.   :693.0

(ad<-filter(tugas[1:2],tugas$Admit=="1"))#Partisi hanya jenis admit

## # A tibble: 31 x 2
##      GPA  GMAT
##    <dbl> <dbl>
##  1  2.96   596
##  2  3.14   473
##  3  3.22   482
##  4  3.29   527
##  5  3.69   505
##  6  3.46   693
##  7  3.03   626
##  8  3.19   663
##  9  3.63   447
## 10  3.59   588
## # ... with 21 more rows

(doadmit<-filter(tugas[1:2],tugas$Admit=="2"))#Partisi hanya jenis do not admit

## # A tibble: 28 x 2
##      GPA  GMAT
##    <dbl> <dbl>
##  1  2.54   446
##  2  2.43   425
##  3  2.2    474
##  4  2.36   531
##  5  2.57   542
##  6  2.35   406
##  7  2.51   412
##  8  2.51   458
##  9  2.36   399
## 10  2.36   482
## # ... with 18 more rows

(boarderline<-filter(tugas[1:2],tugas$Admit=="3"))#Partisi hanya jenis borderline

## # A tibble: 26 x 2
##      GPA  GMAT
##    <dbl> <dbl>
##  1  2.86   494
##  2  2.85   496
##  3  3.14   419
##  4  3.28   371
##  5  2.89   447
##  6  3.15   313
##  7  3.5    402
##  8  2.89   485
##  9  2.8    444
## 10  3.13   416
## # ... with 16 more rows

(n1 <- nrow(ad))#Jumlah Obs dari admit

## [1] 31

(n2 <-nrow(doadmit))#Jumlah Obs dari do nt admid

## [1] 28

(n3 <- nrow(boarderline))#Jumlah Obs dari borderline

## [1] 26

Uji Asumsi Multivariate Normal

library(MVN) # Uji multivariate normal jika skala datanya numerik 

mvn(data = tugas[,c(1:2)],mvnTest = "mardia",multivariatePlot = "adj",showOutliers = F)

## $multivariateNormality
##              Test         Statistic           p value Result
## 1 Mardia Skewness  3.44380148125992 0.486474738991557    YES
## 2 Mardia Kurtosis -1.52773196129016 0.126579101220688    YES
## 3             MVN              <NA>              <NA>    YES
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling    GPA       0.5759    0.1307    YES   
## 2 Anderson-Darling   GMAT       0.4648    0.2483    YES   
## 
## $Descriptives
##       n       Mean    Std.Dev Median    Min   Max  25th  75th        Skew
## GPA  85   2.974588  0.4289954   3.01   2.13   3.8   2.6   3.3 -0.04826749
## GMAT 85 488.447059 81.5223466 482.00 313.00 693.0 425.0 538.0  0.38849916
##        Kurtosis
## GPA  -1.0091148
## GMAT -0.1362411

Dari grafik Chi-Square QQ Plot diatas, dapat dilihat bahwasanya secara umum terbentuk garis linear X = Y, maka dapat dikatakan bahwa data berdistribusi multivariate normal.

Matriks Ragam-peragam antar kategori sama

library(biotools)
#melakukan uji statistik Box’s M di R dapat menggunakan fungsi boxM.
boxM(tugas[,c(1:2)],tugas$Admit)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  tugas[, c(1:2)]
## Chi-Sq (approx.) = 16.074, df = 6, p-value = 0.01336

(wilk <- manova(formula = cbind(tugas$GPA, tugas$GMAT)~tugas$Admit))

## Call:
##    manova(formula = cbind(tugas$GPA, tugas$GMAT) ~ tugas$Admit)
## 
## Terms:
##                 tugas$Admit Residuals
## resp 1                12.50      2.96
## resp 2             258471.1  299783.9
## Deg. of Freedom           2        82
## 
## Residual standard errors: 0.1899156 60.46405
## Estimated effects may be unbalanced

summary(object = wilk, test = "Wilks")

##             Df   Wilks approx F num Df den Df    Pr(>F)    
## tugas$Admit  2 0.12638   73.426      4    162 < 2.2e-16 ***
## Residuals   82                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

utput diatas menunjukkan bahwa dengan tingkat signifikansi 5%, didapat keputusan untuk menolak hipotesis nol atau dengan kata lain terdapat perbedaan rata-rata (nilai variabel prediktor) antar kategori

Selain itu, untuk memeriksa korelasi dan multikolinieritas antar variabel independen dapat menjalankan perintah berikut

#Correlation dan multikolinearity check
kor<-cor(tugas[,1:2])
kor

##            GPA      GMAT
## GPA  1.0000000 0.4606332
## GMAT 0.4606332 1.0000000

Kemudian dilakukan canonical discriminant analysis (analisis diskriminan kanonik) untuk melihat hubungan antara variabel diskriminan dengan variabel independen secara multivariat

library(candisc)
#Canonical Correlation (CC)
cc<-candisc(wilk)
cc

## 
## Canonical Discriminant Analysis for tugas$Admit:
## 
##    CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.84953    5.64604     5.4554 96.7342     96.734
## 2 0.16010    0.19061     5.4554  3.2658    100.000
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F numDF denDF   Pr(> F)    
## 1      0.12638   73.426     4   162 < 2.2e-16 ***
## 2      0.83990   15.630     1    82 0.0001626 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#Eigenvalue untuk fungsi deskriminan yang pertama menunjukan ada tidaknya multikolinearitas antar variabel. Didapatkan nilai Eigenvalue = 5,6404.Dikarenakan nilai eigenvalue tidak mendekati 0, maka tidak terjadi multikolinearitas antar variabel independen.

untuk fungsi deskriminan yang kedua menunjukan ada multikolinearitas antar variabel. Didapatkan nilai Eigenvalue = 0.19061.Dikarenakan nilai eigenvalue mendekati 0, maka tidak terjadi multikolinearitas antar variabel independen.

#Canonical Corellation Canonical Correlation pertama (CC)=0.8495, artinya bahwa hubungan antara variabel diskriminan dengan variabel independen secara multivariat sebesar 0.8495 atau besarnya (CC)²= (0.8495)²= 0,721. #Jadi dapat disimpulkan bahwa 72,1% variasi antara ketiga jenis Admid dapat dijelaskan oleh variabel GPA dan GMAT ).

Memulai Analisis Diskriminan

Analisis Diskriminan Dengan Menggunakan Package membentuk fungsi diskriminan

Linear Discriminant An. with Jacknifed Prediction

library(MASS)
fit <- lda(Admit~., data =  tugas)
fit

## Call:
## lda(Admit ~ ., data = tugas)
## 
## Prior probabilities of groups:
##         1         2         3 
## 0.3647059 0.3294118 0.3058824 
## 
## Group means:
##        GPA     GMAT
## 1 3.403871 561.2258
## 2 2.482500 447.0714
## 3 2.992692 446.2308
## 
## Coefficients of linear discriminants:
##               LD1         LD2
## GPA  -5.008766354  1.87668220
## GMAT -0.008568593 -0.01445106
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9673 0.0327

berdasarkan hasil output, maka persamaan yang terbentuk adalah terbentuk satu persamaan : \[ ycap1 = -5.008766354GPA-0.008568593GMAT\\ ycap2 = 1.87668220GPA-0.01445106GMAT\\ \] Kemudian akan dilakukan prediksi untuk menentukan kelas peluang anggota setiap kelas dan diskriminan linear, serta membandingkan data asli dengan data prediksi. dapat dijalankan dengan syntax sebagai berikut :

#membuat prediksi 
predik <-predict(fit,tugas)
predik

## $class
##  [1] 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 2 2 2 2 2 2 2
## [39] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 1 3 3 3 3 3 3 3 3 2 3
## [77] 3 3 3 3 3 3 3 3 3
## Levels: 1 2 3
## 
## $posterior
##               1            2            3
## 1  6.274441e-01 4.805886e-03 3.677500e-01
## 2  1.400290e-01 2.125682e-03 8.578453e-01
## 3  4.070245e-01 4.261541e-04 5.925494e-01
## 4  9.111531e-01 1.618774e-05 8.883071e-02
## 5  9.989964e-01 6.361722e-10 1.003605e-03
## 6  9.999849e-01 6.281005e-11 1.506768e-05
## 7  9.339968e-01 2.448575e-04 6.575834e-02
## 8  9.981064e-01 5.117961e-07 1.893124e-03
## 9  9.787312e-01 5.091012e-08 2.126871e-02
## 10 9.998460e-01 2.192548e-10 1.540215e-04
## 11 9.794875e-01 2.437537e-06 2.051009e-02
## 12 9.921127e-01 2.357944e-07 7.887049e-03
## 13 9.990257e-01 5.840399e-09 9.743436e-04
## 14 9.999896e-01 9.039119e-13 1.038230e-05
## 15 9.999794e-01 1.156024e-10 2.055507e-05
## 16 9.929468e-01 7.984816e-08 7.053134e-03
## 17 9.968155e-01 3.457110e-08 3.184444e-03
## 18 9.463457e-01 4.304984e-06 5.365001e-02
## 19 9.867345e-01 4.961726e-07 1.326496e-02
## 20 8.845399e-01 2.511215e-05 1.154350e-01
## 21 7.954080e-01 1.168806e-04 2.044752e-01
## 22 9.995508e-01 8.926617e-10 4.492087e-04
## 23 9.872747e-01 9.612200e-07 1.272435e-02
## 24 5.205375e-01 3.722376e-05 4.794253e-01
## 25 9.986259e-01 3.661300e-08 1.374081e-03
## 26 9.992954e-01 6.809424e-08 7.045741e-04
## 27 9.999407e-01 6.195388e-11 5.934438e-05
## 28 9.906305e-01 4.139368e-07 9.369097e-03
## 29 9.998791e-01 1.359825e-11 1.209210e-04
## 30 9.999984e-01 1.211841e-13 1.598827e-06
## 31 3.342864e-01 4.019391e-04 6.653117e-01
## 32 8.119618e-07 9.507696e-01 4.922959e-02
## 33 1.430843e-08 9.912438e-01 8.756170e-03
## 34 2.195048e-10 9.995474e-01 4.526164e-04
## 35 2.840709e-07 9.927593e-01 7.240407e-03
## 36 1.514891e-04 8.546137e-01 1.452348e-01
## 37 6.221596e-10 9.976347e-01 2.365289e-03
## 38 7.291784e-08 9.750166e-01 2.498334e-02
## 39 6.203497e-07 9.645788e-01 3.542058e-02
## 40 5.941739e-10 9.974098e-01 2.590231e-03
## 41 2.877452e-08 9.950535e-01 4.946453e-03
## 42 5.977416e-06 8.044211e-01 1.955729e-01
## 43 7.521694e-06 7.630422e-01 2.369502e-01
## 44 8.847701e-06 9.591041e-01 4.088704e-02
## 45 1.665217e-06 9.742825e-01 2.571580e-02
## 46 1.246193e-04 7.674623e-01 2.324131e-01
## 47 2.961470e-10 9.949313e-01 5.068668e-03
## 48 1.395205e-12 9.999027e-01 9.734658e-05
## 49 6.376893e-08 9.907758e-01 9.224144e-03
## 50 7.452464e-05 8.903745e-01 1.095510e-01
## 51 2.068197e-08 9.971380e-01 2.862023e-03
## 52 1.622653e-07 9.892303e-01 1.076950e-02
## 53 8.697548e-12 9.997610e-01 2.390311e-04
## 54 1.168557e-11 9.987821e-01 1.217885e-03
## 55 3.753738e-07 9.236190e-01 7.638061e-02
## 56 4.707168e-05 8.978003e-01 1.021526e-01
## 57 1.008537e-05 6.689518e-01 3.310382e-01
## 58 6.544192e-05 2.589159e-01 7.410186e-01
## 59 1.680421e-04 1.413231e-01 8.585089e-01
## 60 7.284975e-03 1.101171e-01 8.825979e-01
## 61 6.774234e-03 1.235863e-01 8.696395e-01
## 62 1.945808e-02 3.692965e-03 9.768489e-01
## 63 2.012813e-02 6.995316e-04 9.791723e-01
## 64 1.777925e-03 1.040385e-01 8.941836e-01
## 65 3.644312e-04 7.429134e-03 9.922064e-01
## 66 5.770402e-01 9.584604e-06 4.229502e-01
## 67 7.977472e-03 7.891988e-02 9.131027e-01
## 68 3.611077e-04 3.066731e-01 6.929658e-01
## 69 1.516921e-02 4.390308e-03 9.804405e-01
## 70 2.478183e-02 1.607342e-02 9.591448e-01
## 71 2.007239e-03 2.627083e-01 7.352845e-01
## 72 9.410361e-04 1.163918e-01 8.826672e-01
## 73 2.301748e-03 8.035770e-02 9.173406e-01
## 74 9.842772e-03 2.894799e-01 7.006773e-01
## 75 2.434555e-04 5.063012e-01 4.934553e-01
## 76 7.746504e-02 3.297301e-03 9.192377e-01
## 77 1.944494e-02 7.487706e-03 9.730674e-01
## 78 4.370318e-03 1.837502e-02 9.772547e-01
## 79 8.900878e-02 1.293484e-02 8.980564e-01
## 80 9.140399e-03 1.580157e-02 9.750580e-01
## 81 2.641752e-03 1.611669e-02 9.812416e-01
## 82 4.040008e-03 1.353891e-01 8.605709e-01
## 83 1.241687e-02 1.869183e-02 9.688913e-01
## 84 3.597029e-03 1.910875e-02 9.772942e-01
## 85 1.426707e-02 1.280087e-02 9.729321e-01
## 
## $x
##            LD1         LD2
## 1  -0.84850831 -1.58163149
## 2  -0.69614932  0.53365169
## 3  -1.17396797  0.55372673
## 4  -1.91016829  0.03479678
## 5  -3.72516579  1.10339298
## 6  -4.18404500 -2.04504321
## 7  -1.45617974 -1.88379553
## 8  -2.57462030 -2.11821560
## 9  -2.92766142  1.82895353
## 10 -3.93548237 -0.28371322
## 11 -2.26872530 -0.46667456
## 12 -2.68391601 -0.13449574
## 13 -3.34759591 -0.22139766
## 14 -4.91285376  0.02950322
## 15 -4.07530108 -2.06812579
## 16 -2.87040249  0.37691534
## 17 -3.02596106  0.01132308
## 18 -2.15071413  0.24855513
## 19 -2.54814242 -0.00875196
## 20 -1.82580626  0.07383420
## 21 -1.53517277 -0.15869098
## 22 -3.67974848  0.04434540
## 23 -2.43612548 -0.43927621
## 24 -1.63854767  1.62853358
## 25 -3.02930752 -0.92348450
## 26 -2.93380254 -2.00129891
## 27 -4.16551049 -0.56841866
## 28 -2.58506458 -0.27750256
## 29 -4.41322758  1.07861106
## 30 -5.28395104 -0.80283872
## 31 -1.14561440  0.80802627
## 32  2.54046250 -0.20217901
## 33  3.27136725 -0.10514180
## 34  4.00352246 -1.24488064
## 35  2.71371004 -1.76832191
## 36  1.56761459 -1.53318031
## 37  3.83487182  0.01929377
## 38  2.98205765  0.23285656
## 39  2.58790237 -0.43189220
## 40  3.84476431  0.13921801
## 41  3.13357110 -1.06021997
## 42  2.16219395  0.39875041
## 43  2.11343018  0.52299042
## 44  2.09552089 -1.57202217
## 45  2.40134245 -1.26273037
## 46  1.59269514 -0.87143910
## 47  3.98388436  1.19976937
## 48  4.91966324 -0.42247844
## 49  2.99452449 -0.77852208
## 50  1.70206428 -1.51290971
## 51  3.18693178 -1.48642846
## 52  2.82315263 -1.06754328
## 53  4.59343148 -0.35323068
## 54  4.56320223  1.24763387
## 55  2.68550335  0.66187704
## 56  1.78775021 -1.36839911
## 57  2.04160842  0.81482360
## 58  1.54470347  1.31891137
## 59  1.26855938  1.36939230
## 60  0.52636480 -0.29529159
## 61  0.55931528 -0.34296053
## 62 -0.23344530  1.31400893
## 63 -0.52338013  2.27039532
## 64  0.77882568  0.44020870
## 65  0.62473789  2.86458811
## 66 -1.89093511  2.23528254
## 67  0.45321915 -0.10893158
## 68  1.25532043  0.31466048
## 69 -0.15765186  1.33859529
## 70 -0.02787251  0.31858512
## 71  0.91125282 -0.36885510
## 72  0.91592317  0.67142566
## 73  0.68721895  0.49219340
## 74  0.63176227 -1.25318175
## 75  1.40885644 -0.14908165
## 76 -0.51028807  0.64062865
## 77 -0.11285977  0.89793574
## 78  0.31751900  1.10757389
## 79 -0.30339138 -0.24932198
## 80  0.15471573  0.83300375
## 81  0.38871553  1.43412873
## 82  0.67070699 -0.15509675
## 83  0.12636216  0.57870420
## 84  0.36036196  1.17982919
## 85  0.03607932  0.73616209

Selanjutnya, kita bisa melihat seberapa banyak kesalahan klasifikasi yang dilakukan oleh model yang terbentuk.

#confussion table 
conf <-table(tugas$Admit,predict(fit)$class)
conf

##    
##      1  2  3
##   1 28  0  3
##   2  0 26  2
##   3  1  1 24

#hit Ratio

#cek akurasinya dengan Hit Ratio   
hit <- sum(diag(prop.table(conf)))
hit

## [1] 0.9176471

#APER

aper= 1-hit
aper

## [1] 0.08235294

Penerapan 7 Diskriminan Fisher 3 Populasi

x1<-matrix(c(-2,0,-1,5,3,1),3,2)
x1

##      [,1] [,2]
## [1,]   -2    5
## [2,]    0    3
## [3,]   -1    1

x2<-matrix(c(0,2,1,6,4,2),3,2)
x2

##      [,1] [,2]
## [1,]    0    6
## [2,]    2    4
## [3,]    1    2

x3<-matrix(c(1,0,-1,-2,0,-4),3,2)
x3

##      [,1] [,2]
## [1,]    1   -2
## [2,]    0    0
## [3,]   -1   -4

xbar1<-matrix(c(colMeans(x1)),2,1)
xbar1

##      [,1]
## [1,]   -1
## [2,]    3

xbar2<-matrix(c(colMeans(x2)),2,1)
xbar2

##      [,1]
## [1,]    1
## [2,]    4

xbar3<-matrix(c(colMeans(x3)),2,1)
xbar3

##      [,1]
## [1,]    0
## [2,]   -2

xbar<-(xbar1+xbar2+xbar3)/3
xbar

##          [,1]
## [1,] 0.000000
## [2,] 1.666667

B<-(xbar1-xbar)%*%t(xbar1-xbar)+(xbar2-xbar)%*%t(xbar2-xbar)+(xbar3-xbar)%*%t(xbar3-xbar)
B

##      [,1]     [,2]
## [1,]    2  1.00000
## [2,]    1 20.66667

W<-2*(var(x1)+var(x2)+var(x3))
W

##      [,1] [,2]
## [1,]    6   -2
## [2,]   -2   24

E<-solve(W)%*%B
E

##            [,1]      [,2]
## [1,] 0.35714286 0.4666667
## [2,] 0.07142857 0.9000000

Spl = W/(3+3+3-3)
eigen(E)

## eigen() decomposition
## $values
## [1] 0.9556904 0.3014525
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.6148676 -0.9929546
## [2,] -0.7886303  0.1184957

I = diag(2)
g.1 = E-eigen(E)$values[1]*I
g.2 = E-eigen(E)$values[2]*I
ve1.g.1 = eigen(g.1)$vectors[,1]
ve2.g.1 = eigen(g.1)$vectors[,2]
ve1.g.2 = eigen(g.2)$vectors[,1]
ve2.g.2 = eigen(g.2)$vectors[,2]
# Normalisasi Vektor
library(ppls)
a = normalize.vector(ve1.g.1)
b = normalize.vector(ve2.g.1)
c = normalize.vector(ve1.g.2)
d = normalize.vector(ve2.g.2)
a

## [1] -0.9929546  0.1184957

## [1] -0.6148676 -0.7886303

## [1] -0.6148676 -0.7886303

## [1] -0.9929546  0.1184957

Deteksi Kuadratik Separate Kak Nadhifan STIS 54

Soal PPT Pertemuan 13

SOAL 1

x1 = matrix(c(3,2,4,7,4,7),3,2)
x2 = matrix(c(6,5,4,9,7,8),3,2)
xbar.1 = matrix(c(colMeans(x1)),2,1)
xbar.2 = matrix(c(colMeans(x2)),2,1)
xbar.1

##      [,1]
## [1,]    3
## [2,]    6

xbar.2

##      [,1]
## [1,]    5
## [2,]    8

Spl = matrix(c(1,1,1,2),2,2)

# Jawaban A
a=t(xbar.1-xbar.2)%*%solve(Spl)  #fungsi diskrim linier = -2x1


# Jawaban B
x.o = matrix(c(2,7),2,1)
ybar.1 = a%*%xbar.1
ybar.1

##      [,1]
## [1,]   -6

ybar.2 = a%*%xbar.2
ybar.2

##      [,1]
## [1,]  -10

# cutting score
m.kep = 0.5*(ybar.1+ybar.2)
m.kep

##      [,1]
## [1,]   -8

# Cek
a%*%x.o  # > m.kep --> masuk pi.1

##      [,1]
## [1,]   -4

Nomor 2

# Diolah secara ghaib diperoleh :
xbar.1 = matrix(c(109.475,20.267),2,1)
xbar.1

##         [,1]
## [1,] 109.475
## [2,]  20.267

xbar.2 = matrix(c(87.4,17.633),2,1)
xbar.2

##        [,1]
## [1,] 87.400
## [2,] 17.633

s1 = matrix(c(352.644,-11.818,-11.818,4.082),2,2)
s2 = matrix(c(200.705, -2.589,-2.589,4.464),2,2)
s1

##         [,1]    [,2]
## [1,] 352.644 -11.818
## [2,] -11.818   4.082

s2

##         [,1]   [,2]
## [1,] 200.705 -2.589
## [2,]  -2.589  4.464

Spl = matrix(c(276.675,-7.204,-7.204,4.273),2,2)
Spl

##         [,1]   [,2]
## [1,] 276.675 -7.204
## [2,]  -7.204  4.273

# Bagian A
a=t(xbar.1-xbar.2)%*%solve(Spl)
a  #fungsi klasifikasi linier

##           [,1]      [,2]
## [1,] 0.1002374 0.7854225

#Cutting score
m.kep = 0.5*(a%*%xbar.1+a%*%xbar.2)
m.kep

##          [,1]
## [1,] 24.75088

Soal UAS (Pertemuan 14)

Dokumentasi Praktikum

Dokumentasi MVN

Terkadang Ambisi Tinggi dapat Menjatuhkan ke Jurang Terdalam, Lakukan versi terbaik dan Tetaplah Bersyukur Matahari yang terbit di pertemuan 1 dan berjalan di setiap pertemuan yang hangat, kini akan tenggelam dalam senja di sebuah perpisahan yang hangat. Waktu yang menjadi saksi atas bertahanya kerja keras kita selama ini. Semoga di malam hari, bisa tidur nyenyak tanpa ada rasa dendam. Semoga sukses ke depannya !!!