masaha03

n <- 1000000 # ケース数
df <- data.frame(v1=runif(n),
                 v2=runif(n),
                 v3=runif(n),
                 v4=runif(n),
                 v5=runif(n),
                 v6=runif(n))
gc()
system.time(
lm(v1 ~ v2 + v3 + v4 + v5 + v6, data=df)

masaha03

# 2012/12/14
# 実装をシンプルに
# vcdパッケージがあればCramer's Vを出力
# summary.CrossTableにlatexオプション（ただし表組み部分のみ）

# 2012/05/15 罫線を修正

crossTable <- function(..., deparse.level = 2){
  table <- table(..., deparse.level = deparse.level)
  class(table) <- c("CrossTable", "table")

masaha03

# endog.variancesは足りない誤差項を自動的に追加するオプション
model.b0 <- specifyEquations(endog.variances=TRUE)
教育年数 = b1*年齢
職業威信 = c1*年齢 + c2*教育年数

# モデルの結合
model.b2 <- combineModels(model.a2, model.b0)
# 重複しているパスの削除
model.b2 <- removeRedundantPaths(model.b2)

masaha03

data(iris)

model.iris <- cfa()
F1: Sepal.Length, Sepal.Width, Petal.Length, Petal.Width

# 特殊記号はエスケープする必要がある
model.iris2 <- update(model.iris)
replace, lam\[Petal.Width:F1\], b2

model.iris3 <- update(model.iris)

masaha03

cor.3 <- read.moments(diag=FALSE,
names=c("年齢", "個人年収", "教育年数", "職業威信",
        "両親に服従", "権威に敬意", "指導者に頼る", "伝統重視"))
 0.14196474
-0.26729787  0.22629459
-0.02846850  0.39823207  0.27631486
 0.16409535  0.04589625 -0.09924375 -0.04642942
 0.05328158 -0.03399866 -0.04022705 -0.06889775  0.48134221
 0.18186547 -0.06105599 -0.09025803 -0.09414690  0.47089471  0.53716791
 0.16529548 -0.03123394 -0.07494290 -0.10559612  0.39698678  0.42569827  0.44360292

masaha03

cor.2 <- readMoments(diag=FALSE,
    names=c("両親に服従", "権威に敬意", "指導者に頼る", "伝統重視"))
0.5063549
0.4457519 0.5200434
0.3975495 0.4201535 0.4834760

# 従来の記法
model.c1 <- specifyModel()
権威主義  -> 両親に服従   , NA, 1
権威主義  -> 権威に敬意   , a2, NA

masaha03

model.b1 <- update(model.a2)
add, 教育年数 <-  年齢    , b1, NA
add, 職業威信 <-  年齢    , c1, NA
add, 職業威信 <-  教育年数, c2, NA
add, 教育年数 <-> 教育年数, e2, NA
add, 職業威信 <-> 職業威信, e3, NA

sem.b <- sem(model.b1, cor.1, N=685, fixed.x="年齢")
summary(sem.b)

masaha03

library(sem)

cor.1 <- readMoments(diag=FALSE,
    names=c("年齢", "個人年収", "教育年数", "職業威信"))
 0.12107847
-0.28075478 0.2329769
-0.04963399 0.3954726  0.2921118

# 従来の記法
model.a1 <- specifyModel()

masaha03

N.lv <- 5
x <- factor(rep(LETTERS[1:N.lv], each = 20))
y <- rnorm(100)
grp <- factor(rep(c("G1", "G2"), 50)) # Group 1 and 2

library(ggplot2)
qplot(x, y, geom="boxplot", fill=grp)

masaha03

f <- function(i){
  l <- c("Fizz", "Buzz", "FizzBuzz")
  ii <- (0==i%%5)*2 + (0==i%%3)*1
  print(l[ii])
  print(rep(l[ii], 10)[i])
  return(ifelse(0<ii, l[ii], i))
}