Use R

Size: px

Start display at page:

Download "Use R"

しょうすけさんきち
5 years ago
Views:

1 Use R! 2008/05/23( ) Index Introduction (GLM) ( ) R. Introduction R,, PLS,,, etc. 2. Correlation coefficient (Pearson s product moment correlation) r = Sxy Sxx Syy :, Sxy, Sxx= X, Syy Y % R

2 cor(x, Y, method = c( pearson )) #X, Y # R x <- rnorm(150) y <- rnorm(150) # 150 # x y plot(x, y) # # # cor(x, y, method="pearson") # R 0 p p p p cor.test(x, y, method= pearson ) # p Pearson's product-moment correlation data: x and y t = , df = 148, p-value = alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: sample estimates: cor t-value df Degree of freedom ( 2) p-value( ) t-value t = r ( n 2) (1 r 2 ) t t-value, r n

3 95% 0 0 ( ) R # cor(x, y, method = spearman ) cor(x, y, method = kendall ) # # # p cor.test(x, y, method= spearman ) cor.test(x, y, method= kendall ) method= ( ) X, Y R ( ) ( ) ( ) R airquality airquality cor( airquality, method= pearson ) cor( airquality, method= pearson, use= pairwise ) NA

4 3. hist ( ) airquality airquality attach(airquality) names(airquality) R y β + α + ε ε y~x ~ 0 0 Y~X+0 # (Wind Ozone ) plot(wind, Ozone) # plot(x, Y ) result <- lm(ozone~wind) # result abline(result) # result # Call: lm(formula = Ozone ~ Wind) Coefficients: (Intercept) Wind # Ozone=-5.551Wind Summary summary(result) Call: lm(formula = Ozone ~ Wind) Residuals: # ( ( ) 1(3) ) Min 1Q Median 3Q Max Coefficients: # Estimate Std. Error t value Pr(> t ) (Intercept) < 2e-16 *** # 0 Wind e-13 *** # 0 S ignif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: on 114 degrees of freedom Multiple R-Squared: , Adjusted R-squared: #R 2 ajsr 2 F-statistic: on 1 and 114 DF, p-value: 9.272e-13 # p

5 new <- data.frame(wind = seq(min(wind), max(wind), by = 0.1)) # cline <-predict(result, new, interval="confidence") # cline plot(wind, Ozone, xlab = "Wind", ylab = "Ozone") # matplot(new, cline, lty=c(1,2,2), type="l", add=t) # ( ) (X) ( ) (x) Ozone Wind Ozone (NO) y x Wind Ozone Ozone Wind 2 (Major axis regression) 4. R R 2 ( ) attach(airquality)

6 lm (Linear Model ) summary summary(result) Call: lm(formula = Ozone ~ Wind + Temp) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) # (Intercept) ** Wind e-05 *** Temp e-11 *** --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: on 113 degrees of freedom Multiple R-Squared: , Adjusted R-squared: F-statistic: 74.5 on 2 and 113 DF, p-value: < 2.2e-16 Ozone = 1.84 * Temp 3.06*Wind 71.0 R 2 =0.57 Multi-co linearity - - ( 0) Wind Temp 1 ( 0 ) 0

7 VIF ( Variance Inflation factor ) VIF 10 1 VIF = 1 2 R j Rj 2 ; Xj ( N-1 ) R DAAG VIF R DAAG # attach(airquality) result <- lm(ozone~wind+temp+solar.r) vif(result) # 3 # VIF Wind Temp Solar.R detach(airquality) #detach 10 VIF Ozone = α*wind + β*temp + γ + e Wind Temp Wind Temp Ozone = Wind + Temp + Wind Temp = (1 + Temp)Wind + Temp Wind Temp Wind Wind Temp Wind Temp

8 R attach(airquality) Wind2 <- (Wind ave(wind)) # ( Temp2 <- (Temp ave(temp)) result <- lm(ozone~wind+temp+wind2*temp2) summary(result) # Call: lm(formula = Ozone ~ Wind + Temp + Wind2 * Temp2) Coefficients: (2 not defined because of singularities) Estimate Std. Error t value Pr(> t ) (Intercept) *** Wind e-06 *** Temp e-12 *** Wind2 NA NA NA NA Temp2 NA NA NA NA Wind2:Temp e-05 *** # --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: on 112 degrees of freedom Multiple R-Squared: , Adjusted R-squared: F-statistic: on 3 and 112 DF, p-value: < 2.2e-16 detach(airquality) 2 (R 2 ) 1.00 (Over fitting ) F AIC AIC AIC (GLM)

