R: Kaplan-Meier曲線を作成する

今回はRを使って、論文用のきれいなKaplan-Meier曲線を作成するコマンドを紹介します。最終的にはこんなグラフが出来上がります。

quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
axis(1, at=XLIM*365.25, labels=XLABEL)
axis(2, las=1)
title(main="  a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan15.png", type="png")

StataでのKaplan-Meir曲線の作成方法は、「Stata: Kaplan-Meier曲線を作成する」をご覧下さい。

サンプルデータの読み込み

今回はsurvivalパッケージのサンプルデータflchainを利用します。data()コマンドを使って、flchainデータを読みこみます。さらにヘルプコマンドを使って、flchainデータの詳細を確認します。

> library(survival)
> data("flchain", package="survival"
> ?flchain

flchain                package:survival                R Documentation

Assay of serum free light chain for 7874 subjects.

Description:

     This is a stratified random sample containing 1/2 of the subjects
     from a study of the relationship between serum free light chain
     (FLC) and mortality.  The original sample contains samples on
     approximately 2/3 of the residents of Olmsted County aged 50 or
     greater.

Usage:

     data(flchain)
     
Format:

     A data frame with 7874 persons containing the following variables.

     'age ' age in years

     'sex' F=female, M=male

     'sample.yr' the calendar year in which a blood sample was obtained
Usage:

     data(flchain)
     
Format:

     A data frame with 7874 persons containing the following variables.

     'age ' age in years

     'sex' F=female, M=male

     'sample.yr' the calendar year in which a blood sample was obtained

     'kappa' serum free light chain, kappa portion

     'lambda' serum free light chain, lambda portion

     'flc.grp' the FLC group for the subject, as used in the original
          analysis

     'creatinine' serum creatinine

     'mgus' 1 if the subject had been diagnosed with monoclonal
          gammapothy (MGUS)

     'futime' days from enrollment until death.  Note that there are 3
          subjects whose sample was obtained on their death date.

     'death' 0=alive at last contact date, 1=dead

     'chapter' for those who died, a grouping of their primary cause of
          death by chapter headings of the International Code of
          Diseases ICD-9

Details:


     'sample.yr' the calendar year in which a blood sample was obtained

     'kappa' serum free light chain, kappa portion

     'lambda' serum free light chain, lambda portion

     'flc.grp' the FLC group for the subject, as used in the original
          analysis

     'creatinine' serum creatinine

     'mgus' 1 if the subject had been diagnosed with monoclonal
          gammapothy (MGUS)

     'futime' days from enrollment until death.  Note that there are 3
     only found by chance from a test (serum electrophoresis) which is
     ordered for other causes.  Later work suggested that one component
     of immunoglobulin production, the serum free light chain, might be
     a possible marker for immune disregulation.  In 2010 Dr. Angela
     Dispenzieri and colleagues assayed FLC levels on those samples
     from the original study for which they had patient permission and
     from which sufficient material remained for further testing.  They
     found that elevated FLC levels were indeed associated with higher
     death rates.

     Patients were recruited when they came to the clinic for other
     appointments, with a final random sample of those who had not yet
     had a visit since the study began.  An interesting side question
     is whether there are differences between early, mid, and late
     recruits.

     This data set contains an age and sex stratified random sample
     that includes 7874 of the original 15759 subjects.  The original
     subject identifiers and dates have been removed to protect patient
     identity.  Subsampling was done to further protect this
     Proceedings 87:512-523.

     R Kyle, T Therneau, SV Rajkumar, D Larson, M Plevak, J Offord, A
     Dispenzieri, J Katzmann, and LJ Melton, III, 2006, Prevalence of
     monoclonal gammopathy of undetermined significance, New England J
     Medicine 354:1362-1369.

Examples:

     data(flchain)
     age.grp <-  cut(flchain$age, c(49,54, 59,64, 69,74,79, 89, 110),
                    labels= paste(c(50,55,60,65,70,75,80,90),
                                  c(54,59,64,69,74,79,89,109), sep='-'))
     table(flchain$sex, age.grp)

qキーをクリックすれば、いつもの入力画面に復帰できます。とりあえず、head()、summary()、str()を使って、flchainデータセットの中身を確認します。

> head(flchain)
  age sex sample.yr kappa lambda flc.grp creatinine mgus futime death
1  97   F      1997  5.70  4.860      10        1.7    0     85     1
2  92   F      2000  0.87  0.683       1        0.9    0   1281     1
3  94   F      1997  4.36  3.850      10        1.4    0     69     1
4  92   F      1996  2.42  2.220       9        1.0    0    115     1
5  93   F      1996  1.32  1.690       6        1.1    0   1039     1
6  90   F      1997  2.01  1.860       9        1.0    0   1355     1
      chapter
1 Circulatory
2   Neoplasms
3 Circulatory
4 Circulatory
5 Circulatory
6      Mental

> summary(flchain)
      age         sex        sample.yr        kappa            lambda      
 Min.   : 50.00   F:4350   Min.   :1995   Min.   : 0.010   Min.   : 0.040  
 1st Qu.: 55.00   M:3524   1st Qu.:1996   1st Qu.: 0.960   1st Qu.: 1.200  
 Median : 63.00            Median :1996   Median : 1.270   Median : 1.510  
 Mean   : 64.29            Mean   :1997   Mean   : 1.431   Mean   : 1.703  
 3rd Qu.: 72.00            3rd Qu.:1997   3rd Qu.: 1.680   3rd Qu.: 1.920  
 Max.   :101.00            Max.   :2003   Max.   :20.500   Max.   :26.600  
                                                                           
    flc.grp         creatinine          mgus            futime    
 Min.   : 1.000   Min.   : 0.400   Min.   :0.0000   Min.   :   0  
 1st Qu.: 3.000   1st Qu.: 0.900   1st Qu.:0.0000   1st Qu.:2852  
 Median : 5.000   Median : 1.000   Median :0.0000   Median :4302  
 Mean   : 5.471   Mean   : 1.093   Mean   :0.0146   Mean   :3661  
 3rd Qu.: 8.000   3rd Qu.: 1.200   3rd Qu.:0.0000   3rd Qu.:4773  
 Max.   :10.000   Max.   :10.800   Max.   :1.0000   Max.   :5215  
                  NA's   :1350                                    
     death               chapter    
 Min.   :0.0000   Circulatory: 745  
 1st Qu.:0.0000   Neoplasms  : 567  
 Median :0.0000   Respiratory: 245  
 Mean   :0.2755   Mental     : 144  
 3rd Qu.:1.0000   Nervous    : 130  
 Max.   :1.0000   (Other)    : 338  
                  NA's       :5705  

> str(flchain)
'data.frame':	7874 obs. of  11 variables:
 $ age       : num  97 92 94 92 93 90 90 90 93 91 ...
 $ sex       : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
 $ sample.yr : num  1997 2000 1997 1996 1996 ...
 $ kappa     : num  5.7 0.87 4.36 2.42 1.32 2.01 0.43 2.47 1.91 0.791 ...
 $ lambda    : num  4.86 0.683 3.85 2.22 1.69 1.86 0.88 2.7 2.18 2.22 ...
 $ flc.grp   : num  10 1 10 9 6 9 1 10 9 6 ...
 $ creatinine: num  1.7 0.9 1.4 1 1.1 1 0.8 1.2 1.2 0.8 ...
 $ mgus      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ futime    : int  85 1281 69 115 1039 1355 2851 372 3309 1326 ...
 $ death     : num  1 1 1 1 1 1 1 1 1 1 ...
 $ chapter   : Factor w/ 16 levels "Blood","Circulatory",..: 2 13 2 2 2 11 11 14 15 2 ...

生存時間とアウトカムの設定

まずは、追跡期間とアウトカムをSurv()で設定します。アウトカム変数は、0=打ち切り、1=アウトカム発症で定義しておくのが無難です。

> s1 <- Surv(flchain$futime, flchain$death)

暴露因子としてflc.grpを設定します。

> m1 <- survfit(s1 ~ flc.grp, data=flchain)

m1オブジェクトにはいろいろな情報が詰まっています。str()で確認できます。

List of 18
 $ n         : int [1:10] 769 811 820 786 791 791 806 730 803 767
 $ time      : num [1:6191] 1 3 6 9 15 27 30 37 63 88 ...
 $ n.risk    : num [1:6191] 769 768 767 766 765 763 762 761 760 759 ...
 $ n.event   : num [1:6191] 0 0 0 1 2 1 0 1 0 1 ...
 $ n.censor  : num [1:6191] 1 1 1 0 0 0 1 0 1 0 ...
 $ surv      : num [1:6191] 1 1 1 0.999 0.996 ...
 $ std.err   : num [1:6191] 0 0 0 0.00131 0.00227 ...
 $ cumhaz    : num [1:6191] 0 0 0 0.00131 0.00392 ...
 $ std.chaz  : num [1:6191] 0 0 0 0.00131 0.00226 ...
 $ start.time: num 0
 $ strata    : Named int [1:10] 530 585 609 611 596 615 651 631 672 691
  ..- attr(*, "names")= chr [1:10] "flc.grp=1" "flc.grp=2" "flc.grp=3" "flc.grp=4" ...
 $ type      : chr "right"
 $ logse     : logi TRUE
 $ conf.int  : num 0.95
 $ conf.type : chr "log"
 $ lower     : num [1:6191] 1 1 1 0.996 0.992 ...
 $ upper     : num [1:6191] 1 1 1 1 1 ...
 $ call      : language survfit(formula = s1 ~ flc.grp, data = flchain)
 - attr(*, "class")= chr "survfit"

累積アウトカム発症率の算出

暴露因子の各群における累積生存率を算出しましょう。研究開始から1年ごとに3年時までの累積生存率を算出します。

> summary(m1, time=seq(0, 3, 1)*365.25)
Call: survfit(formula = s1 ~ flc.grp, data = flchain)

                flc.grp=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    769       0    1.000 0.00000        1.000        1.000
  365    752      11    0.986 0.00431        0.977        0.994
  730    744       3    0.982 0.00485        0.972        0.991
 1096    730       9    0.970 0.00621        0.958        0.982

                flc.grp=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    811       0    1.000 0.00000        1.000        1.000
  365    797       9    0.989 0.00369        0.982        0.996
  730    789       7    0.980 0.00491        0.971        0.990
 1096    778       7    0.971 0.00587        0.960        0.983

                flc.grp=3 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    820       0    1.000 0.00000        1.000        1.000
  365    805       9    0.989 0.00365        0.982        0.996
  730    791      11    0.975 0.00543        0.965        0.986
 1096    778      11    0.962 0.00672        0.949        0.975

                flc.grp=4 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    786       0    1.000 0.00000        1.000        1.000
  365    768      14    0.982 0.00473        0.973        0.991
  730    756       7    0.973 0.00578        0.962        0.985
 1096    745       8    0.963 0.00677        0.950        0.976

                flc.grp=5 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    791       0    1.000 0.00000        1.000        1.000
  365    768      19    0.976 0.00546        0.965        0.987
  730    758       5    0.970 0.00613        0.958        0.982
 1096    741      13    0.953 0.00757        0.938        0.968


                flc.grp=6 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    791       0    1.000 0.00000        1.000        1.000
  365    774      12    0.985 0.00437        0.976        0.993
  730    752      14    0.967 0.00640        0.954        0.979
 1096    735      14    0.949 0.00789        0.933        0.964

                flc.grp=7 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    806       0    1.000 0.00000        1.000        1.000
  365    778      18    0.977 0.00524        0.967        0.988
  730    758      17    0.956 0.00726        0.942        0.970
 1096    740      16    0.936 0.00869        0.919        0.953

                flc.grp=8 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    730       0    1.000 0.00000        1.000        1.000
  365    696      27    0.963 0.00703        0.949        0.977
  730    675      15    0.942 0.00870        0.925        0.959
 1096    658      13    0.924 0.00988        0.905        0.943

                flc.grp=9 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    803       0    1.000  0.0000        1.000        1.000
  365    760      40    0.950  0.0077        0.935        0.965
  730    728      31    0.911  0.0101        0.892        0.931
 1096    700      27    0.877  0.0116        0.855        0.900

                flc.grp=10 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    767       3    0.996 0.00225        0.992        1.000
  365    651     105    0.858 0.01262        0.834        0.884
  730    583      62    0.776 0.01512        0.747        0.807
 1096    532      49    0.711 0.01648        0.679        0.744

各暴露群の累積発症率のデータフレームを作成してみましょう。names()コマンドでオブジェクトcp1 (cumulative probability)内のどこに累積生存率、95%信頼区間等が保存されているかを確認します。

> cp1 <- summary(m1, time=seq(0, 3, 1)*365.25)
> names(cp1)
 [1] "n"             "time"          "n.risk"        "n.event"      
 [5] "n.censor"      "surv"          "std.err"       "cumhaz"       
 [9] "std.chaz"      "start.time"    "strata"        "type"         
[13] "logse"         "conf.int"      "conf.type"     "lower"        
[17] "upper"         "call"          "table"         "rmean.endtime"

いろいろ確認すると、必要な情報は下記に保管されています。

  • 累積生存率:surv
  • 95%信頼区間:lower, upper
  • 暴露群のカテゴリ名:strata
  • 時間:time

上記の変数をまとめて、データフレームdfを作成します。

> df <- data.frame(group=cp1$strata, day=cp1$time, cp=cp1$surv, cpl=cp1$lower, cpu=cp1$upper)
> head(df)
      group     day        cp       cpl       cpu
1 flc.grp=1    0.00 1.0000000 1.0000000 1.0000000
2 flc.grp=1  365.25 0.9856156 0.9772125 0.9940910
3 flc.grp=1  730.50 0.9816731 0.9722073 0.9912312
4 flc.grp=1 1095.75 0.9697498 0.9576508 0.9820017
5 flc.grp=2    0.00 1.0000000 1.0000000 1.0000000
6 flc.grp=2  365.25 0.9888582 0.9816461 0.9961233

累積生存率を累積アウトカム発症率に変換したければ、1からcp、cpl、cpuを引き算すればOKです。ただし、cplとcpuが逆になってしまうので、ご注意ください。

> df$cp <- 1 - df$cp
> df$cpl <- 1 - df$cpl
> df$cpu <- 1 - df$cpu
> head(df)
      group     day         cp        cpl         cpu
1 flc.grp=1    0.00 0.00000000 0.00000000 0.000000000
2 flc.grp=1  365.25 0.01438439 0.02278754 0.005908989
3 flc.grp=1  730.50 0.01832687 0.02779273 0.008768846
4 flc.grp=1 1095.75 0.03025016 0.04234917 0.017998293
5 flc.grp=2    0.00 0.00000000 0.00000000 0.000000000
6 flc.grp=2  365.25 0.01114179 0.01835391 0.003876691

最後にdayをyearに変換します。

> df$yr <- df$day/365.25
> head(df)
      group     day         cp        cpl         cpu yr
1 flc.grp=1    0.00 0.00000000 0.00000000 0.000000000  0
2 flc.grp=1  365.25 0.01438439 0.02278754 0.005908989  1
3 flc.grp=1  730.50 0.01832687 0.02779273 0.008768846  2
4 flc.grp=1 1095.75 0.03025016 0.04234917 0.017998293  3
5 flc.grp=2    0.00 0.00000000 0.00000000 0.000000000  0
6 flc.grp=2  365.25 0.01114179 0.01835391 0.003876691  1

Kaplan-Meier曲線の作成

とりあえず、デフォルト設定のままでKaplan-Meier曲線を作成します。3.7インチの正方形のデバイスの上にKM曲線を作成し、png形式で出力します。

> quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
> plot(m1)
> quartz.save(file="kaplan1.png", type="png")

累積生存率ではなく、累積アウトカム発症率のKM曲線を作成します。

> plot(m1, fun="event")

累積アウトカム発症率の最大値に計算して、小数点1桁レベルで切り上げてしまい、Y軸の最大値を調整します。累積アウトカム発症率の最大値が0.69なので、Y軸の最大値 YLIMを0.70に設定します。

> YLIM <- ceiling( (1 - min(m1$surv))*10 )/10
> plot(m1, fun="event", ylim=c(0, YLIM))

曲線がすべて黒実践なので、区別がつきません。Rでは下記の6種類のタイプを選択できます。

  • lty=1: solid
  • lty=2: dashed
  • lty=3: dotted
  • lty=4: dotdash
  • lty=5: longdash
  • lty=6: twodash

flc.grpのカテゴリ数は10なので、線の種類が足りません。そこで、2種類の色 (black = #000000とgray = #999999) を使うことによって、10種類の曲線を区別することにします。

%%の戻り値は割り算の余りであり (5%%3 = 2)、%/%の戻り値は商です (5%/%3 = 1) 。

> exname <- names(m1$strata)
> ex     <- 1:length(exname)
> N <- 6
> exlty  <- ex %% N
> exlty <- ifelse(exlty==0, N, exlty)
>
> excol <- (ex - 1) %/% N
> excol <- ifelse(excol==0, "#000000", excol)
> excol <- ifelse(excol==1, "#999999", excol)
>
> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)

凡例の追加

凡例を追加します。

> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)
> legend("topleft", exname, lty=exlty, col=excol)

