Pertemuan 5 - Bu Sarni
Qorinul Huda
2022-09-19
#Lihat di File dasar Teori #Paired Observation Test (Univariate)
Orang.ke=c(seq(1:10))
BB.Sebelum=c(57,69,56,67,55,56,62,67,67,56)
BB.Sesudah=c(55,70,56,65,54,55,64,65,67,54)
data=data.frame(Orang.ke,BB.Sebelum,BB.Sesudah)
t.test(x=data$BB.Sebelum, y=data$BB.Sesudah,
alternative = "greater",
mu = 0.5, paired = TRUE, var.equal = TRUE,
conf.level = 0.95)##
## Paired t-test
##
## data: data$BB.Sebelum and data$BB.Sesudah
## t = 0.44598, df = 9, p-value = 0.3331
## alternative hypothesis: true mean difference is greater than 0.5
## 95 percent confidence interval:
## -0.1220671 Inf
## sample estimates:
## mean difference
## 0.7#Paired Observation Test (multivariate) #Dataset
library(readxl)
soal1 <- read_excel("F:/APG_Praktikum/APG_Bu Sarni/Pertemuan 5/dataset 5.xlsx",
sheet = "No1")
soal1## # A tibble: 15 x 4
## kata1 katakerja1 kata2 katakerja2
## <dbl> <dbl> <dbl> <dbl>
## 1 148 20 137 15
## 2 159 24 164 25
## 3 144 19 224 27
## 4 103 18 208 33
## 5 121 17 178 24
## 6 89 11 128 20
## 7 119 17 154 18
## 8 123 13 158 16
## 9 76 16 102 21
## 10 217 29 214 25
## 11 148 22 209 24
## 12 151 21 151 16
## 13 83 7 123 13
## 14 135 20 161 22
## 15 178 15 175 23#memakai fungsi Hotelling pertemuan kita sebelumnya. Bedanya mdata yang dimasukkan adalah difference
Hotelling<-function(mdata,Ho,alfa=0.05){
n<-nrow(mdata)
p<-ncol(mdata)
#Membentuk vektor rata-rata sampel
xbar_mdata=apply(mdata, 2, mean)
varcovar=cov(mdata)
#T hotelling
Tsq<-n*(t(xbar_mdata-Ho)%*%solve(varcovar)%*%(xbar_mdata-Ho))
#wilayah kritis
Ftab<-qf(1-alfa,p,n-p)
Fmod<-(p*(n-1)/(n-p))*Ftab
uji<-ifelse(Tsq>Fmod, "Tolak Ho","Gagal Tolak Ho")
library("DescTools")
auto<-HotellingsT2Test(mdata, mu=Ho, test = 'chi')
hasil<-list(xbar_mdata, varcovar, Tsq, Ftab, Fmod, uji,auto)
names(hasil)<-c("vektor rata-rata","matriks var-covar","nilai T2 Hotelling",
"F tabel","F modifikasi","Keputusan","Hasil uji otomatis")
return(hasil)
}#Penggunaan
#Buat matriks selisih
d = matrix(c(soal1$kata1-soal1$kata2,soal1$katakerja1-soal1$katakerja2),15,2)
d## [,1] [,2]
## [1,] 11 5
## [2,] -5 -1
## [3,] -80 -8
## [4,] -105 -15
## [5,] -57 -7
## [6,] -39 -9
## [7,] -35 -1
## [8,] -35 -3
## [9,] -26 -5
## [10,] 3 4
## [11,] -61 -2
## [12,] 0 5
## [13,] -40 -6
## [14,] -26 -2
## [15,] 3 -8#Hipotesis
Ho = matrix(c(0,0),2,1)
Hotelling(d,Ho)## $`vektor rata-rata`
## [1] -32.800000 -3.533333
##
## $`matriks var-covar`
## [,1] [,2]
## [1,] 1096.029 139.90000
## [2,] 139.900 31.55238
##
## $`nilai T2 Hotelling`
## [,1]
## [1,] 15.19123
##
## $`F tabel`
## [1] 3.805565
##
## $`F modifikasi`
## [1] 8.196602
##
## $Keputusan
## [,1]
## [1,] "Tolak Ho"
##
## $`Hasil uji otomatis`
##
## Hotelling's one sample T2-test
##
## data: mdata
## T.2 = 15.191, df = 2, p-value = 0.0005026
## alternative hypothesis: true location is not equal to c(0,0)#Otomatis
#Pisahkan dua populasi
pop1 = matrix(c(soal1$kata1,soal1$katakerja1),15,2)
pop2 = matrix(c(soal1$kata2,soal1$katakerja2),15,2)
library(MVTests)##
## Attaching package: 'MVTests'## The following object is masked from 'package:datasets':
##
## iris#Cara 1
Mpaired(pop1, pop2)## $HT2
## [,1]
## [1,] 15.19123
##
## $F
## [,1]
## [1,] 7.053073
##
## $df
## [1] 2 13
##
## $p.value
## [,1]
## [1,] 0.008427354
##
## $Descriptive1
## [,1] [,2]
## Means 132.93333 17.933333
## Sd 37.49375 5.351457
##
## $Descriptive2
## [,1] [,2]
## Means 165.73333 21.466667
## Sd 36.05247 5.289702
##
## $Descriptive.Difference
## [,1] [,2]
## Means 32.80000 3.533333
## Sd 33.10632 5.617151
##
## $Test
## [1] "Mpaired"
##
## attr(,"class")
## [1] "MVTests" "list"#Cara 2
OneSampleHT2(d, Ho, alpha = 0.05)## $HT2
## [,1]
## [1,] 15.19123
##
## $F
## [,1]
## [1,] 7.053073
##
## $df
## [1] 2 13
##
## $p.value
## [,1]
## [1,] 0.008427354
##
## $CI
## Lower Upper Mu0 Important Variables?
## 1 -57.27272 -8.3272804 0 *TRUE*
## 2 -7.68562 0.6189537 0 FALSE
##
## $alpha
## [1] 0.05
##
## $Descriptive
## [,1] [,2]
## N 15.00000 15.000000
## Means -32.80000 -3.533333
## Sd 33.10632 5.617151
##
## $Test
## [1] "OneSampleHT2"
##
## attr(,"class")
## [1] "MVTests" "list"#Selang Simultan Paired Data
selang_paired<-function(d,alfa=0.05){
#rata-rata d
rd = matrix(c(colMeans(d)),2,1)
n = nrow(d)
p = ncol(d)
sd=var(d)
#Definisi variabel loop
ci_lower=matrix()
ci_upper=matrix()
ci=matrix()
for (i in 1 : p){
ci_lower=rd[i,]-sqrt((n-1)*p*qf(1-alfa,p,n-p)/(n-p))*sqrt(sd[i,i]/n)
ci_upper=rd[i,]+sqrt((n-1)*p*qf(1-alfa,p,n-p)/(n-p))*sqrt(sd[i,i]/n)
ci<-c(ci_lower, ci_upper)
print(ci)
}
}
selang_paired(d)## [1] -57.27272 -8.32728
## [1] -7.6856203 0.6189537#Coba di PPT tuh
difference<-matrix(c(-19,-22,-18,-27,-4,-10,-14,17,9,4,-19,
+ 12,10,42,15,-1,11,-4,60,-2,10,-7),11,2)
#Hasilnya sama tuh selangnya#Repeated Measure () #Dataset
library(readxl)
soal2 <- read_excel("F:/APG_Praktikum/APG_Bu Sarni/Pertemuan 5/dataset 5.xlsx",
sheet = "No2")
soal2## # A tibble: 19 x 4
## `1` `2` `3` `4`
## <dbl> <dbl> <dbl> <dbl>
## 1 426 609 556 600
## 2 253 236 392 395
## 3 359 433 349 357
## 4 432 431 522 600
## 5 405 426 513 513
## 6 324 438 507 539
## 7 310 312 410 456
## 8 326 326 350 504
## 9 375 447 547 548
## 10 286 286 403 422
## 11 349 382 473 497
## 12 429 410 488 547
## 13 348 377 447 514
## 14 412 473 472 446
## 15 347 326 455 468
## 16 434 458 637 524
## 17 364 367 432 469
## 18 420 395 508 531
## 19 397 556 645 625#Membuat matriks kontras
Cn = matrix(c(-1,1,1,-1,-1,-1,1,1,-1,1,-1,1),3,4)
Cn## [,1] [,2] [,3] [,4]
## [1,] -1 -1 1 1
## [2,] 1 -1 1 -1
## [3,] 1 -1 -1 1mperlakuan<-function(mdt,C,alfa=0.05){
n = nrow(mdt)
q = ncol(mdt)
xbar=colMeans(mdt)
varcovar=cov(mdt)
cxbar=C%*%xbar
cSc=C%*%varcovar%*%t(C)
#Statistik Uji
tsq=n*t(cxbar)%*%solve(cSc)%*%cxbar
Fmod=((n-1)*(q-1))/(n-q+1)
Ftab=qf(1-alfa,q-1,n-q+1)
Fstat=Fmod*Ftab
uji<-ifelse(tsq>Fstat, "Tolak Ho","Gagal Tolak Ho")
hasil<-list(xbar,varcovar,cxbar,cSc,tsq,Fmod,Ftab,Fstat,uji)
names(hasil)<-c("Vektor rata-rata","Var-Covar","C*xbar",
"CSC'","T-square","F mod","F tabel",
"Statistik uji","Keputusan")
return(hasil)
}
mperlakuan(soal2,Cn)## $`Vektor rata-rata`
## 1 2 3 4
## 368.2105 404.6316 479.2632 502.8947
##
## $`Var-Covar`
## 1 2 3 4
## 1 2819.287 3568.415 2943.497 2295.357
## 2 3568.415 7963.135 5303.991 4065.459
## 3 2943.497 5303.991 6851.316 4499.640
## 4 2295.357 4065.459 4499.640 4878.988
##
## $`C*xbar`
## [,1]
## [1,] 209.31579
## [2,] -60.05263
## [3,] -12.78947
##
## $`CSC'`
## [,1] [,2] [,3]
## [1,] 9432.2281 1098.9064 927.5965
## [2,] 1098.9064 5195.8304 914.5673
## [3,] 927.5965 914.5673 7557.3977
##
## $`T-square`
## [,1]
## [1,] 116.0163
##
## $`F mod`
## [1] 3.375
##
## $`F tabel`
## [1] 3.238872
##
## $`Statistik uji`
## [1] 10.93119
##
## $Keputusan
## [,1]
## [1,] "Tolak Ho"#Selang data treatment (cxbar[i]) ± sqrt(Fstat)*sqrt(cSc[i,i]/n)
selang_treatment<-function(mdt,C,alfa=0.05){
#identifikasi dataset
n = nrow(mdt)
q = ncol(mdt)
xbar=colMeans(mdt)
varcovar=cov(mdt)
cxbar=C%*%xbar
cSc=C%*%varcovar%*%t(C)
#Statistik F
Fmod=((n-1)*(q-1))/(n-q+1)
Ftab=qf(1-alfa,q-1,n-q+1)
Fstat=Fmod*Ftab
#Definisi variabel loop
ci_lower=matrix()
ci_upper=matrix()
ci=matrix()
for (i in 1 : nrow(C)){
ci_lower=cxbar[i,]-sqrt(Fstat)*sqrt(cSc[i,i]/n)
ci_upper=cxbar[i,]+sqrt(Fstat)*sqrt(cSc[i,i]/n)
ci<-c(ci_lower, ci_upper)
print(ci)
}
}
selang_treatment(soal2,Cn)## [1] 135.6503 282.9813
## [1] -114.72708 -5.37818
## [1] -78.72858 53.14964#One-way Anova
#Kita mencoba dengan dataset iris
iris=data.frame(iris)
head(iris,10)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 Setosa
## 2 4.9 3.0 1.4 0.2 Setosa
## 3 4.7 3.2 1.3 0.2 Setosa
## 4 4.6 3.1 1.5 0.2 Setosa
## 5 5.0 3.6 1.4 0.2 Setosa
## 6 5.4 3.9 1.7 0.4 Setosa
## 7 4.6 3.4 1.4 0.3 Setosa
## 8 5.0 3.4 1.5 0.2 Setosa
## 9 4.4 2.9 1.4 0.2 Setosa
## 10 4.9 3.1 1.5 0.1 Setosa#Kita akan mengambil kelompok virginica dan setosa dari dataset iris
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionvirginica <- filter(iris, Species == "virginica")
setosa <- filter(iris, Species == "setosa")
#Misalkan kita ambil variabel Sepal length
anov = data.frame(log(iris$Sepal.Length),iris$Species)
names(anov)=c("sepal.length","spesies")
head(anov,10)## sepal.length spesies
## 1 1.629241 Setosa
## 2 1.589235 Setosa
## 3 1.547563 Setosa
## 4 1.526056 Setosa
## 5 1.609438 Setosa
## 6 1.686399 Setosa
## 7 1.526056 Setosa
## 8 1.609438 Setosa
## 9 1.481605 Setosa
## 10 1.589235 Setosasummary(anov)## sepal.length spesies
## Min. :1.459 Setosa :50
## 1st Qu.:1.629 Versicolor:50
## Median :1.758 Virginica :50
## Mean :1.755
## 3rd Qu.:1.856
## Max. :2.067#Uji kehomogenan
#Levine Test
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:DescTools':
##
## RecodeleveneTest(sepal.length ~ spesies, data = anov) ## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 2.2006 0.1144
## 147#Uji normalitas
shapiro.test(anov$sepal.length)##
## Shapiro-Wilk normality test
##
## data: anov$sepal.length
## W = 0.98253, p-value = 0.05388qqPlot(anov$sepal.length)
## [1] 132 14# lm: kode instruksi untuk penyusunan model linier dalam analisa data
model = lm(sepal.length ~ spesies, data = anov)
# anova: kode instruksi untuk analisa ragam satu jalur (one way)
anova(model)## Analysis of Variance Table
##
## Response: sepal.length
## Df Sum Sq Mean Sq F value Pr(>F)
## spesies 2 1.8918 0.94588 128.93 < 2.2e-16 ***
## Residuals 147 1.0784 0.00734
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Manova one way
presentasi <- read.csv("F:/APG_Praktikum/APG_Bu Sarni/Pertemuan 5/Manova_data.csv")
presentasi## Group X1 X2 X3
## 1 1 19.6 5.15 9.5
## 2 1 15.4 5.75 9.1
## 3 1 22.3 4.35 3.3
## 4 1 24.3 7.55 5.0
## 5 1 22.5 8.50 6.0
## 6 1 20.5 10.25 5.0
## 7 1 14.1 5.95 18.8
## 8 1 13.0 6.30 16.5
## 9 1 14.1 5.45 8.9
## 10 1 16.7 3.75 6.0
## 11 1 16.8 5.10 7.4
## 12 2 17.1 9.00 7.5
## 13 2 15.7 5.30 8.5
## 14 2 14.9 9.85 6.0
## 15 2 19.7 3.60 2.9
## 16 2 17.2 4.05 0.2
## 17 2 16.0 4.40 2.6
## 18 2 12.8 7.15 7.0
## 19 2 13.6 7.25 3.2
## 20 2 14.2 5.30 6.2
## 21 2 13.1 3.10 5.5
## 22 2 16.5 2.40 6.6
## 23 3 16.0 4.55 2.9
## 24 3 12.5 2.65 0.7
## 25 3 18.5 6.50 5.3
## 26 3 19.2 4.85 8.3
## 27 3 12.0 8.75 9.0
## 28 3 13.0 5.20 10.3
## 29 3 11.9 4.75 8.5
## 30 3 12.0 5.85 9.5
## 31 3 19.8 2.85 2.3
## 32 3 16.5 6.55 3.3
## 33 3 17.4 6.60 1.9str(presentasi)## 'data.frame': 33 obs. of 4 variables:
## $ Group: int 1 1 1 1 1 1 1 1 1 1 ...
## $ X1 : num 19.6 15.4 22.3 24.3 22.5 20.5 14.1 13 14.1 16.7 ...
## $ X2 : num 5.15 5.75 4.35 7.55 8.5 ...
## $ X3 : num 9.5 9.1 3.3 5 6 5 18.8 16.5 8.9 6 ...library(dplyr)
grup1 <- filter(presentasi, Group == "1")
grup2 <- filter(presentasi, Group == "2")
grup3 <- filter(presentasi, Group == "3")
#Sample mean vector
mean_1 = matrix(c(colMeans(grup1[,c(2:4)])),3,1)
mean_2 =matrix(c(colMeans(grup2[,c(2:4)])),3,1)
mean_3 = matrix(c(colMeans(grup3[,c(2:4)])),3,1)
#jumlah n dan g
n1 = nrow(grup1)
n2 = nrow(grup2)
n3 = nrow(grup3)
#Sample covariance matrix
var1 = var(grup1[,c(2:4)])
var2 = var(grup2[,c(2:4)])
var3 = var(grup3[,c(2:4)])
#Pooled W
W = (n1-1)*var1+(n2-1)*var2+(n3-1)*var3
W## X1 X2 X3
## X1 293.965455 6.550909 -207.77727
## X2 6.550909 126.287273 34.18591
## X3 -207.777273 34.185909 426.37091xbar=matrix(c((n1*mean_1+n2*mean_2+n3*mean_3)/(n1+n2+n3)),3,1)
xbar## [,1]
## [1,] 16.330303
## [2,] 5.715152
## [3,] 6.475758#B
B = n1*(mean_1-xbar)%*%t(mean_1-xbar)+n2*(mean_2-xbar)%*%t(mean_2-xbar)+n3*(mean_3-xbar)%*%t(mean_3-xbar)
B## [,1] [,2] [,3]
## [1,] 52.92424 14.243939 64.55152
## [2,] 14.24394 3.975152 16.71121
## [3,] 64.55152 16.711212 81.82970#Wilks Lambda
W = det(W)/(det(B+W))
W## [1] 0.5257884