1.1 Det, Inv, tr, t, diag

x <- matrix(c(3,2,-1,
              1,6,3,
              2,-4,0),nrow=3, byrow=T);x

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

t(x) #Transpose

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

det(x) #Determinan

## [1] 64

install_load('MASS')
fractions(solve(x)) #Invers

##      [,1]  [,2]  [,3] 
## [1,]  3/16  1/16  3/16
## [2,]  3/32  1/32 -5/32
## [3,]  -1/4   1/4   1/4

install_load('psych')
tr(x) #Teras

## [1] 9

diag(x) #mengambil diagonal matriks

## [1] 3 6 0

diag(c(1,2,3)) #Membuat matriks diagonal

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

#Baris 3 = baris 1
x <- matrix(c(3,2,-1,
              6,4,-2,
              3,2,-1),nrow=3, byrow=T);x

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

det(x)

## [1] 0

fractions(solve(x)) #Invers

## Error in solve.default(x): Lapack routine dgesv: system is exactly singular: U[2,2] = 0

#Baris 2 = baris 1 dikali 2
x <- matrix(c(3,2,-1,
              6,4,-2,
              2,-4,0),nrow=3, byrow=T);x

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

det(x)

## [1] 0

fractions(solve(x)) #Invers

## Error in solve.default(x): Lapack routine dgesv: system is exactly singular: U[3,3] = 0

#Baris 3 = baris 1 ditambah baris 2
x <- matrix(c(3,2,-1,
              1,1,1,
              4,3,0),nrow=3, byrow=T);x

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

det(x)

## [1] 0

fractions(solve(x)) #Invers

## Error in solve.default(x): Lapack routine dgesv: system is exactly singular: U[3,3] = 0

#Kolom 1 = kolom 2 dikurang kolom 3
x <- matrix(c(-3,-2,1,
              0, 1,1,
              3, 3,0),nrow=3, byrow=T);x

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

det(x)

## [1] 0

fractions(solve(x)) #Invers

## Error in solve.default(x): Lapack routine dgesv: system is exactly singular: U[3,3] = 0

1.2 Ortogonal

x <- matrix(c( 1/3, -2/3,  2/3,
               2/3,  2/3,  1/3,
              -2/3,  1/3,  2/3),nrow=3, byrow=T);fractions(x)

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

y <- matrix(c( 1/3, -2/3,  2/3,
               2/3,  2/3,  1/3,
              -2/3,  1/3,  5/3),nrow=3, byrow=T);fractions(y)

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

fractions(solve(x))

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

fractions(t(x))

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

t(x) %*% x

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

x %*% t(x)

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

orto <- function(x){
  if( all(t(x) %*% x == diag(nrow(x))) ){
    cat("Matriks Ortgonal")
  } else cat('Bukan Matriks Ortogonal')
}
orto(x)

## Matriks Ortgonal

a <- c(5, 9, 5); a

## [1] 5 9 5

b <- c(-1, 0, 1); b

## [1] -1  0  1

c <- c(-1, 1, 1); c

## [1] -1  1  1

t(a) %*% b; t(b) %*% a

##      [,1]
## [1,]    0

##      [,1]
## [1,]    0

orto <- function(x, y=NA){
  if(is.matrix(x)){
    if( all(t(x) %*% x == diag(nrow(x))) ){
      cat("\nMatriks x Ortgonal\n\nkarena, t(x)\n")
      print(fractions(t(x)))
      
      cat("\nDikali dengan x\n")
      print(fractions(x))
      
      cat("\nSama dengan I\n")
      print(t(x) %*% x)
      
    } else {
      cat('\nBukan Matriks Ortogonal\n\nkarena, t(x)\n')
      print(fractions(t(x)))
      
      cat("\nDikali dengan x\n")
      print(fractions(x))
      
      cat("\nTidak Sama dengan I, melainkan:\n")
      print(fractions(t(x) %*% x))
    }
  } else{
      if( (t(x) %*% y | t(y) %*% x) == 0 ){
      cat("\nVektor Ortgonal\n\nkarena, t(x)\n")
      print(x)
      
      cat("\nDikali dengan y\n")
      print(y)
      
      cat("\nSama dengan\n")
      print(t(x) %*% y)
      } else {
        cat('\nBukan Vektor Ortogonal\n\nkarena, t(x)\n')
        print(x)
        
        cat("\nDikali dengan y\n")
        print(y)
        
        cat("\nTidak Sama dengan nol, melainkan\n")
        print(t(x) %*% y)
    }
  }
}
orto(y)

## 
## Bukan Matriks Ortogonal
## 
## karena, t(x)
##      [,1] [,2] [,3]
## [1,]  1/3  2/3 -2/3
## [2,] -2/3  2/3  1/3
## [3,]  2/3  1/3  5/3
## 
## Dikali dengan x
##      [,1] [,2] [,3]
## [1,]  1/3 -2/3  2/3
## [2,]  2/3  2/3  1/3
## [3,] -2/3  1/3  5/3
## 
## Tidak Sama dengan I, melainkan:
##      [,1] [,2] [,3]
## [1,]    1    0 -2/3
## [2,]    0    1  1/3
## [3,] -2/3  1/3 10/3

orto(a, c)

