http://d.hatena.ne.jp/matsuou1/20110418/1303144024
https://to-kei.net/r-beginner/r-3/#lm-2
2018-07-26_09-34
/home/manda
manda@calypso
$ R
> data=read.table("SCAT_SST.RAIN.ACCUM.AAVE.2018-07-23_120-129_24_34.txt", header=FALSE, sep="\t")
> data
V1 V2 V3 V4
1 UKMO 571 32 27.085
2 DR_B 565 31 26.704
3 OI 549 29 26.760
4 MGD 532 32 27.281
5 RTG5 527 34 26.898
6 HIM 522 36 27.454
7 JPL 513 40 27.461
8 NAVO 513 41 27.418
9 ABOM 462 39 27.314
10 RTG8 421 35 26.990
> mean(data[,4])
[1] 27.1365
> mean(data[,3])
[1] 34.9
> sd(data[,3])
[1] 4.067486
> sd(data[,4])
[1] 0.2882014
[1] 0.2882014
> y<-(data[,3])
> x<-(data[,4])
> y
[1] 32 31 29 32 34 36 40 41 39 35
> x
[1] 27.085 26.704 26.760 27.281 26.898 27.454 27.461 27.418 27.314 26.990
> ans <- lm(y~x)
> ans
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-265.94 11.09
> s.ans<-summary(ans)
> s.ans
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-4.502 -2.178 1.199 1.739 2.979
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -265.945 83.798 -3.174 0.01312 *
x 11.086 3.088 3.590 0.00708 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.67 on 8 degrees of freedom
Multiple R-squared: 0.617, Adjusted R-squared: 0.5692
F-statistic: 12.89 on 1 and 8 DF, p-value: 0.007081
Estimate その回帰分析による、回帰係数の推定値（推定値の実現値）。
Multiple R-squared 決定係数R^2であり、単回帰分析の場合は、xとyの相関係数の2乗がこれを表している。
Pr(>|t|) p値。回帰係数が0（説明変数は目的変数に影響を与えない）であるという帰無仮説に対する仮設検定結果のp値。
Residual standard error 誤差項の標準偏差の推定値。
std.Error 各回帰係数の推定量の標準誤差。
t value t値。回帰係数=0（説明変数は目的変数に影響を与えない）であるという帰無仮説に対するt検定によって計算される値。
> coe <- s.ans$coefficient
> aic <- AIC(ans)
> N <- length(y)
> result <- cbind(coe,aic,N)
> result[2,5:6] <- ""
coe <- s.ans$coefficient#回帰係数を抽出
aic <- AIC(ans)#AICを抽出
N <- length(y)#解析したデータの総数を抽出
result <- cbind(coe,aic,N) #結果をまとめる
result[2,5:6] <- "" #2行目の５列目と６列目は、上の行と同じなので空にする。
ここでは回帰係数とAICとデータの総数を抽出して変数に格納したあと、resultという変数に全てまとめています。そこで”cbind()”という関数を使っています。この関数は、行数が同じ行列同士を横にしてくっつけて新しい行列を作るというものです。ここではcoe,aic,Nをくっつけています。
> write.table(matrix(c("",colnames(result)),nrow=1),"SCAT_SST.RAIN.ACCUM.AAVE.2018-07-23_120-129_24_34.csv",append=T,quote=F,sep=",",row.names=F,col.names=F)
> write.table(result,"SCAT_SST.RAIN.ACCUM.AAVE.2018-07-23_120-129_24_34.csv",append=T,quote=F,sep=",",row.names=T,col.names=F)
> q()
2018-07-26_10-28
/work05/manda/WRF.POST/GrADS/Kyushu.170704.Bosai/BAR.CHART_RAIN.ACCUM.AAVE
manda@calypso
$ cat SCAT_SST.RAIN.ACCUM.AAVE.2018-07-23_120-129_24_34.csv
,Estimate,Std. Error,t value,Pr(>|t|),aic,N
(Intercept),-265.944895694614,83.7980818970961,-3.17363941600948,0.0131225298375126,51.7872607744775,10
x,11.0863558563048,3.08786427397356,3.5902989486123,0.00708141044422536,,
GMTで回帰直線を描く場合のサンプル
(このままでは動かない)
echo
echo PREPROCESS FOR PLOTTING LINEAR REGRESSION LINE
echo
in2=$(basename $in .txt).csv
if [ ! -f $in2 ]; then
echo Error in $0 : No such file, $in2
exit 1
fi
echo
echo READ LINEAR REGRESSION
ls -lh $in2
echo
b=$(awk -F, '{if(NR==2) print $2}' $in2)
a=$(awk -F, '{if(NR==3) print $2}' $in2)
echo "a=$a"
echo "b=$b"
echo
echo SET X AND Y OF LINEAR REGRESSION
echo
xmin=$(cat $in|awk 'BEGIN{m=100000}{if($1!="#" && m>$4) m=$4} END{print m}')
xmax=$(cat $in|awk '{if($1!="#" && m<$4) m=$4} END{print m}')
y1=$(echo "$a*$xmin + $b;scale=3"|bc)
y2=$(echo "$a*$xmax + $b;scale=3"|bc)
dx=$(echo "$xmax - $xmin ;scale=3"|bc)
dy=$(echo "$y2 - $y1 ;scale=3"|bc)
echo
echo "xmin = $xmin y1 = $y1"
echo "xmax = $xmax y2 = $y2"
echo "dx = $dx dy = $dy"
echo
echo DONE.
echo PLOT LINEAR REGRESSION LINE
psxy <<EOF -R -JX -W1/${CLR} -K -O >>$out
$xmin $y1
$xmax $y2
EOF
psbasemap -R$range -JX$size -B$anot -O -K >>$out