An R function: OLS/Robust scaled mass index

R code

	# Inspired by Peig & Green
	# "New perspectives for estimating body condition from mass/length data:
	# the scaled mass index as an alternative method"
	# Oikos 118: 1883-1891, 2009
	# Author: Chen-Pan Liao
	scaledMassIndex <-
	function(x, y, x.0 = mean(x)) {
	require(smatr)
	require(magrittr)
	require(MASS)
	require(data.table)
	logM.ols <- lm(log(y) ~ log(x))
	logM.rob <- rlm(log(y) ~ log(x), method = "M")
	b.msa.ols <- coef(sma(log(y) ~ log(x)))[2]
	b.msa.rob <- coef(sma(log(y) ~ log(x), robust = T))[2]
	SMI.ols <- y * (x.0 / x) ^ b.msa.ols
	SMI.rob <- y * (x.0 / x) ^ b.msa.rob
	res <- data.frame(SMI.ols, SMI.rob, x, y)
	pred.DT <-
	data.table(x = seq(min(x), max(x), length = 100)) %>%
	.[, y.ols := predict(logM.ols, newdata = .) %>% exp] %>%
	.[, y.rob := predict(logM.rob, newdata = .) %>% exp]
	attr(res, "b.msa") <- c(ols = b.msa.ols, rob = b.msa.rob)
	return(res)
	}

view raw robust-scaledMassIndex.R hosted with ❤ by GitHub

R 繪圖：柱狀圖的低階繪圖範例

R 繪圖時的柱狀圖如何加上誤差線？如何使用線條填色？如何把圖示放在外部？本文提供一個例子。

以下是產生上圖的原始碼。

	# dat 為輸入的二因子資料
	set.seed(1234)
	dat <- data.frame(
	y = rnorm(80, 5, 2),
	x1 = gl(2, 40, labels = c("A","B")),
	x2 = gl(4, 10, labels = c("1","2","3","4"))
	)
	# 利用 tapply 求出各組之 mean 和 sd
	y.mean <- with(dat, tapply(y, paste(x1, x2, sep="."), mean))
	y.sd <- with(dat, tapply(y, paste(x1, x2, sep="."), sd))
	# 為了方便之後畫圖，把 y.mean 和 y.sd 轉成矩陣
	# 矩陣各列表示 X2，各欄表示 X1
	Y <- matrix(y.mean, 4)
	E <- matrix(y.sd, 4)
	colnames(Y) <- colnames(E) <- c("A", "B")
	rownames(Y) <- rownames(E) <- c("1", "2", "3", "4")
	# 接下來就開始畫圖
	# 過程小小複雜，若有不了解可以逐步進行並觀看圖怎麼變化
	# 先開一個畫布
	# 因為我想把圖示畫在右側外部，所以 mar 和 xpd 要自行設定
	dev.new(width = 5, height = 5)
	par(mar=c(4, 4, 1, 10) + 0.1, xpd = T)
	# 接下來是畫上柱狀圖主體
	# 因為 X2 == "4" 柱體要畫成同時有「斜線填色」和「顏色填色」
	# 所以 barplot 要畫二次
	# 以下是第一次，只畫主要文字、框線和顏色填色
	# 我習慣先自行算好 ylim
	# 注意 barplot 的回傳值被我定義為變數 bp（之後有用）
	bp <- barplot(
	Y, beside = T,
	xaxt = "n", ylim = c(0,10),
	ylab = "Title for y-axis (dat$y)", xlab = "Title for x-axis (dat$x1)",
	col = c(0,3,0,"yellow")
	)
	# 接下來是疊上 X2 == "4" 柱體的「斜線填色」
	# 上一步已經畫過的就不要再畫一次
	# 只要畫上 X2 == "4" 柱體的「斜線填色」即可
	# 注意到 add = T 就是疊圖的作用
	barplot(
	Y, beside = T,
	xaxt = "n", ylim = c(0,10), yaxt = "n",
	ylab = "", xlab = "",
	col = c(0,1,1,1), density = c(0,0,25,15), angle = c(NULL,NULL,45,135),
	border = c(0,0,0,0),
	add = T
	)
	# 這裡只是示範怎麼用 axis()
	axis(1, colMeans(bp), c("A", "B"), tick = F)
	# 這裡是畫「只向上方延伸的 error bar」
	# 這時候 bp 的值（也就是每個柱的橫軸位置）就很有用了
	arrows(bp, Y, bp, Y+E, length = 0.1, angle = 90)
	# 當然也可以再補畫上「向下方延伸的 error bar」
	# 但這個例子不太適當就不做了
	# 有需要的話請再加入
	# arrows(bp, Y, bp, Y-E, length = 0.1, angle = 90)
	# 接下來是畫圖示
	# 這裡示範如何把圖示畫在右側外部
	# 需要與前面 par() 配合才行
	# 也因為 X2 == "4" 柱體的同時有「斜線填色」和「顏色填色」
	# 所以 legend() 一樣要畫兩次
	# 第一次負責文字、線框、顏色填色
	legend(
	"right",
	legend = c("1", "2", "3", "4"),
	fill = c(0,3,0,"yellow"),
	title = "Row name (dat$x2)",
	text.width = 5.5,
	box.lwd = 1, col = c(1,1,1,1),
	inset=c(-0.9,0)
	)
	# 第二次負責 X2 == "4" 的斜線填色就好
	# 重覆的東西都不要再畫
	# 注意兩次畫 legend 要完全重疊才會成功
	# 因此我定了 text.width 之類看似多餘的參數
	legend(
	"right",
	legend = c("", "", "", ""),
	fill = c(0,0,1,1),
	border = c(0,0,0,0),
	density = c(0,0,25,15), angle = c(NULL,NULL,45,135),
	title = "",
	text.width = 5.5,
	box.lwd = 0, col = c(0,0,0,0),
	inset=c(-0.9,0)
	)
	# 畫個框框
	box()