## 
## Bukan Vektor Ortogonal
## 
## karena, t(x)
## [1] 5 9 5
## 
## Dikali dengan y
## [1] -1  1  1
## 
## Tidak Sama dengan nol, melainkan
##      [,1]
## [1,]    9

1.3 Rank

qr(x)$rank #pake ini

## [1] 3

rank <- function(X){ #Kalo pake function
  if (qr(X)$rank == ncol(X)){ 
      cat("Matriks Model Penuh, dengan rank :", qr(X)$rank)
  } else {
      cat("Matriks BUKAN Model Penuh, dengan rank :", qr(X)$rank) 
  }
}
rank(x)

## Matriks Model Penuh, dengan rank : 3

1.4 Eigen (Akar Ciri)

A <- matrix(c(2, 0, 1,
              0, 1, 1,
              1, 1, 1), nrow=3, byrow=T)
eigen(A)

## eigen() decomposition
## $values
## [1]  2.8019377  1.4450419 -0.2469796
## 
## $vectors
##           [,1]       [,2]       [,3]
## [1,] 0.7369762  0.5910090 -0.3279853
## [2,] 0.3279853 -0.7369762 -0.5910090
## [3,] 0.5910090 -0.3279853  0.7369762

x <- eigen(A)$values
t(x) %*% A %*% x

##          [,1]
## [1,] 15.75302

sum(eigen(A)$values)

## [1] 4

tr(A)

## [1] 4

prod(eigen(A)$values)

## [1] -1

det(A)

## [1] -1

A <- diag(c(2, 4, 3)); A

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

x <- eigen(A)$values; x

## [1] 4 3 2

v <- eigen(A)$vectors; fractions(v)

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

1.5 Leading & Free Variabel

# Misalnya kita memiliki sistem persamaan linear berikut:
# 2x + 3y + z = 5
# 4x + 7y + 2z = 10
# 6x + 10y + 3z = 15

# Membuat matriks augmented dari sistem persamaan di atas
x <- matrix(c(2, 3, 1, 5, 
              4, 7, 2, 10, 
              6, 10, 3, 15), nrow = 3, byrow = TRUE)

lead.free <- function(x){
  l.var <- NA; f.var <- NA
  #menghitung bentuk eselon baris dari matriks augmented
  install_load('pracma')
  es <- rref(x)[,-ncol(x)]
  # Mencari leading variable dan free variable
  for (i in 1:nrow(x)) {
    pivot <- which(es[i, ] == 1)[1]  #Mencari indeks kolom pertama yang bernilai 1
    
    if (!is.na(pivot)) {
      l.var <- na.omit(c(l.var, pivot)) 
    } else {
      f.var <- na.omit(c(1:ncol(es))[-l.var])
    }
  }
  
  print("Bentuk eselon:")
  print(fractions(rref(x)))
  print("Leading Variable:")
  print(l.var)
  
  print("Free Variable:")
  print(f.var)
  if(!is.na(all(f.var))) cat(" Karena Memiliki Free Variable,\n",
                        'Maka Persamaan ini memiliki Solusi tak terbatas')
}
lead.free(x)

## [1] "Bentuk eselon:"
##      [,1] [,2] [,3] [,4]
## [1,]   1    0  1/2  5/2 
## [2,]   0    1    0    0 
## [3,]   0    0    0    0 
## [1] "Leading Variable:"
## [1] 1 2
## [1] "Free Variable:"
## [1] 3
##  Karena Memiliki Free Variable,
##  Maka Persamaan ini memiliki Solusi tak terbatas

#Contoh lain 
x <- matrix(c(1, 1, 1, 1, 6, 
              1, 0, 1, 1, 4, 
              1, 0, 1, 1, 2), nrow = 3, byrow = TRUE)
lead.free(x)

## [1] "Bentuk eselon:"
##      [,1] [,2] [,3] [,4] [,5]
## [1,] 1    0    1    1    0   
## [2,] 0    1    0    0    0   
## [3,] 0    0    0    0    1   
## [1] "Leading Variable:"
## [1] 1 2
## [1] "Free Variable:"
## [1] 3 4
##  Karena Memiliki Free Variable,
##  Maka Persamaan ini memiliki Solusi tak terbatas

#Contoh lain 
x <- matrix(c(1, 1, -3,  5, 8, 10, 
              0, 0,  1, -3, 4, 15, 
              0, 0,  0,  1, 3, 4,
              0, 0,  0,  0, 0, 0), nrow = 4, byrow = TRUE)
lead.free(x) 

## [1] "Bentuk eselon:"
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  1    1    0    0   32   71  
## [2,]  0    0    1    0   13   27  
## [3,]  0    0    0    1    3    4  
## [4,]  0    0    0    0    0    0  
## [1] "Leading Variable:"
## [1] 1 3 4
## [1] "Free Variable:"
## [1] 2 5
##  Karena Memiliki Free Variable,
##  Maka Persamaan ini memiliki Solusi tak terbatas

#Contoh lain 
x <- matrix(c(1, 1, -3,  5, 
              0, 1, -1, 15,
              0, 0,  1, 12), nrow = 3, byrow = TRUE)
lead.free(x) 

## [1] "Bentuk eselon:"
##      [,1] [,2] [,3] [,4]
## [1,]  1    0    0   14  
## [2,]  0    1    0   27  
## [3,]  0    0    1   12  
## [1] "Leading Variable:"
## [1] 1 2 3
## [1] "Free Variable:"
## [1] NA

1.6 Solusi Pertamaan Linier

Referensi : link

Misal persamaannya :

\(2x + 3y = 84\)

\(x - 2y = 2\)

A <- matrix(c(2, 3, 
              4, -2), nrow = 2, byrow = TRUE)
b <- c(8, 2)
solve(A, b)

## [1] 1.375 1.750

A <- matrix(c(2, 3, 
              4, -2), nrow = 2, byrow = TRUE)
b <- c(8, 2)
solusi <- function(A, b=NA){
  if(missing(b)) {
    b <- A[,ncol(A)] #Jika tidak ada b, maka ambil dari A
    A <- A[,-ncol(A)]
  }
  n <- ncol(A)
  install_load('matlib')
  cat('Persamaan :\n')
  showEqn(A, b)
  
  cat('\nRank Matriks :',R(A),
      '\nRank Gabung  :', R(cbind(A,b)) )
  
  plotEqn(A,b)
  
  if(all.equal( R(A), R(cbind(A,b)) )){
    if(all.equal( R(A), R(cbind(A,b)), n)){
      cat('\nPersamaan Konsisten dan Solusi Unik\n')
    } else{
      cat('\nPersamaan Konsisten dan Solusi Tak Terbatas\n')
    }
    Solve(A, b, fractions = TRUE)
  } else{
      cat('\nPersamaan Tidak Konsisten (Tidak memiliki Solusi)\n')
    }
}
solusi(A,b)

## Persamaan :
## 2*x1 + 3*x2  =  8 
## 4*x1 - 2*x2  =  2 
## 
## Rank Matriks : 2 
## Rank Gabung  : 22*x[1] + 3*x[2]  =  8 
## 4*x[1] - 2*x[2]  =  2

## 
## Persamaan Konsisten dan Solusi Unik
## x1    =  11/8 
##   x2  =   7/4

#Contoh Jika A dan b digabung : 
A <- matrix(c(2, 3, 8,
              4, -2, 2), nrow = 2, byrow = TRUE)
solusi(A)

## Persamaan :
## 2*x1 + 3*x2  =  8 
## 4*x1 - 2*x2  =  2 
## 
## Rank Matriks : 2 
## Rank Gabung  : 22*x[1] + 3*x[2]  =  8 
## 4*x[1] - 2*x[2]  =  2

## 
## Persamaan Konsisten dan Solusi Unik
## x1    =  11/8 
##   x2  =   7/4

1.7 Matriks Korelasi dan Matriks Kovarian

A <- matrix(c(25, -2, 4,
              -2,  4, 1,
               4,  1, 9), nrow=3, byrow= T)
#Varians
fractions(var(A))

##      [,1]  [,2]  [,3] 
## [1,]   201 -81/2     4
## [2,] -81/2     9  -9/2
## [3,]     4  -9/2  49/3

#Sama dengan
satu <- rep(1,nrow(A))

fractions(1/(nrow(A)-1) * t(A) %*% 
  (diag(satu) - 1/nrow(A) * satu %*% t(satu)) %*% A)

##      [,1]  [,2]  [,3] 
## [1,]   201 -81/2     4
## [2,] -81/2     9  -9/2
## [3,]     4  -9/2  49/3

#Sama dengan
#Kovarians
fractions(cov(A))

##      [,1]  [,2]  [,3] 
## [1,]   201 -81/2     4
## [2,] -81/2     9  -9/2
## [3,]     4  -9/2  49/3

D <- diag(sqrt(diag(A))); D

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

fractions(cov2cor(A))

##      [,1] [,2] [,3]
## [1,]    1 -1/5 4/15
## [2,] -1/5    1  1/6
## [3,] 4/15  1/6    1

#Sama dengan
fractions( solve(D) %*% A %*% solve(D) )

##      [,1] [,2] [,3]
## [1,]    1 -1/5 4/15
## [2,] -1/5    1  1/6
## [3,] 4/15  1/6    1

rowMeans(A)

## [1] 9.000000 1.000000 4.666667

#sama dengan
1/nrow(A) * t(A) %*% satu

##          [,1]
## [1,] 9.000000
## [2,] 1.000000
## [3,] 4.666667

1.8 Euclid

a <- c(2, 1, 2)
fractions(a / sqrt(t(a) %*% a))

## [1] 2/3 1/3 2/3

a <- c(5, 3, 2); a

## [1] 5 3 2

b <- c(6, 1, 4); b

## [1] 6 1 4

sqrt(t(a - b) %*% (a - b))

##      [,1]
## [1,]    3

1.9 Perkalian Kronecker

C <- matrix(c(1, 0, 3, 4,
              0, 4, 1, -1,
              1, 1, -3, 2), nrow=3, byrow = T)
D <- c(1, 3, 7)
kronecker(C, D)

##       [,1] [,2] [,3] [,4]
##  [1,]    1    0    3    4
##  [2,]    3    0    9   12
##  [3,]    7    0   21   28
##  [4,]    0    4    1   -1
##  [5,]    0   12    3   -3
##  [6,]    0   28    7   -7
##  [7,]    1    1   -3    2
##  [8,]    3    3   -9    6
##  [9,]    7    7  -21   14

1.10 Bentuk Kuadrat

