## Example: Percentage of body fat and body measurements ------------- ## Authors: Cesar Augusto Taconeli and Wagner Hugo Bonat LEG/UFPR ---- ## Date: 07/05/2019 -------------------------------------------------- rm(list=ls()) ## Loading extra packages require(bbmle) require(faraway) ## Loading extra functions source("Taconeli_Bonat.R") ## Loading data set data(fat) head(fat) data_full <- fat[which(fat$brozek != 0),] par(mfrow = c(1,3)) hist(data_full$brozek) with(data_full, plot(brozek ~ neck)) with(data_full, plot(brozek ~ adipos)) ## Sampling schemes seed <- 30 set.seed(seed) m <- 6 n <- 5 N <- m*n ## Case 1: weak correlation ------------------------------------------ ## Ranked set sampling (RSS) data <- data_full[, c('brozek', 'neck')] sample_rss <- rss_sampling_imp(data = data, m = m, n = n) # Simple random sample (SRS) sample_srs <- sample(data$brozek, m*n, replace = TRUE) ## Initial values p10 <- log(3) p20 <- log(20) ## Fitting whole data set (population) pop_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = data$brozek), optimizer="optim") param_srs <- exp(coef(pop_srs)) ## Fitting sample data ----------------------------------------------- ## MLE-RSS mle_rss <- mle2(LikBisaRSS, start = list(p1 = p10, p2 = p20), data = list(y = sample_rss, n = n), optimizer="optim") est_mle_rss <- exp(coef(mle_rss)) ## MLE-SRS mle_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = sample_srs, n = n), optimizer="optim") est_mle_srs <- exp(coef(mle_srs)) ## MPS y <- sort(sample_rss[,2]) mps <- mle2(f_mps, start = list(p1 = p10, p2 = p20), data = list(y = y, N = N), optimizer="optim") est_mps <- exp(coef(mps)) ## OLS ind <- 1:(m*n) ols <- mle2(f_ols, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_ols <- exp(coef(ols)) ## WLS wls <- mle2(f_wls, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_wls <- exp(coef(wls)) ## CVM cvm <- mle2(f_cvm, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_cvm <- exp(coef(cvm)) ## AND and <- mle2(f_and, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_and <- exp(coef(and)) ## RAND rand <- mle2(f_rand, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_rand <- exp(coef(rand)) ## Comparing methods weak <- rbind(my_gof(y = data[,1], p1 = est_mle_rss[1], p2 = est_mle_rss[2]), my_gof(y = data[,1], p1 = est_mps[1], p2 = est_mps[2]), my_gof(y = data[,1], p1 = est_ols[1], p2 = est_ols[2]), my_gof(y = data[,1], p1 = est_wls[1], p2 = est_wls[2]), my_gof(y = data[,1], p1 = est_cvm[1], p2 = est_cvm[2]), my_gof(y = data[,1], p1 = est_and[1], p2 = est_and[2]), my_gof(y = data[,1], p1 = est_rand[1], p2 = est_rand[2]), my_gof(y = data[,1], p1 = est_mle_srs[1], p2 = est_mle_srs[2]) ) rownames(weak) <- c("MLE_RSS","MPS","OLS","WLS","CVM","AND","RAND","MLE_SRS") round(weak,4) ## Case 2: strong correlation ---------------------------------------- rm(data) ## Ranked set sampling (RSS) data <- data_full[, c('brozek', 'adipos')] set.seed(seed) sample_rss <- rss_sampling_imp(data = data, m = m, n = n) # Simple random sample (SRS) sample_srs <- sample(data$brozek, m*n, replace = TRUE) ## Initial values p10 <- log(3) p20 <- log(20) ## Fitting whole data set (population) pop_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = data$brozek), optimizer="optim") param_srs <- exp(coef(pop_srs)) ## Fitting sample data ----------------------------------------------- ## MLE-RSS mle_rss <- mle2(LikBisaRSS, start = list(p1 = p10, p2 = p20), data = list(y = sample_rss, n = n), optimizer="optim") est_mle_rss <- exp(coef(mle_rss)) ## MLE-SRS mle_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = sample_srs, n = n), optimizer="optim") est_mle_srs <- exp(coef(mle_srs)) ## MPS y <- sort(sample_rss[,2]) mps <- mle2(f_mps, start = list(p1 = p10, p2 = p20), data = list(y = y, N = N), optimizer="optim") est_mps <- exp(coef(mps)) ## OLS ind <- 1:(m*n) ols <- mle2(f_ols, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_ols <- exp(coef(ols)) ## WLS wls <- mle2(f_wls, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_wls <- exp(coef(wls)) ## CVM cvm <- mle2(f_cvm, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_cvm <- exp(coef(cvm)) ## AND and <- mle2(f_and, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_and <- exp(coef(and)) ## RAND rand <- mle2(f_rand, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_rand <- exp(coef(rand)) ## Comparing methods strong <- rbind(my_gof(y = data[,1], p1 = est_mle_rss[1], p2 = est_mle_rss[2]), my_gof(y = data[,1], p1 = est_mps[1], p2 = est_mps[2]), my_gof(y = data[,1], p1 = est_ols[1], p2 = est_ols[2]), my_gof(y = data[,1], p1 = est_wls[1], p2 = est_wls[2]), my_gof(y = data[,1], p1 = est_cvm[1], p2 = est_cvm[2]), my_gof(y = data[,1], p1 = est_and[1], p2 = est_and[2]), my_gof(y = data[,1], p1 = est_rand[1], p2 = est_rand[2]), my_gof(y = data[,1], p1 = est_mle_srs[1], p2 = est_mle_srs[2]) ) rownames(strong) <- c("MLE_RSS","MPS","OLS","WLS","CVM","AND","RAND","MLE_SRS") round(strong,4) ## Case 2: Perfectly ranking ----------------------------------------- rm(data) ## Ranked set sampling (RSS) data <- data_full[, c('brozek', 'brozek')] set.seed(seed) sample_rss <- rss_sampling_imp(data = data, m = m, n = n) # Simple random sample (SRS) sample_srs <- sample(data$brozek, m*n, replace = TRUE) ## Initial values p10 <- log(3) p20 <- log(20) ## Fitting whole data set (population) pop_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = data$brozek), optimizer="optim") param_srs <- exp(coef(pop_srs)) ## Fitting sample data ----------------------------------------------- ## MLE-RSS mle_rss <- mle2(LikBisaRSS, start = list(p1 = p10, p2 = p20), data = list(y = sample_rss, n = n), optimizer="optim") est_mle_rss <- exp(coef(mle_rss)) ## MLE-SRS mle_srs <- mle2(LikBisaSRS, start = list(p1 = p10, p2 = p20), data = list(y = sample_srs, n = n), optimizer="optim") est_mle_srs <- exp(coef(mle_srs)) ## MPS y <- sort(sample_rss[,2]) mps <- mle2(f_mps, start = list(p1 = p10, p2 = p20), data = list(y = y, N = N), optimizer="optim") est_mps <- exp(coef(mps)) ## OLS ind <- 1:(m*n) ols <- mle2(f_ols, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_ols <- exp(coef(ols)) ## WLS wls <- mle2(f_wls, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_wls <- exp(coef(wls)) ## CVM cvm <- mle2(f_cvm, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_cvm <- exp(coef(cvm)) ## AND and <- mle2(f_and, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_and <- exp(coef(and)) ## RAND rand <- mle2(f_rand, start = list(p1 = p10, p2 = p20), data = list(y = y, i = ind, N = N), optimizer="optim") est_rand <- exp(coef(rand)) ## Comparing methods perf <- rbind(my_gof(y = data[,1], p1 = est_mle_rss[1], p2 = est_mle_rss[2]), my_gof(y = data[,1], p1 = est_mps[1], p2 = est_mps[2]), my_gof(y = data[,1], p1 = est_ols[1], p2 = est_ols[2]), my_gof(y = data[,1], p1 = est_wls[1], p2 = est_wls[2]), my_gof(y = data[,1], p1 = est_cvm[1], p2 = est_cvm[2]), my_gof(y = data[,1], p1 = est_and[1], p2 = est_and[2]), my_gof(y = data[,1], p1 = est_rand[1], p2 = est_rand[2]), my_gof(y = data[,1], p1 = est_mle_srs[1], p2 = est_mle_srs[2]) ) rownames(perf) <- c("MLE_RSS","MPS","OLS","WLS","CVM","AND","RAND","MLE_SRS") round(perf,4) round(cbind(weak, strong, perf),4) # END ----------------------------------------------------------------