9 AIC(Akaike s Information Criteria) AIC = 2*p 2 ( 1/2 *n*log(2 2 )) p () 2 = n 2 ( e ) n i= 1 i / / AIC ( ) ( ) AIC p 2 P AIC wle AIC mle.aic HP AIC # Wind Temp 2 attach(airquality) Wind2 <- (Wind ave(wind)) # Temp2 <- (Temp ave(temp) ) WinTem <- (Wind2 * Temp2) # result <- lm(ozone~wind+temp+wintem) # select <- mle.aic(result) # select summary(select) # select Call: mle.aic(formula = Ozone ~ Wind + Temp + WinTem) Akaike Information Criterion (AIC): (Intercept) Wind Temp WinTem aic # (1) (0) [1,] [2,] [3,] [4,] [5,] [6,] [7,] [8,] [9,] [10,] # Printed the first 15 best models AIC

10 5. (GLM, Generalized Linear Model) GLM lm glm lm GLM (Binomial) (Poisson) (Gamma) ( ) ( ) HP lm 6 nls Non-linear Least Square ( ) ( ) ( ) (GAM) ( ) nls nls lm 2 R nls SSlogis ( ) ( ( )) a f ( x i ) = + ε cx i i 1 + b * e

11 a, b, c: ( ) c >0, ε: 0 ( ) ( ) ( ) dy = dx α β ( y ) y β a= α/β, α/β β α y( ) α/β 0 y # X <- (1:20) Y <- c(3, 16, 54, 139, 263, 420, 611, 750, 860, 900, 940, 950, 960, 930, 980, 970, 989, 980, 975, 980) Y <- Y + rnorm(20, 40, 20) # plot(y~x) # plot(x, Y) # result <- nls (Y ~ a / (1+b*exp(-c*X)), start = c(a=1000, b=100, c=1), trace=t) summary(result) Formula: Y ~ a/(1 + b * exp(-c * X)) Parameters: Estimate Std. Error t value Pr(> t ) a < 2e-16 *** b e-07 *** c < 2e-16 *** # --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: on 17 degrees of freedom # AIC AIC(result) lines(x, fitted(result)) # (1000 a=1000 )

12 @ glm 0 X( ) X ( ) ( ) glm GLM GLMM GLM 7. R R help! help example( ) example( ) help HP PC 8. The R-tips,2005, 9-ten (Web ) Statistics an introduction using R M.J. Crawley, 2005, John Wiley & Sons,, HP # HP # HP HP

13 x y y= x= ( ) Smatr plot(airquality$wind, airquality$ozone, xlab = "Wind", ylab = "Ozone") result <- lm(airquality$ozone~airquality$wind) # abline(result) result2 <- lm(airquality$wind~airquality$ozone) # abline(-result2$coef[1]/result2$coef[2], 1/result2$coef[2], lty=2) # #SMA SMA <- f unction(x, y) { # 2 dt <- data.frame(x, y) dt2 <- na.omit(dt) x1 <- dt2$x y1 <- dt2$y slope <- sign(cor(x1, y1))*sqrt(var(y1)/v ar(x1)) # intercept <- mean(y1)-slope*mean(x1) # return(list(slope=slope, Intercept=intercept)) } result3 <- SMA(airquality$Wind, airquality$ozone) result3 abline(result3$intercept, result3$slope, col="red") # (SMA) # # # smart library(smatr) # res <- line.cis(airquality$ozone, airquality$wind) # (SMA) res # coef(sma) lower limit upper limit # 95 ( ) elevation slope # slope.test(airquality$ozone, airquality$wind, test.value = -5.55) $r [1] $p [1] e-10 # $test.value [1] #

(lm) lm AIC 2 / 1

W707 s-taiji@is.titech.ac.jp 1 / 1 (lm) lm AIC 2 / 1 : y = β 1 x 1 + β 2 x 2 + + β d x d + β d+1 + ϵ (ϵ N(0, σ 2 )) y R: x R d : β i (i = 1,..., d):, β d+1 : ( ) (d = 1) y = β 1 x 1 + β 2 + ϵ (d > 1) y