Diskusi 1

Bangkitkan dua gugus data berikut:

Data dari sebaran simetris
Data campuran:

50% dari sebaran normal + 50% dari sebaran chi square
50% sebaran chisquare dengan parameter a + 50% sebaran chisquare dengan parameter b
25% sebaran chisquare dengan parameter a + 25% sebaran chisquare dengan parametr b + 25% sebaran normal dengan parameter a + 25% sebaran normal dengan parameter b

Ambil sampel dengan n=4,12,20,60,100
Pada n berapa sebaran rataan dari masing-masing data mulai simetris/mendekati sebaran normal? Apakah data yang simetris lebih cepat mendekati sebaran normal? (note: lebih cepat berarti dengna n yang tidak terlalu banyak sudah mulai mendekati sebaran normal).

Jawaban:

Untuk sebaran simetris, sebaran yang akan digunakan adalah sebaran Uniform(0,1).

N <- 10^6

Untuk data campuran, a dan b untuk sebaran chisquare adalah 1 dan 2. Sendangkan untuk sebaran normal, parameter a dan b yang digunakan adalah μ=0,σ2=1 dan μ=1,σ2=2

N <- 10^6
uniform_pop <- runif(N)
dataGab1 <- c(rnorm(0.5*N, 0,1), rchisq(0.5*N, df = 1))
dataGab2 <- c(rchisq(0.5*N, 1), rchisq(0.5*N, 2))
dataGab3 <- c(rchisq(0.25*N, 1), rchisq(0.25*N, 2), rnorm(0.5*N), rnorm(0.5*N, 1, 2))

Sampel sebanyak n diambil dari masing-masing populasi.

set.seed(1234)

n <- c(4, 12, 20, 60, 100)
reps <- 10^2

sampSim <- list(NA)
sampPop1 <- list(NA)
sampPop2 <- list(NA)
sampPop3 <- list(NA)

for (i in 1:length(n)) {
    sampSim[[i]] <- matrix(sample(uniform_pop, size = n[i]*reps, replace = TRUE), ncol = n[i], byrow = TRUE)
    sampPop1[[i]] <- matrix(sample(dataGab1, size = n[i]*reps, replace = TRUE), ncol = n[i], byrow = TRUE)
    sampPop2[[i]] <- matrix(sample(dataGab2, size = n[i]*reps, replace = FALSE), ncol = n[i], byrow = TRUE)
    sampPop3[[i]] <- matrix(sample(dataGab3, size = n[i]*reps, replace = FALSE), ncol = n[i], byrow = TRUE)
}

apply(sampSim[[1]], 1, mean)

##   [1] 0.1455706 0.4408534 0.5055819 0.5467021 0.5247918 0.6130597 0.4619818
##   [8] 0.4511359 0.3727174 0.6013238 0.6460694 0.7657099 0.5537328 0.3252087
##  [15] 0.7168864 0.6897833 0.3889895 0.5688048 0.5659902 0.4658998 0.5257760
##  [22] 0.3415701 0.6632682 0.8130515 0.1846734 0.3761575 0.5333111 0.2877637
##  [29] 0.6535122 0.4961465 0.4191465 0.4922921 0.6266047 0.7496451 0.6004322
##  [36] 0.3063924 0.5534378 0.5811441 0.5097919 0.2468361 0.3782586 0.6754070
##  [43] 0.3251770 0.3885657 0.4866480 0.5461808 0.3388354 0.4972481 0.6069386
##  [50] 0.1734082 0.5841854 0.1825670 0.7321562 0.7072252 0.2068596 0.4986379
##  [57] 0.3789680 0.3753707 0.4223732 0.4829185 0.5166820 0.6233485 0.4267682
##  [64] 0.2508832 0.2424451 0.4546804 0.5686296 0.5338262 0.3117600 0.1285049
##  [71] 0.5537374 0.3859767 0.3467693 0.3102380 0.4984951 0.3995398 0.5189746
##  [78] 0.6514056 0.4870990 0.5505915 0.6088603 0.5588382 0.2577329 0.4148124
##  [85] 0.6661753 0.7100346 0.6420909 0.5625488 0.6653349 0.6656333 0.4754505
##  [92] 0.4785773 0.5716569 0.5399865 0.4206195 0.4603906 0.7398491 0.6044728
##  [99] 0.3499299 0.3627396

Ambil rata-rata untuk masing-masing baris dalam sample tersebut

meanSampSim <- list(NA)
meanSampPop1 <- list(NA)
meanSampPop2 <- list(NA)
meanSampPop3 <- list(NA)

for (i in 1:length(n)) {
    meanSampSim[[i]] <- apply(sampSim[[i]], 1, mean)
    meanSampPop1[[i]] <- apply(sampPop1[[i]], 1, mean)
    meanSampPop2[[i]] <- apply(sampPop2[[i]], 1, mean)
    meanSampPop3[[i]] <- apply(sampPop3[[i]], 1, mean)
}

str(meanSampSim)

## List of 5
##  $ : num [1:100] 0.146 0.441 0.506 0.547 0.525 ...
##  $ : num [1:100] 0.64 0.507 0.615 0.501 0.485 ...
##  $ : num [1:100] 0.472 0.435 0.526 0.516 0.519 ...
##  $ : num [1:100] 0.487 0.44 0.485 0.487 0.501 ...
##  $ : num [1:100] 0.511 0.496 0.509 0.492 0.497 ...

Buat histogram untuk masing-masing sebaran

par(mfrow = c(2,3))

# Untuk sebarn simetris
hist(meanSampSim[[1]], breaks = 12, 
     main = "Rataan Sebaran Simetris untuk n = 4", xlab = "")
hist(meanSampSim[[2]], breaks = 12, 
     main = "Rataan Sebaran Simetris untuk n = 12", xlab = "")
hist(meanSampSim[[3]], breaks = 12, 
     main = "Rataan Sebaran Simetris untuk n = 20", xlab = "")
hist(meanSampSim[[4]], breaks = 12, 
     main = "Rataan Sebaran Simetris untuk n = 60", xlab = "")
hist(meanSampSim[[5]], breaks = 12, 
     main = "Rataan Sebaran Simetris untuk n = 100", xlab = "")

Pada sebaran simetris, sebaran rataan sampel terlihat mulai simetris dan menyebar normal dari n=20.

par(mfrow = c(2,3))

# Untuk sebarn campuran 1
hist(meanSampPop1[[1]], breaks = 12, 
     main = "Rataan Sebaran Campuran 1 untuk n = 4", xlab = "")
hist(meanSampPop1[[2]], breaks = 12, 
     main = "Rataan Sebaran Campuran 1 untuk n = 12", xlab = "")
hist(meanSampPop1[[3]], breaks = 12, 
     main = "Rataan Sebaran Campuran 1 untuk n = 20", xlab = "")
hist(meanSampPop1[[4]], breaks = 12, 
     main = "Rataan Sebaran Campuran 1 untuk n = 60", xlab = "")
hist(meanSampPop1[[5]], breaks = 12, 
     main = "Rataan Sebaran Campuran 1 untuk n = 100", xlab = "")

Pada sebaran campuran normal baku dan chisquare (1), sebaran rataan sampel terlihat mulai simetris dan menyebar normal dari n=12.

par(mfrow = c(2,3))

# Untuk sebarn campuran 1
hist(meanSampPop2[[1]], breaks = 12, 
     main = "Rataan Sebaran Campuran 2 untuk n = 4", xlab = "")
hist(meanSampPop2[[2]], breaks = 12, 
     main = "Rataan Sebaran Campuran 2 untuk n = 12", xlab = "")
hist(meanSampPop2[[3]], breaks = 12, 
     main = "Rataan Sebaran Campuran 2 untuk n = 20", xlab = "")
hist(meanSampPop2[[4]], breaks = 12, 
     main = "Rataan Sebaran Campuran 2 untuk n = 60", xlab = "")
hist(meanSampPop2[[5]], breaks = 12, 
     main = "Rataan Sebaran Campuran 2 untuk n = 100", xlab = "")

Pada sebaran campuran chisquare (1) dan chisquare (2), sebaran rataan sampel terlihat mulai simetris dan menyebar normal dari n=100, meskipun pada n=100 masih terlihat menjulur ke kanan.

par(mfrow = c(2,3))

# Untuk sebarn campuran 1
hist(meanSampPop3[[1]], breaks = 12, 
     main = "Rataan Sebaran Campuran 3 untuk n = 4", xlab = "")
hist(meanSampPop3[[2]], breaks = 12, 
     main = "Rataan Sebaran Campuran 3 untuk n = 12", xlab = "")
hist(meanSampPop3[[3]], breaks = 12, 
     main = "Rataan Sebaran Campuran 3 untuk n = 20", xlab = "")
hist(meanSampPop3[[4]], breaks = 12, 
     main = "Rataan Sebaran Campuran 3 untuk n = 60", xlab = "")
hist(meanSampPop3[[5]], breaks = 12, 
     main = "Rataan Sebaran Campuran 3 untuk n = 100", xlab = "")

Pada sebaran campuran chisquare (1), chisquare (2), normal (0,1), dan normal (1, 2) sebaran rataan sampel terlihat mulai simetris dan menyebar normal dari n=12 dan merupakan sebaran rataan data yang paling cepat menyerupai sebaran normal.

Diskusi 2

Latihan 2

Seorang peneliti ingin melihat bagaimana nilai koefisien korelasi dari metode korelasi pearson dan metode korelasi spearman untuk satu populasi pada berbagai ukuran sampel jika datanya terdapat pencilan. Bangkitan data pencilan dengan cara berikut: a.Bangkitkan data berdasarkan distribusi multivariate Normal dengan vektor rataan (10,12) dan korelasi 0.85 serta vektor simpangan baku (2,3) b.Bangkitan outlier pada kedua peubah dengan berdasarkan distribusi N(25,2) sebanyak 20% dari ukuran sampel pada langkah a. c.Ganti amatan dari hasil bangkitan pada langkah a dengan hasil bangkitan pada langkah b secara acak. Bangkitkan data dengan ulangan 1000 dan ukuran sampel yang diteliti adalah 10,20,30,40,50,60,70,80,90,100. Berikan kesimpulan dari hasil simulasi!