凡例の枠線を消します。また、凡例をrev()コマンドを使って降順ソートします。

> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol)
> legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")

軸の調整

軸の表示方法を調整します。まずは軸無しのグラフを作成します。

> plot(km1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
> legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")

X軸の単位を年にします。追跡期間はm1$timeに保存されており、最大値は5215日=  14.3年です。X軸のラベルは0〜14にします。

> XLIM <- 0:floor(max(m1$time) / 365.25)
> axis(1, at=XLIM*365.25, labels=XLIM)

Y軸を追加します。軸の数値を水平に表示します。

> axis(2, las=1)

デフォルト設定のまま、メインタイトルととxy軸のタイトルを追加するとこうなります。

title(main="a", xlab="Observational period (year)", ylab="Cumulative probability of the incidence of outcome")

メインタイトルをグラフの左上に記載するためには、figure regionの外側であるouter region (outer=TRUE) に、左揃え(adj=0) で、1行だけ内側に記載する設定が必要である。

title(main="a", line=-1, adj=0, outer=TRUE)

X軸のタイトルは、X軸から少し離れすぎですので、もう少し軸に近づけましょう。X軸から2行分だけ外向けに配置します。

title(xlab="Observational period (year)", line=2)

Y軸のラベルは長すぎるので途中で改行します。改行した場所に”\n”を挿入します。

> title(ylab="Cumulative probability of\nthe incidence of outcome", line=2.5)

余白の設定

figure regionの上下右の余白が大きいので、調整します。上下左右の余白のデフォルト値はpar()$marで確認できます。

> quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
> par()$mar
[1] 5.1 4.1 4.1 2.1

下5.1、左4.1、上4.1、右2.1に設定されています。

いろいろ微調整してみましたが、mar=c(3, 3.5, 1.5, 0.2)くらいが良さそうです。

quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
axis(1, at=XLIM*365.25, labels=XLIM)
axis(2, las=1)
title(main="  a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan14.png", type="png")
graphics.off()

左右の余白が小さくなった分、figure regionが大きくなり、X軸のメモリの表示がおかしなことになっていますので、微調整します。偶数年だけ表示するようにX軸のラベルを設定します。2で割って余りが0ならば表示するように設定すればいいのです。

XLABEL <- ifelse(XLIM%%2==0, XLIM, "")

このXLABELをX軸のラベルに設定します。

quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))
plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
XLIM <- 0:floor(max(m1$time) / 365.25)
XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
axis(1, at=XLIM*365.25, labels=XLABEL)
axis(2, las=1)
title(main="  a", line=-1, adj=0, outer=TRUE)
title(xlab="Observational period (year)", line=2)
title(ylab="Cumulative probability of the incidence of outcome", line=2.5)
quartz.save(file="kaplan15.png", type="png")

