R: glmnetでlassoをやってみたメモ [統計]
Rのglmnetパッケージでlassoをつかってみたのでメモです。
昔かいた論文のネタから、京都市内で採取したスギの年輪幅データと、当年および前年の月平均気温・月降水量との関係の関係を解析してみます。スギの年輪はデータは19本分あって時系列データなので，状態空間モデルにして個体ランダム効果をつけたいところですが、今回は、個体ごとに直線回帰からの残差を標準化して、それを平均したものをインデックスとして使用しました。
今回のコードは下記のとおりです。データはこのあたりにおいておきます。→ Nenrin.csv, Tenki.csv
library(ggplot2)
library(reshape2)
library(glmnet)
nenrin <- read.csv("Nenrin.csv")
tenki <- read.csv("Tenki.csv", colClasses = rep("numeric", 4))
tenki$YM <- tenki$Year + (tenki$Month - 1) / 12
p <- ggplot(tenki)
p + geom_line(aes(x = YM, y = Temperture)) + xlab("Year")
p + geom_bar(aes(x = YM, y = Precipitation), stat = "identity") +
xlab("Year")
tenki2 <- cbind(1969:1990,
t(matrix(tenki$Temperture, nrow = 12)),
t(matrix(tenki$Precipitation, nrow = 12)))
colnames(tenki2) <- c("Year", paste("T", 1:12, sep = ""),
paste("P", 1:12, sep = ""))
## 1970-1990
nenrin2 <- nenrin[(nrow(nenrin) - 20):nrow(nenrin), ]
nenrin3 <- melt(nenrin2, id = "Year", na.rm = TRUE)
p <- ggplot(nenrin3, aes(x = Year, y = value))
p + geom_line(aes(colour = variable)) + ylab("Width (mm)")
## Standardizing
standardize <- function(y) {
x <- seq_along(y)
r <- residuals(lm(y ~ x))
r / sd(r)
}
nenrin.std <- apply(nenrin2[, -1], 2, standardize)
nenrin.m <- apply(nenrin.std, 1, mean)
p <- ggplot(data.frame(Year = 1970:1990, Width = nenrin.m))
p + geom_line(aes(x = Year, y = Width)) + ylab("Width (mm)")
## Current Year (1970-1990)
tenki.c <- tenki2[-1, -1]
## Previous Year (1969-1989)
tenki.p <- tenki2[-21, -1]
colnames(tenki.p) <- c(paste("PT", 1:12, sep = ""),
paste("PP", 1:12, sep = ""))
## Lasso
fit <- glmnet(cbind(tenki.c, tenki.p),
nenrin.m, family = "gaussian",
alpha = 1)
plot(fit, las = 1)
月平均気温（日平均気温の月平均）
月降水量
スギ19本の年輪幅の推移
年輪幅のインデックス
glmnetの結果のプロット
罰則の値をかえて、係数を表示。
> coef(fit, s = 0.01)
49 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 3.0758007606
T1 0.0712534570
T2 .
T3 .
T4 .
T5 .
T6 .
T7 -0.0291844521
T8 .
T9 .
T10 .
T11 .
T12 .
P1 .
P2 .
P3 0.0009881278
P4 .
P5 .
P6 0.0006566845
P7 .
P8 0.0010069106
P9 0.0016475067
P10 -0.0007586308
P11 0.0031013661
P12 .
PT1 .
PT2 .
PT3 .
PT4 -0.3025697658
PT5 0.0534530616
PT6 .
PT7 -0.0394780998
PT8 .
PT9 .
PT10 .
PT11 .
PT12 .
PP1 -0.0039431618
PP2 .
PP3 -0.0007690416
PP4 .
PP5 .
PP6 -0.0003101598
PP7 0.0016621623
PP8 0.0011100792
PP9 0.0014363442
PP10 0.0021611376
PP11 .
PP12 .
> coef(fit, s = 0.05)
49 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 3.264440e+00
T1 7.992246e-02
T2 .
T3 .
T4 .
T5 .
T6 .
T7 .
T8 .
T9 .
T10 .
T11 .
T12 .
P1 .
P2 .
P3 7.148670e-04
P4 .
P5 .
P6 .
P7 .
P8 .
P9 1.271184e-03
P10 -1.747889e-05
P11 2.260789e-03
P12 .
PT1 .
PT2 .
PT3 .
PT4 -2.494856e-01
PT5 .
PT6 .
PT7 -4.685454e-02
PT8 .
PT9 .
PT10 .
PT11 .
PT12 .
PP1 -1.794279e-03
PP2 .
PP3 .
PP4 .
PP5 .
PP6 .
PP7 6.633849e-04
PP8 1.262963e-03
PP9 6.545825e-04
PP10 1.809208e-03
PP11 .
PP12 .
> coef(fit, s = 0.1)
49 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 3.3668797916
T1 0.0332514184
T2 .
T3 .
T4 .
T5 .
T6 .
T7 .
T8 .
T9 .
T10 .
T11 .
T12 .
P1 .
P2 .
P3 .
P4 .
P5 .
P6 .
P7 .
P8 .
P9 0.0010024302
P10 .
P11 0.0016638911
P12 .
PT1 .
PT2 .
PT3 .
PT4 -0.2371797705
PT5 .
PT6 .
PT7 -0.0311172873
PT8 .
PT9 .
PT10 .
PT11 .
PT12 .
PP1 .
PP2 .
PP3 .
PP4 .
PP5 .
PP6 .
PP7 0.0001790516
PP8 0.0009766224
PP9 0.0002329433
PP10 0.0004096242
PP11 .
PP12 .
> coef(fit, s = 0.3)
49 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 1.277959
T1 .
T2 .
T3 .
T4 .
T5 .
T6 .
T7 .
T8 .
T9 .
T10 .
T11 .
T12 .
P1 .
P2 .
P3 .
P4 .
P5 .
P6 .
P7 .
P8 .
P9 .
P10 .
P11 .
P12 .
PT1 .
PT2 .
PT3 .
PT4 -0.091908
PT5 .
PT6 .
PT7 .
PT8 .
PT9 .
PT10 .
PT11 .
PT12 .
PP1 .
PP2 .
PP3 .
PP4 .
PP5 .
PP6 .
PP7 .
PP8 .
PP9 .
PP10 .
PP11 .
PP12 .
前年4月の気温がもっとも影響を与えていそうという結果に。
コメント 0