view raw R-bar-plot-example.R hosted with ❤ by GitHub

「在 95% 信心水準，抽樣誤差於正負 3.1% 以內」到底是什麼意思？

如果你對統計學有點興趣，可以往下看看這篇文。不知不覺寫得有點長。

故事：「某民調訪問 1000 名受訪者，有 40% 比例的人支持；在 95% 信心水準，抽樣誤差於正負 3.1% 以內。」現實中，民調常常出現相似的句子。我相信很多人看不懂這句話到底是什麼意思，或是誤解了這句話。所以現在來聊聊這到底是什麼鬼好了。

我們先把上面這個故事改成另一種情境好了。

例子：「你有一個大袋子。袋子裡有無數顆球。球只可能是黑色或白色。你希望知道一件事：到底黑球在大袋子中的個數比例是多少。現在，你抽出 1000 顆球，然後你發現有 400 顆是黑球（40%）。在 95% 信心水準，抽樣誤差於正負 3.1% 以內，也就是黑球比例的區間估計是 40% − 3.1% 到 40% + 3.1%。」

如果你可以看出這個例子其實和民調結果的故事是同一件事，那就可以往下看了。

在例子中，每一次調查就是抽出 n 顆球（例子中的 n = 1000）。袋子中黑球的比例，也就是我們希望知道的未知數，我們叫 p 好了。因為球實在太多了，我們不可能知道 p 的大小，只有抽出 n 顆球來估 p 的大小。

在統計學理，我們可以在還沒抽出球之前就預知「每次調查有多準」。什麼叫「有多準」呢？統計學的辦法是決定一個機率和設計一個區間，並估計 p 有多大的機率會座落在這個設計出來的區間之內。如果你每次調查抽出的球數 n 很大很大，那這個區間會很容易計算。

有點難懂？我試著換個方式說。例如，如果我可以做 100 次相同調查，我將會得到 100 個黑球比例和 100 個區間。我可以保證，這 100 次調查中，約有 95 個區間會包括了那個我們不可能知道的 p。這就是「在 95% 信心水準」的實際意義。

你也許會問，可不可能算 100 個區間有 99 個區間會包括 p 的調查？可以，但每個區間大小會變寬（由 ± 3.1% 變成 ± 4.1%），但結果也變得更值得信賴（因為區間將有 99% 的機率可以包括了真實的 p）。

如果你每次抽出更多顆球，那每個區間也都會變窄（因為抽出越多球就準嘛）。例如，每次抽出 5000 顆球，那這個區間大小會由 ± 3.1% 變成 ± 1.4%。

回到原本的例子。如果以相同的調查方法（一樣訪問 1000 名受訪者的情況下），並且重覆很多很多次（當然，實際上只進行了一次）。我們可以保證，那個我們想估計的支持度將座落在 95% 次所算出的區間之內。在故事中，這次調查結果的區間是 40% − 3.1% 到 40% + 3.1%。我們 *永遠不可能* 知道這次得到的區間會不會包括了真實的支持度 p，但它很可能會包括，因為這個調查方法就保證了進行 100 次約有 95 次是成立的。

有些人會把 40% − 3.1% 到 40% + 3.1% 理解成因為調查進行所造成的誤差，例如打錯資料、拒絕受訪、無效結果之類的。這都不是信賴區間的意義。

我並不想把這中間的數學式寫出來啦。只想說個概念。