XY軸の目盛り線がちょっと長いのが気になりますが、今回はこれで完成とします。

今回のコマンドをplotKaplanMeier()にまとめてみました。

s1 <- Surv(flchain$futime, flchain$death)
m1 <- survfit(s1 ~ flc.grp, data=flchain)

quartz("km", height=3.7, width=3.7, pointsize=9, bg="white")
par(mar=c(3, 3.5, 1.5, 0.2))

MAIN   <- "  a"
XTITLE <- "Observational period (year)"
YTITLE <- "Cumulative probability of the incidence of outcome"

plotKaplanMeier <- function(m1, MAIN, XTITLE, YTITLE){
  YLIM   <- ceiling( (1 - min(m1$surv))*10 )/10
  exname <- names(m1$strata)
  ex     <- 1:length(exname)
  N      <- 6
  exlty  <- ex %% N
  exlty  <- ifelse(exlty==0, N, exlty)
  excol  <- (ex - 1) %/% N
  excol  <- ifelse(excol==0, "#000000", excol)
  excol  <- ifelse(excol==1, "#999999", excol)

  plot(m1, fun="event", ylim=c(0, YLIM), lty=exlty, col=excol, axes=FALSE)
  
  legend("topleft", rev(exname), lty=rev(exlty), col=rev(excol), bty="n")
  XLIM   <- 0:floor(max(m1$time) / 365.25)
  XLABEL <- ifelse(XLIM%%2==0, XLIM, "")
  axis(1, at=XLIM*365.25, labels=XLABEL)
  axis(2, las=1)
  title(main=MAIN, line=-1, adj=0, outer=TRUE)
  title(xlab=XTITLE, line=2)
  title(ylab=YTITLE, line=2.5)
}

plotKaplanMeir(m1, MAIN, XTITLE, YTITLE)

quartz.save(file=FILENAME, type="png")

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください