串本・浦神の潮位差と黒潮流路緯度 Rで線形モデル

前回の続き。今回はrpy2を通してRを使い、串本・浦神の潮位差と黒潮流路緯度を線形モデルとして評価してみる。
(1) データをRに渡し、scatter plotを描く。
(2) 潮位差(rlevel)と黒潮緯度(rlat)の線形モデルを考える。
(3) 線形モデルをもとに、regression line の直線を引く。
(4) 線形モデルの結果を打ち出す。
missing valueもmasked arrayから自動的に選別して除いてくれているようだ。

#read module
#for data
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

#for R
import rpy2.robjects as robjects

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()   


#for R
rlat=robjects.FloatVector(kuroshio_lat)
rlevel=robjects.FloatVector(kushimoto_uragami)
robjects.globalEnv["rlat"] = rlat
robjects.globalEnv["rlevel"] = rlevel
r=robjects.r
r('png(\'kuroshio_reg_R.png\')')    
r('plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\")') #scatter plot 
r('fm<-lm(rlat ~ rlevel)')  # linear model
r('abline(fm)')             # plot regression line
r('dev.off()')            
summary=r('sm<-summary(fm)')
print summary

図

 

summaryの出力

Call:
lm(formula = rlat ~ rlevel)

Residuals:
     Min       1Q   Median       3Q      Max
-3.17592 -0.32259  0.08407  0.40907  1.30406

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 31.955937   0.038044  839.98   <2e-16 ***
rlevel       0.066666   0.003672   18.15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6225 on 532 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared: 0.3825,    Adjusted R-squared: 0.3813
F-statistic: 329.5 on 1 and 532 DF,  p-value: < 2.2e-16

線形モデルの係数は0.066666、切片は31.955937。
ともにp値は小さく(2e-16)、0であるという帰無仮説は棄却できる。
線形モデルのR^2値は約38%で、前回求めた相関係数0.62の2乗になっている。

係数などは以下のようにとりだせる。

print summary.names

 [1] "call"          "terms"         "residuals"     "coefficients"
 [5] "aliased"       "sigma"         "df"            "r.squared"  
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"  "na.action"  

cf=r('cf<-sm$coefficients')
print cf.names

[[1]]
[1] "(Intercept)" "rlevel"   

[[2]]
[1] "Estimate"   "Std. Error" "t value"    "Pr(>|t|)"

for i in range(8):
    print cf[i]

31.9559366577
0.066665747972
0.0380437732903
0.00367237036394
839.978106637
18.1533291486
0.0
1.13483941203e-57

参考：単回帰http://www.geocities.jp/kashiwabara_fact/statistics-tankaiki.htm

関連ポスト
北極振動　Rで統計
北極振動　Rで統計グラフ編