quad <- function(x){
  for(i in 1:(nrow(x))){
    cat(paste0(x[i,i],"(X",i,")^2"))
    for(j in (i+1):(nrow(x))) 
      if((i != j) & (j <= (nrow(x)) )) 
        cat(paste0(" + ",(x[i,j] + x[j,i]), "(X",i,j,")"))
    cat("\n")
  }
    cat("\n")
}

A <- matrix(c(1, 2, 3,
              4, 2, 1,
              1, 1, 1), nrow=3, byrow = T)
B <- matrix(c(1, 2, 3, 4,
              5, 6, 7, 8,
              9, 10, 11, 12,
              13, 14, 15, 16), nrow=4, byrow=TRUE)
C <- matrix(c(1, 2,
              4, 2), nrow=2, byrow = T)

quad(A)

## 1(X1)^2 + 6(X12) + 4(X13)
## 2(X2)^2 + 2(X23)
## 1(X3)^2

quad(B)

## 1(X1)^2 + 7(X12) + 12(X13) + 17(X14)
## 6(X2)^2 + 17(X23) + 22(X24)
## 11(X3)^2 + 27(X34)
## 16(X4)^2

quad(C)

## 1(X1)^2 + 6(X12)
## 2(X2)^2

1.11 Yang blm paham

2.1 Uji normalitas ganda - Mardia’s Skewness

data <- matrix(c(8.8, 8.5, 7.7, 4.9, 9.6, 10, 11.5, 11.6, 11.2, 10.7, 10, 6.8, 
  2589, 1186, 291, 1276, 6633, 12125, 36717, 43319, 10530, 3931, 1536, 61400),
  nrow=12, ncol=2)

install_load('MVN')

## Warning: package 'MVN' was built under R version 4.2.3

mardia <- mvn(data, mvnTest = c("mardia"), covariance = TRUE, 
              multivariatePlot = "qq")

mardia

## $multivariateNormality
##              Test         Statistic            p value Result
## 1 Mardia Skewness  11.1215382968049 0.0252315192005968     NO
## 2 Mardia Kurtosis 0.590006258740303  0.555186453489953    YES
## 3             MVN              <NA>               <NA>     NO
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling  Column1     0.3122    0.5053    YES   
## 2 Anderson-Darling  Column2     1.3825    0.0008    NO    
## 
## $Descriptives
##    n      Mean   Std.Dev Median   Min     Max   25th      75th       Skew
## 1 12     9.275     2.041    9.8   4.9    11.6    8.3    10.825 -0.6760453
## 2 12 15127.750 20400.481 5282.0 291.0 61400.0 1471.0 18273.000  1.1413585
##     Kurtosis
## 1 -0.6824234
## 2 -0.2880875

Dari hasil R diatas akan dilakukan pengujian hipotesis terhadap mardia skewness test sebagai berikut :
Ho : Peubah ganda mengikuti distribusi normal
H1 : Peubah ganda tidak mengikuti distribusi normal
Sehingga dikesimpulan p-value = 0.02523 < α = 0.05 yaitu menolak Ho. Artinya bahwa peubah ganda tidak mengikuti distribusi normal. Dari Q-Q plot yang dihasilkan dari output R dibawah ini juga menunjukan bahwa sebaran dari data tidak mengikuti distribusi normal.

2.2 Uji normalitas ganda - Henze-Zikler Test

henze <- mvn(data, mvnTest = c("hz"), covariance = TRUE, multivariatePlot = "none")
henze

## $multivariateNormality
##            Test        HZ     p value MVN
## 1 Henze-Zirkler 0.9995781 0.004125476  NO
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling  Column1     0.3122    0.5053    YES   
## 2 Anderson-Darling  Column2     1.3825    0.0008    NO    
## 
## $Descriptives
##    n      Mean   Std.Dev Median   Min     Max   25th      75th       Skew
## 1 12     9.275     2.041    9.8   4.9    11.6    8.3    10.825 -0.6760453
## 2 12 15127.750 20400.481 5282.0 291.0 61400.0 1471.0 18273.000  1.1413585
##     Kurtosis
## 1 -0.6824234
## 2 -0.2880875

Dari Henze-Zirkler’s Multivariate Normality Test menghasilkan nilai p-value = 0.004125 < α = 0.05. Hal ini berarti data mendukung untuk menolak Ho, dengan demikian dapat dikatakan bahwa peubah-peubah tersebut tidak mengikuti distribusi normal ganda.

2.3 Uji normalitas ganda - Royston

royston <- mvn(data, mvnTest = c("royston"), covariance = TRUE, multivariatePlot = "persp")

royston

## $multivariateNormality
##      Test        H     p value MVN
## 1 Royston 10.10069 0.006407917  NO
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling  Column1     0.3122    0.5053    YES   
## 2 Anderson-Darling  Column2     1.3825    0.0008    NO    
## 
## $Descriptives
##    n      Mean   Std.Dev Median   Min     Max   25th      75th       Skew
## 1 12     9.275     2.041    9.8   4.9    11.6    8.3    10.825 -0.6760453
## 2 12 15127.750 20400.481 5282.0 291.0 61400.0 1471.0 18273.000  1.1413585
##     Kurtosis
## 1 -0.6824234
## 2 -0.2880875

