options(digits = 4, show.signif.stars = FALSE)載入package MASS並檢視cats中前六筆資料。
library(MASS)
head(dta <- cats)## Sex Bwt Hwt
## 1 F 2.0 7.0
## 2 F 2.0 7.4
## 3 F 2.0 9.5
## 4 F 2.1 7.2
## 5 F 2.1 7.3
## 6 F 2.1 7.6採dummy coding方式將性別類別變項轉為虛擬變數。
coef(lm(Bwt ~ Sex, data = cats))## (Intercept) SexM
## 2.3596 0.5404截距為對照組(Sex=F)的平均值，迴歸係數為實驗組(Sex=M)之平均數減截距。
採effect coding將性別類別變項轉為虛擬變數。
coef(lm(Bwt ~ Sex, data = cats, contrasts = list(Sex = "contr.sum")))## (Intercept) Sex1
## 2.6298 -0.2702截距為未加權總平均數，迴歸係數為截距減實驗組(Sex=M)之細格平均數。
載入並檢視資料。
dta <- read.table("readingtimes.txt", header = T)
head(dta)## Snt Sp Wrds New S01 S02 S03 S04 S05 S06 S07 S08 S09
## 1 1 1 13 1 3.429 2.795 4.161 3.071 3.625 3.161 3.232 7.161 1.536
## 2 2 2 16 3 6.482 5.411 4.491 5.063 9.295 5.643 8.357 4.313 2.946
## 3 3 3 9 2 1.714 2.339 3.018 2.464 6.045 2.455 4.920 3.366 1.375
## 4 4 4 9 2 3.679 3.714 2.866 2.732 4.205 6.241 3.723 6.330 1.152
## 5 5 5 10 3 4.000 2.902 2.991 2.670 3.884 3.223 3.143 6.143 2.759
## 6 6 6 18 4 6.973 8.018 6.625 7.571 8.795 13.188 11.170 6.071 7.964
## S10
## 1 4.063
## 2 6.652
## 3 2.179
## 4 3.661
## 5 3.330
## 6 7.866rank(colMeans(dta[,5:14]))## S01 S02 S03 S04 S05 S06 S07 S08 S09 S10
## 5 4 2 3 10 9 8 7 1 6colSums(dta[,5:14])/sum(dta$Wrds)## S01 S02 S03 S04 S05 S06 S07 S08 S09 S10
## 0.3569 0.3325 0.3262 0.3266 0.4857 0.4642 0.4518 0.4331 0.2369 0.3883getAnywhere(t.test.default)## A single object matching 't.test.default' was found
## It was found in the following places
## registered S3 method for t.test from namespace stats
## namespace:stats
## with value
##
## function (x, y = NULL, alternative = c("two.sided", "less", "greater"),
## mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95,
## ...)
## {
## alternative <- match.arg(alternative)
## if (!missing(mu) && (length(mu) != 1 || is.na(mu)))
## stop("'mu' must be a single number")
## if (!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) ||
## conf.level < 0 || conf.level > 1))
## stop("'conf.level' must be a single number between 0 and 1")
## if (!is.null(y)) {
## dname <- paste(deparse(substitute(x)), "and", deparse(substitute(y)))
## if (paired)
## xok <- yok <- complete.cases(x, y)
## else {
## yok <- !is.na(y)
## xok <- !is.na(x)
## }
## y <- y[yok]
## }
## else {
## dname <- deparse(substitute(x))
## if (paired)
## stop("'y' is missing for paired test")
## xok <- !is.na(x)
## yok <- NULL
## }
## x <- x[xok]
## if (paired) {
## x <- x - y
## y <- NULL
## }
## nx <- length(x)
## mx <- mean(x)
## vx <- var(x)
## if (is.null(y)) {
## if (nx < 2)
## stop("not enough 'x' observations")
## df <- nx - 1
## stderr <- sqrt(vx/nx)
## if (stderr < 10 * .Machine$double.eps * abs(mx))
## stop("data are essentially constant")
## tstat <- (mx - mu)/stderr
## method <- if (paired)
## "Paired t-test"
## else "One Sample t-test"
## estimate <- setNames(mx, if (paired)
## "mean of the differences"
## else "mean of x")
## }
## else {
## ny <- length(y)
## if (nx < 1 || (!var.equal && nx < 2))
## stop("not enough 'x' observations")
## if (ny < 1 || (!var.equal && ny < 2))
## stop("not enough 'y' observations")
## if (var.equal && nx + ny < 3)
## stop("not enough observations")
## my <- mean(y)
## vy <- var(y)
## method <- paste(if (!var.equal)
## "Welch", "Two Sample t-test")
## estimate <- c(mx, my)
## names(estimate) <- c("mean of x", "mean of y")
## if (var.equal) {
## df <- nx + ny - 2
## v <- 0
## if (nx > 1)
## v <- v + (nx - 1) * vx
## if (ny > 1)
## v <- v + (ny - 1) * vy
## v <- v/df
## stderr <- sqrt(v * (1/nx + 1/ny))
## }
## else {
## stderrx <- sqrt(vx/nx)
## stderry <- sqrt(vy/ny)
## stderr <- sqrt(stderrx^2 + stderry^2)
## df <- stderr^4/(stderrx^4/(nx - 1) + stderry^4/(ny -
## 1))
## }
## if (stderr < 10 * .Machine$double.eps * max(abs(mx),
## abs(my)))
## stop("data are essentially constant")
## tstat <- (mx - my - mu)/stderr
## }
## if (alternative == "less") {
## pval <- pt(tstat, df)
## cint <- c(-Inf, tstat + qt(conf.level, df))
## }
## else if (alternative == "greater") {
## pval <- pt(tstat, df, lower.tail = FALSE)
## cint <- c(tstat - qt(conf.level, df), Inf)
## }
## else {
## pval <- 2 * pt(-abs(tstat), df)
## alpha <- 1 - conf.level
## cint <- qt(1 - alpha/2, df)
## cint <- tstat + c(-cint, cint)
## }
## cint <- mu + cint * stderr
## names(tstat) <- "t"
## names(df) <- "df"
## names(mu) <- if (paired || !is.null(y))
## "difference in means"
## else "mean"
## attr(cint, "conf.level") <- conf.level
## rval <- list(statistic = tstat, parameter = df, p.value = pval,
## conf.int = cint, estimate = estimate, null.value = mu,
## alternative = alternative, method = method, data.name = dname)
## class(rval) <- "htest"
## return(rval)
## }
## <bytecode: 0x0000000009175038>
## <environment: namespace:stats>