Jawaban 2

Langkah 1 abstraksi bangkitan data langkah a

library(MBESS)
library(MASS)
library(mvtnorm)
library(ggplot2)
set.seed(123)
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# membangkitkan data berdistribus multivariate normal
x_mvn <- rmvnorm(n=40, mean=c(10,12), sigma=cov_mat)

langkah b

set.seed(123)
#membangkitkan data berdasarkan distribusi normal
x1_out <- rnorm(8,17,2)
y1_out <- rnorm(8,17,2)
x1y1_out <-cbind(x1_out,y1_out)

langkah c

# memilih secara acak amatan yang akan diganti dengan pencilan
index_out <- sample.int(nrow(x_mvn),8)
# mengganti amatan dari langkah a dengan amatan dari langkah b
x_mvn[index_out,] <- x1y1_out

qplot(x = V1,y=V2,data = as.data.frame(x_mvn))

## Warning: `qplot()` was deprecated in ggplot2 3.4.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Langkah 2 ada mendefinisikan parameter-parameter yang digunakan untuk simulasi

# mendefinisikan ukuran sampel yang dikaji
n <- seq(10,100,10)
# mendefinisikan jumlah ulangan
num_rep2 <- 1000
# mendefinisikan vektor simpangan baku
sd0 <- c(2,3)
# mendefinisikan matrix korelasi
cor_mat <- matrix(c(1,0.85,0.85,1),nrow=2)
# mengubah matrix korelasi ke kovarians
cov_mat <- MBESS::cor2cov(cor_mat,sd0)
# mendefinisikan vektor rataan
mu_mvn <- c(10,12)
# mendefinisikan mu dan sigma untuk data pencilan
mu1 <- 25
sigma1 <- 2

Langkah 3 mendefinisikan fungsi-fungsi yang dibutuhkan

# generate multivariate normal dengan pencilan
gen_mvn_out <- function(n0,mu0,sigma0,n1,mu1,sigma1){
x_mvn <- rmvnorm(n=n0, mean=mu0, sigma=sigma0)
x1_out <- rnorm(n1,mu1,sigma1)
y1_out <- rnorm(n1,mu1,sigma1)
x1y1_out <-cbind(x1_out,y1_out)
index_out <- sample.int(nrow(x_mvn),n1)
x_mvn[index_out,] <- x1y1_out
return(x_mvn)
}

# extract koefisien korelasi
cor_res <- function(x_mat,method){
  cor(x_mat[,1],x_mat[,2],method = method)
}

Langkah 4 Menjalankan simulasi

set.seed(123)
# membuat objek data.frame dummy untuk menyimpan hasil simulasi
result_sim2 <- data.frame(n=0,correlation=0,type="")
for(i in seq_along(n)){
# membangkitan data berdistribusi multivariate normal dengan ulangan 1rb
sim_mc_mvn <- replicate(num_rep2,                     gen_mvn_out(n0 = n[i],mu_mvn,cov_mat,n1=0.2*n[i],mu1,sigma1),
                        simplify = F)
# menghitung keofisien korelasi pearson dari 1rb data hasil bangkitan
pearson_res <- sapply(sim_mc_mvn,cor_res,method="pearson")
# menghitung keofisien spearman pearson dari 1rb data hasil bangkitan
spearman_res <- sapply(sim_mc_mvn,cor_res,method="spearman")
# menghitung rata-rata koefisien
pearson_sim <- mean(pearson_res)
spearman_sim <- mean(spearman_res)
# menyimpan hasil simulasi
result_sim2 <- rbind(result_sim2,list(n[i],pearson_sim,"pearson"),
                  list(n[i],spearman_sim,"spearman")
)
}
result_sim2 <- result_sim2[-1,]

Langkah 5 Membuat Grafik berdsarkan hasil simulasi

ggplot(data = result_sim2,aes(n,correlation))+
  geom_point(aes(color=type))+geom_hline(yintercept = 0.85)+scale_y_continuous(limits = c(0.5,1))

Berdasarkan perbandingan kinerja antara korelasi Pearson dan Spearman di atas, terlihat bahwa metode korelasi Spearman lebih unggul dalam menangani data yang mengandung outlier. Hal ini dapat dilihat dari nilai korelasi Spearman yang mendekati 0.85, yaitu nilai korelasi yang telah ditetapkan. Dengan demikian, dapat disimpulkan bahwa metode korelasi Spearman lebih cocok digunakan dalam situasi di mana terdapat outlier pada data yang dianalisis.

KEL 11 MSR

Kelompok 11 P2

2023-05-08