Dari Royston Test menghasilkan nilai p-value = 0.0064 < α = 0.05. Hasil uji ini juga menunjukkan data mendukung untuk menolak Ho, dengan demikian dapat dikatakan bahwa peubah-peubah tersebut tidak mengikuti distribusi normal ganda.

royston <- mvn(data, mvnTest = c("royston"), covariance = TRUE,
               multivariatePlot = "contour")

Berdasarkan beberapa uji normalitas ganda di atas, dapat disimpulkan bahwa data tidak menyebar multivariat normal. Sehingga dapat dilakukan penangan lanjutan terhadap data tersebut, dapat menggunakan transformasi normal atau metode lainnya.

3.1 Data

install_load('readxl')
data <- read_excel("Data (3).xlsx")
install_load('DT')
datatable(data, filter = 'top', 
          options = list(pageLength = 10))

3.2 Vektor Rataan & Mariks Kovarians

xbar <- apply(data[,2:3], 2, mean)
xbar

## matematika     fisika 
##      62.58      66.60

cov_m <- cov(data[,2:3])
cov_m

##            matematika   fisika
## matematika  164.93289 36.00889
## fisika       36.00889 66.96444

3.3 Uji T2 Hotelling Satu Populasi

\[H_{0}:\mu_{0}=\begin{bmatrix} 55\\ 60 \end{bmatrix}\]

\[H_{1}:\mu_{0}\neq \begin{bmatrix} 55\\ 60 \end{bmatrix}\]

install_load('MVTests')

mean0 <- c(55,60)
result <- OneSampleHT2(data[,2:3],mu0=mean0, alpha=0.1)
summary(result)

##        One Sample Hotelling T Square Test 
## 
## Hotelling T Sqaure Statistic = 7.621161 
##  F value = 3.387 , df1 = 2 , df2 = 8 , p-value: 0.086 
## 
##                    Descriptive Statistics
## 
##       matematika    fisika
## N       10.00000 10.000000
## Means   62.58000 66.600000
## Sd      12.84262  8.183181
## 
## 
##                  Detection important variable(s)
## 
##               Lower    Upper Mu0 Important Variables?
## matematika 51.83163 73.32837  55                FALSE
## fisika     59.75125 73.44875  60                FALSE

install_load('ICSNP')
test <- HotellingsT2(data[,2:3],mu=mean0,test="f")
print(test)

## 
##  Hotelling's one sample T2-test
## 
## data:  data[, 2:3]
## T.2 = 3.3872, df1 = 2, df2 = 8, p-value = 0.08597
## alternative hypothesis: true location is not equal to c(55,60)

Kesimpulan : Tolak H0 pada taraf nyata 10% karena p-value < alpha
Artinya : Dengan tingkat kepercayaan 90%, dapat disimpulkan bahwa minimal ada salah satu dari rata-rata ujian matematika dan fisika yang memiliki nilai rata-rata tidak sama dengan nilai 55 atau nilai 60.

3.4 Elips Kepercayaan

install_load('ellipse')

n = 10
p = 2
plot(ellipse(cov_m,centre=xbar,level = 0.90, 
             t=sqrt(((n-1)*p/(n*(n-p)))*qf(0.90,p,n-p))),
     type="l",main = "Ellips Kepercayaan 90%")
points(xbar[1],xbar[2])

#Siswa dengan skor matematika 65 dan fisika 70
points(65,70,col = "red")

Berdasarkan gambar terlihat bahwa nilai tersebut dalam elips, dengan kata lain masuk ke dalam elips kepercayaan.

3.5 Selang Kepercayaan Simultan

T.ci = function(mu, Sigma, n, avec=rep(1,length(mu)), level=0.95){
p = length(mu)
cval = qf(level, p, n-p) * p * (n-1) / (n-p)
zhat = crossprod(avec, mu)
zvar = crossprod(avec, Sigma %*% avec) / n
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}

n = 10

#Matematika
T.ci(mu=xbar, Sigma=cov_m, n=n, avec=c(1,0),level=0.9)

##    lower    upper 
## 51.83163 73.32837

#Fisika
T.ci(mu=xbar, Sigma=cov_m, n=n, avec=c(0,1),level=0.9)

##    lower    upper 
## 59.75125 73.44875

3.6 Selang Kepercayaan Bonferroni

bon = function(mu,S,n,alpha,k){
 p = length(mu)
 lower = mu[k] - sqrt(S[k,k]/n) * abs(qt(alpha/(2*p), df=n-1))
 upper = mu[k] + sqrt(S[k,k]/n) * abs(qt(alpha/(2*p), df=n-1))
 c(lower = lower,upper = upper)
}

#Matematika
bon(xbar, cov_m,10,0.1,1)

## lower.matematika upper.matematika 
##         53.39294         71.76706

#Fisika
bon(xbar, cov_m,10,0.1,2)

## lower.fisika upper.fisika 
##     60.74611     72.45389

5.1 Soal 1

install_load('car')
jagung = read.table(header=T,text =
               "
Varietas Y1 Y2
v1 6 7
v1 5 9
v2 4 6
v2 6 6
v2 4 7
v3 5 4
v3 6 4
")
datatable(jagung, filter = 'top', 
          options = list(pageLength = 10))

mod1 = Manova(lm(cbind(Y1,Y2)~Varietas,data=jagung),type="III")
summary(mod1,multivariate=TRUE)

## 
## Type III MANOVA Tests:
## 
## Sum of squares and products for error:
##           Y1        Y2
## Y1  3.666667 -1.666667
## Y2 -1.666667  2.666667
## 
## ------------------------------------------
##  
## Term: (Intercept) 
## 
## Sum of squares and products for the hypothesis:
##      Y1  Y2
## Y1 60.5  88
## Y2 88.0 128
## 
## Multivariate Tests: (Intercept)
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1   0.99248      198      2      3 0.00065196 ***
## Wilks             1   0.00752      198      2      3 0.00065196 ***
## Hotelling-Lawley  1 132.00000      198      2      3 0.00065196 ***
## Roy               1 132.00000      198      2      3 0.00065196 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: Varietas 
## 
## Sum of squares and products for the hypothesis:
##            Y1         Y2
## Y1  1.1904762 -0.4761905
## Y2 -0.4761905 16.1904762
## 
## Multivariate Tests: Varietas
##                  Df test stat  approx F num Df den Df   Pr(>F)  
## Pillai            2  1.138478  2.642948      4      8 0.112828  
## Wilks             2  0.080460  3.788127      4      6 0.071870 .
## Hotelling-Lawley  2  8.707483  4.353741      4      4 0.091633 .
## Roy               2  8.382882 16.765765      2      4 0.011359 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2 Soal 4

diet = read.table(header=T,text =
               "
Kel P   K   KP
K1  20  5   18
K1  25  9   8
K1  23  15  20
K1  16  9   22
K1  20  6   22
K2  28  7   14
K2  25  14  5
K2  26  9   20
K2  19  15  22
K2  29  14  12
K3  15  6   3
K3  22  8   12
K3  27  9   14
K3  21  10  7
K3  17  9   1
")

datatable(diet, filter = 'top', 
          options = list(pageLength = 10))

mod2 = Manova(lm(cbind(P,K,KP)~Kel,data=diet),type="III")
summary(mod2,multivariate=TRUE)

## 
## Type III MANOVA Tests:
## 
## Sum of squares and products for error:
##        P     K    KP
## P  195.2   6.4 -15.0
## K    6.4 120.8  -7.2
## KP -15.0  -7.2 444.4
## 
## ------------------------------------------
##  
## Term: (Intercept) 
## 
## Sum of squares and products for the hypothesis:
##         P     K   KP
## P  2163.2 915.2 1872
## K   915.2 387.2  792
## KP 1872.0 792.0 1620
## 
## Multivariate Tests: (Intercept)
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1  0.948328 61.17669      3     10 9.7542e-07 ***
## Wilks             1  0.051672 61.17669      3     10 9.7542e-07 ***
## Hotelling-Lawley  1 18.353006 61.17669      3     10 9.7542e-07 ***
## Roy               1 18.353006 61.17669      3     10 9.7542e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------
##  
## Term: Kel 
## 
## Sum of squares and products for the hypothesis:
##       P        K        KP
## P  77.2 51.60000  41.00000
## K  51.6 34.53333  30.86667
## KP 41.0 30.86667 292.93333
## 
## Multivariate Tests: Kel
##                  Df test stat approx F num Df den Df   Pr(>F)  
## Pillai            2 0.7797319 2.342942      6     22 0.066826 .
## Wilks             2 0.3653752 2.181205      6     20 0.088300 .
## Hotelling-Lawley  2 1.3397669 2.009650      6     18 0.117387  
## Roy               2 0.8970362 3.289133      3     11 0.061957 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.1 Data

dataset <- read.csv("Data.csv",sep = ";")
datatable(dataset, filter = 'top', 
          options = list(pageLength = 10))

6.2 Plot

install_load('profileR')
mod <- pbg(data=dataset[,1:4], group=dataset[,5], original.name = TRUE, profile.plot = TRUE)

Berdasarkan grafik terlihat bahwa profil pikun sejajar dengan profil tidak pikun. Selain ini kedua profil tidak berhimpit. Dengan kata lain berdasarkan grafik peningkatan faktor pikun dan tidak pikun terhadap beberapa subtes sama. Sedangkan rata-rata hasil skor pikun berbeda dengan yang tidak pikun, di mana yang tidak pikun lebih tinggi.

6.3 Summary

# rata-rata tiap perlakuan untuk dua kelompok
print(mod)

## 
## Data Summary:
##                Pikun Tidak pikun
## Informasi   8.750000   12.567568
## Similaritas 5.333333    9.486486
## Aritmetik   8.500000   11.513514
## Gambar      4.750000    7.972973

# semua perhitungan analisis profil
profil <- summary(mod)

## Call:
## pbg(data = dataset[, 1:4], group = dataset[, 5], original.names = TRUE, 
##     profile.plot = TRUE)
## 
## Hypothesis Tests:
## $`Ho: Profiles are parallel`
##   Multivariate.Test  Statistic  Approx.F num.df den.df   p.value
## 1             Wilks 0.97461349 0.3907166      3     45 0.7602428
## 2            Pillai 0.02538651 0.3907166      3     45 0.7602428
## 3  Hotelling-Lawley 0.02604778 0.3907166      3     45 0.7602428
## 4               Roy 0.02604778 0.3907166      3     45 0.7602428
## 
## $`Ho: Profiles have equal levels`
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## group        1  114.3  114.31   17.44 0.000128 ***
## Residuals   47  308.1    6.56                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`Ho: Profiles are flat`
##          F df1 df2      p-value
## 1 51.80163   3  45 1.223077e-14

6.4 Uji Kesejajaran (Parallel Test)

\[H_{0}:C\mu_{1}=C\mu_{2}\]

\[H_{1}:C\mu_{1}\neq C\mu_{2}\]

profil$`Ho: Profiles are parallel`

##   Multivariate.Test  Statistic  Approx.F num.df den.df   p.value
## 1             Wilks 0.97461349 0.3907166      3     45 0.7602428
## 2            Pillai 0.02538651 0.3907166      3     45 0.7602428
## 3  Hotelling-Lawley 0.02604778 0.3907166      3     45 0.7602428
## 4               Roy 0.02604778 0.3907166      3     45 0.7602428

Berdasarkan statistik di atas, diperoleh F hitung sebesar 0.39 dengan p-value sebesar 0.76 sehingga pada taraf signifikansi 5%, H0 gagal ditolak. Jadi, dapat disimpulkan bahwa profil antara kelompok Usia (Pikun dan Tidak Pikun) adalah sejajar.

6.5 Uji Keberhimpitan

\[H_{0}:1'\mu_{1}=1'\mu_{2}\]

\[H_{1}:1'\mu_{1}\neq 1'\mu_{2}\]

profil$`Ho: Profiles have equal levels`

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## group        1  114.3  114.31   17.44 0.000128 ***
## Residuals   47  308.1    6.56                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan statistik di atas, diperoleh F hitung sebesar 17.44 dengan p-value 0.000128 sehingga pada taraf signifikansi 5%, H0 ditolak. Jadi, dapat disimpulkan bahwa profil antara kelompok Usia (Pikun dan Tidak Pikun) tidak berhimpit.

6.6 Uji Kesamaan

\[H_{0}:C(\mu_{1}-\mu_{2})=0\] \[H_{1}:C(\mu_{1}-\mu_{2})\neq 0\]

profil$`Ho: Profiles are flat`

##          F df1 df2      p-value
## 1 51.80163   3  45 1.223077e-14

Berdasarkan statistik di atas, diperoleh F hitung sebesar 51.80 dengan p-value 1.223077e-14 sehingga pada taraf signifikansi 5%, H0 ditolak. Jadi, dapat disimpulkan bahwa rata-rata untuk setiap sub tes pada lansia kelompok yang ada faktor kepikunan dan kelompok yang tidak ada faktor kepikunan menunjukkan konstanta yang berbeda.

7.1 Data Pelari Wanita

Berikut adalah data catatan waktu hasil tujuh nomor cabang lari atletik wanita yang berasal dari 55 negara pada salah satu event olimpiade yaitu lari 100 meter, 200 meter, 400 meter, 800 meter, 1500 meter, 3000 meter, dan maraton. Tiga nomor cabang lari pertama dicatat dalam satuan detik, sedangkan empat nomor yang lain dalam menit.

install_load('dplyr')
data_women_records <- read_excel("women_track_records.xlsx") %>% as.data.frame()
datatable(data_women_records, filter = 'top', 
          options = list(pageLength = 10))

rownames(data_women_records) <- data_women_records$country
data_women_records <- data_women_records[,-8]

7.2 Eksplorasi

cor_women <- cor(data_women_records)
ggcorrplot(cor_women,type="lower",lab = TRUE)

7.3 Menerapkan PCA (AKU)

Dalam R, Penerapan PCA ini dapat dilakukan dengan menggunakan fungsi prcomp. Fungsi ini memiliki argumen scale dan center. Jika kedua argumen ini TRUE maka matrix yang digunakan untuk menghitung PCA adalah matrix korelasi. Namun, jika kedua argumen ini FALSE atau scale=FALSE, maka matrix yang digunakan adalah matrix covariance.

pca_women_records <- prcomp(data_women_records,scale=TRUE,center=TRUE)
summary(pca_women_records)

## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4095 0.80848 0.54762 0.35423 0.23198 0.19761 0.14981
## Proportion of Variance 0.8294 0.09338 0.04284 0.01793 0.00769 0.00558 0.00321
## Cumulative Proportion  0.8294 0.92276 0.96560 0.98353 0.99122 0.99679 1.00000

Hasil yang dikeluarkan dari sintaks di atas terdiri dari tiga macam, yaitu Standard deviation, Proportion of Variance dan Cumulative Proportion dari masing-masing Komponen Utama (Principal Component).

• Standard deviaton merupakan akar dari akar ciri (eigenvalue). Dalam hal ini akar ciri berperan sebagai variance dari masing-masing komponen utama.

• Proportion of Variance didapatkan dari akar ciri pada masing-masing komponen dibagi dengan total akar ciri. Proportion of Variance menjelaskan seberapa besar keragaman peubah asal yang dapat dijelaskan oleh masing-masing komponen utama. Semakin besar nilainya berarti semakin baik pula komponen utama tersebut untuk merepresentasikan peubah asal.

• Cumulative Proportion menjelaskan seberapa besar keragaman yang dapat dijelaskan oleh komponen utama secara kumulatif. Misalnya saja dengan menggunakan dua komponen utama saja (PC1 dan PC2), sudah bisa menjelaskan lebih dari 92% keragaman dari data.
  Berdasarkan hal ini, kita akan memilih menggunakan dua komponen utama saja.

fviz_screeplot(pca_women_records,geom="line")

Hal lain yang bisa dilakukan untuk menentukan berapa banyak komponen utama yang digunakan adalah dengan screeplot. Fungsi untuk menampilkan screeplot pada R adalah fviz_screeplot yang didapat dari package factoextra. Banyaknya komponen utama bisa ditentukan dengan screeplot dengan melihat di komponen utama yang mana garisnya berbentuk seperti siku (elbow).
Pada gambar diatas, garis membentuk siku saat berada di komponen utama kedua (dimension kedua) sehingga banyaknya komponen utama yang digunakan sebanyak dua (Komponen Utama 1 dan Komponen Utama 2).

7.4 Interpretasi PCA (AKU)

Interpretasi metode PCA dapat dilakukan dengan menggunakan vektor ciri pada masing-masing komponen utama. Semakin besar vektor ciri pada komponen utama tertentu, maka semakin besar pula kontribusi dari peubah asal untuk membangun komponen utama tersebut. Catatan lain yang perlu diperhatikan adalah nilai negatif pada vektor ciri menandakan peubah asal memberikan kontribusi yang berkembalikan pada pembentukan komponen utama. Dalam konteks vektor ciri negatif, semakin besar nilai peubah asal semakin kecil nilai pada komponen utama.

pca_women_records$rotation

##                PC1        PC2         PC3         PC4         PC5          PC6
## 100m     0.3683561  0.4900597 -0.28601157  0.31938631  0.23116950  0.619825234
## 200m     0.3653642  0.5365800 -0.22981913 -0.08330196  0.04145457 -0.710764580
## 400m     0.3816103  0.2465377  0.51536655 -0.34737748 -0.57217791  0.190945970
## 800m     0.3845592 -0.1554023  0.58452608 -0.04207636  0.62032379 -0.019089032
## 1500m    0.3891040 -0.3604093  0.01291198  0.42953873  0.03026144 -0.231248381
## 3000m    0.3888661 -0.3475394 -0.15272772  0.36311995 -0.46335476  0.009277159
## Marathon 0.3670038 -0.3692076 -0.48437037 -0.67249685  0.13053590  0.142280558
##                  PC7
## 100m      0.05217655
## 200m     -0.10922503
## 400m      0.20849691
## 800m     -0.31520972
## 1500m     0.69256151
## 3000m    -0.59835943
## Marathon  0.06959828

Karena kita hanya menggunakan dua komponen saja, maka vector ciri yang akan dinterpretasikan hanya pada PC1 dan PC2. PC1 memiliki vektor ciri yang relatif sama yaitu berkisar di 0.3 untuk semua cabang lomba. Vektor ciri yang relatif sama ini menandakan bahwa kontribusi peubah asal untuk membangun komponen utama ini relatif sama. Artinya nilai-nilai yang ada di PC1 (score value) dapat menggambarkan waktu lari untuk semua cabang lomba. Oleh karena itu kita dapat dapat menggunakan PC1 untuk menentukan negara mana yang memiliki pelari tercepat untuk semua kategori lomba.
Vektor ciri di PC2 memiliki nilai positif untuk cabang lari jarak dekat (100m -400m) dan nilai negatif untuk cabang lari jarak jauh (800m-marathon). Hal ini berarti semakin besar score value pada PC2 maka waktu lari cabang jarak dekat semakin lambat namun waktu lari untuk cabang jarak jauh semakin cepat. Oleh karena itu, PC2 dapat digunakan untuk menentukan negara mana yang pada cabang lari jarak dekat waktunya mirip seperti cabang lari jarak jauh.

Note: Interpretasi komponen utama memiliki subjektifitas yang tinggi, oleh karena itu setiap orang menginterpretasikanya berbeda.

Hal terakhir yang bisa diinterpretasikan adalah score value pada PC1 dan PC2. Score value merupakan observasi/koordinat baru pada peubah komponen utama. Dalam konteks data pelari diatas, observasinya adalah negara, sehingga kita dapat memberi insight cabang perlombaan lari dari setiap negara. Untuk melihat score value pada komponen utama dapat dilihat dengan menggunakan sintaks berikut.

datatable(pca_women_records$x, filter = 'top', 
          options = list(pageLength = 10))

Agar lebih mudah dalam menginterpretasikan score value maka digunakaan grafik di bawah ini.

fviz_pca_ind(pca_women_records,col.ind = "darkred")

Berdasarkan grafik score value dapat diketahui bahwa negara yang memiliki catatan waktu pelari terlambat untuk semua cabang lomba adalah negara wsamoa. Hal ini dikarenakan wsamoa score value wsakoa untuk PC1 (Dim1) paling besar diantara yang lain. Walaupun negara wsamoa memiliki cabang lari terlama disemua cabang lomba, namun perbedaan waktu terkecil antara pelari jarak jauh dan jarak dekat adalah negara wsamoa. Hal ini berarti pelari untuk lomba jarak dekat sangat lambat karena memiliki waktu yang hampir mirip seperti pelari jarak jauh. Sedangkan negara yang memiliki pelari tercepat untuk semua cabang lomba adalah gdr.