Comments on R script

不要顯示程式碼、程式碼要執行

echo=FALSE
eval=TRUE

環境設定：超過一行的字串用" "、顯示三個位數、控制顯示寬度

options(continue="  ")
options(digits=3)
options(width=72) # narrow output

讀取檔案、使用dplyr、選取變項

ds = read.csv("http://www.amherst.edu/~nhorton/r2/datasets/help.csv")
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

newds = select(ds, cesd, female, i1, i2, id, treat, f1a, f1b, f1c, f1d, f1e, 
  f1f, f1g, f1h, f1i, f1j, f1k, f1l, f1m, f1n, f1o, f1p, f1q, f1r, f1s, f1t)

顯示檔案變項名稱、顯示檔案前十個變項的資料

names(newds)

##  [1] "cesd"   "female" "i1"     "i2"     "id"     "treat"  "f1a"   
##  [8] "f1b"    "f1c"    "f1d"    "f1e"    "f1f"    "f1g"    "f1h"   
## [15] "f1i"    "f1j"    "f1k"    "f1l"    "f1m"    "f1n"    "f1o"   
## [22] "f1p"    "f1q"    "f1r"    "f1s"    "f1t"

str(newds[,1:10]) # structure of the first 10 variables

## 'data.frame':    453 obs. of  10 variables:
##  $ cesd  : int  49 30 39 15 39 6 52 32 50 46 ...
##  $ female: int  0 0 0 1 0 1 1 0 1 0 ...
##  $ i1    : int  13 56 0 5 10 4 13 12 71 20 ...
##  $ i2    : int  26 62 0 5 13 4 20 24 129 27 ...
##  $ id    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ treat : int  1 1 0 0 0 1 0 1 0 1 ...
##  $ f1a   : int  3 3 3 0 3 1 3 1 3 2 ...
##  $ f1b   : int  2 2 2 0 0 0 1 1 2 3 ...
##  $ f1c   : int  3 0 3 1 3 1 3 2 3 3 ...
##  $ f1d   : int  0 3 0 3 3 3 1 3 1 0 ...

描述前十個變項：最大值、最小值、平均數、中位數、第一、三四分位數

summary(newds[,1:10]) # summary of the first 10 variables

##       cesd          female            i1              i2       
##  Min.   : 1.0   Min.   :0.000   Min.   :  0.0   Min.   :  0.0  
##  1st Qu.:25.0   1st Qu.:0.000   1st Qu.:  3.0   1st Qu.:  3.0  
##  Median :34.0   Median :0.000   Median : 13.0   Median : 15.0  
##  Mean   :32.8   Mean   :0.236   Mean   : 17.9   Mean   : 22.6  
##  3rd Qu.:41.0   3rd Qu.:0.000   3rd Qu.: 26.0   3rd Qu.: 32.0  
##  Max.   :60.0   Max.   :1.000   Max.   :142.0   Max.   :184.0  
##        id          treat            f1a            f1b      
##  Min.   :  1   Min.   :0.000   Min.   :0.00   Min.   :0.00  
##  1st Qu.:119   1st Qu.:0.000   1st Qu.:1.00   1st Qu.:0.00  
##  Median :233   Median :0.000   Median :2.00   Median :1.00  
##  Mean   :233   Mean   :0.497   Mean   :1.63   Mean   :1.39  
##  3rd Qu.:348   3rd Qu.:1.000   3rd Qu.:3.00   3rd Qu.:2.00  
##  Max.   :470   Max.   :1.000   Max.   :3.00   Max.   :3.00  
##       f1c            f1d      
##  Min.   :0.00   Min.   :0.00  
##  1st Qu.:1.00   1st Qu.:0.00  
##  Median :2.00   Median :1.00  
##  Mean   :1.92   Mean   :1.56  
##  3rd Qu.:3.00   3rd Qu.:3.00  
##  Max.   :3.00   Max.   :3.00

從頭顯示前三筆資料

head(newds, n=3)

##   cesd female i1 i2 id treat f1a f1b f1c f1d f1e f1f f1g f1h f1i f1j
## 1   49      0 13 26  1     1   3   2   3   0   2   3   3   0   2   3
## 2   30      0 56 62  2     1   3   2   0   3   3   2   0   0   3   0
## 3   39      0  0  0  3     0   3   2   3   0   2   2   1   3   2   3
##   f1k f1l f1m f1n f1o f1p f1q f1r f1s f1t
## 1   3   0   1   2   2   2   2   3   3   2
## 2   3   0   0   3   0   0   0   2   0   0
## 3   1   0   1   3   2   0   0   3   2   0

對newds做解釋與命名，並儲存

comment(newds) = "HELP baseline dataset"
comment(newds)

## [1] "HELP baseline dataset"

save(ds, file="savedfile")

讀取資料

write.csv(ds, file="ds.csv")

載入foreign套件，讀取檔案

library(foreign)
write.foreign(newds, "file.dat", "file.sas", package="SAS")

讀取newds檔案前十筆資料

with(newds, cesd[1:10])

##  [1] 49 30 39 15 39  6 52 32 50 46

with(newds, head(cesd, 10))

##  [1] 49 30 39 15 39  6 52 32 50 46

讀取newds檔案大於56的資料

with(newds, cesd[cesd > 56])

## [1] 57 58 57 60 58 58 57

載入dplyr套件，挑選大於56的資料，並顯示id

library(dplyr)
filter(newds, cesd > 56) %>% select(id, cesd)

##    id cesd
## 1  71   57
## 2 127   58
## 3 200   57
## 4 228   60
## 5 273   58
## 6 351   58
## 7  13   57

排序前四筆資料、最小值在第199個位置

with(newds, sort(cesd)[1:4])

## [1] 1 3 3 4

with(newds, which.min(cesd))

## [1] 199

載入mosaic套件，分組計數f1g的缺失值，對f1g做描述性統計

library(mosaic)

## Loading required package: lattice

## Loading required package: ggformula

## Loading required package: ggplot2

## Loading required package: ggstance

## 
## Attaching package: 'ggstance'

## The following objects are masked from 'package:ggplot2':
## 
##     geom_errorbarh, GeomErrorbarh

## 
## New to ggformula?  Try the tutorials: 
##  learnr::run_tutorial("introduction", package = "ggformula")
##  learnr::run_tutorial("refining", package = "ggformula")

## Loading required package: mosaicData

## Loading required package: Matrix

## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2

## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Note: If you use the Matrix package, be sure to load it BEFORE loading mosaic.
## 
## Have you tried the ggformula package for your plots?

## 
## Attaching package: 'mosaic'

## The following object is masked from 'package:Matrix':
## 
##     mean

## The following object is masked from 'package:ggplot2':
## 
##     stat

## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally

## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median,
##     prop.test, quantile, sd, t.test, var

## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum

tally(~ is.na(f1g), data=newds)

## is.na(f1g)
##  TRUE FALSE 
##     1   452

favstats(~ f1g, data=newds)

##  min Q1 median Q3 max mean  sd   n missing
##    0  1      2  3   3 1.73 1.1 452       1

編碼cesditems

cesditems = with(newds, cbind(f1a, f1b, f1c, (3 - f1d), f1e, f1f, f1g, 
   (3 - f1h), f1i, f1j, f1k, (3 - f1l), f1m, f1n, f1o, (3 - f1p), 
   f1q, f1r, f1s, f1t))

計算na個數

nmisscesd = apply(is.na(cesditems), 1, sum)

將cesditems設定為0，計算扣除後的總數，設定imputemeancesd的計算方式

ncesditems = cesditems
ncesditems[is.na(cesditems)] = 0
newcesd = apply(ncesditems, 1, sum)
imputemeancesd = 20/(20-nmisscesd)*newcesd

列出上述資料表格

data.frame(newcesd, newds$cesd, nmisscesd, imputemeancesd)[nmisscesd>0,]

##     newcesd newds.cesd nmisscesd imputemeancesd
## 4        15         15         1           15.8
## 17       19         19         1           20.0
## 87       44         44         1           46.3
## 101      17         17         1           17.9
## 154      29         29         1           30.5
## 177      44         44         1           46.3
## 229      39         39         1           41.1

不希望跑出訊息，載入套件

message=FALSE
library(dplyr)
library(memisc)

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## 
## Attaching package: 'memisc'

## The following object is masked from 'package:Matrix':
## 
##     as.array

## The following object is masked from 'package:ggplot2':
## 
##     syms

## The following objects are masked from 'package:dplyr':
## 
##     collect, recode, rename, syms

## The following objects are masked from 'package:stats':
## 
##     contr.sum, contr.treatment, contrasts

## The following object is masked from 'package:base':
## 
##     as.array

新增&定義abstinent、moderate、highrisk

newds = mutate(newds, drinkstat= 
  cases(
    "abstinent" = i1==0,
    "moderate" = (i1>0 & i1<=1 & i2<=3 & female==1) |
               (i1>0 & i1<=2 & i2<=4 & female==0),
    "highrisk" = ((i1>1 | i2>3) & female==1) |
               ((i1>2 | i2>4) & female==0)))

不要顯示程式碼，載入套件

echo=FALSE
library(mosaic)

不要顯示程式碼，卸除套件

echo=FALSE
detach(package:memisc)
detach(package:MASS)

載入套件，選擇變項，檢視365~370筆資料

library(dplyr)
tmpds = select(newds, i1, i2, female, drinkstat)
tmpds[365:370,]

##     i1 i2 female drinkstat
## 365  6 24      0  highrisk
## 366  6  6      0  highrisk
## 367  0  0      0 abstinent
## 368  0  0      1 abstinent
## 369  8  8      0  highrisk
## 370 32 32      0  highrisk

載入套件，選擇飲酒程度moderate的女性

library(dplyr)
filter(tmpds, drinkstat=="moderate" & female==1)

##   i1 i2 female drinkstat
## 1  1  1      1  moderate
## 2  1  3      1  moderate
## 3  1  2      1  moderate
## 4  1  1      1  moderate
## 5  1  1      1  moderate
## 6  1  1      1  moderate
## 7  1  1      1  moderate

不希望跑出訊息，載入套件，呈現飲酒程度次數分配表

message=FALSE
library(gmodels)
with(tmpds, CrossTable(drinkstat))

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  453 
## 
##  
##           | abstinent |  moderate |  highrisk | 
##           |-----------|-----------|-----------|
##           |        68 |        28 |       357 | 
##           |     0.150 |     0.062 |     0.788 | 
##           |-----------|-----------|-----------|
## 
## 
## 
##

呈現飲酒程度與性別的次數分配表

with(tmpds, CrossTable(drinkstat, female, 
  prop.t=FALSE, prop.c=FALSE, prop.chisq=FALSE))

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |           N / Row Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  453 
## 
##  
##              | female 
##    drinkstat |         0 |         1 | Row Total | 
## -------------|-----------|-----------|-----------|
##    abstinent |        42 |        26 |        68 | 
##              |     0.618 |     0.382 |     0.150 | 
## -------------|-----------|-----------|-----------|
##     moderate |        21 |         7 |        28 | 
##              |     0.750 |     0.250 |     0.062 | 
## -------------|-----------|-----------|-----------|
##     highrisk |       283 |        74 |       357 | 
##              |     0.793 |     0.207 |     0.788 | 
## -------------|-----------|-----------|-----------|
## Column Total |       346 |       107 |       453 | 
## -------------|-----------|-----------|-----------|
## 
##

創立“Male”,“Female”變項，計算個數

newds = transform(newds, 
  gender=factor(female, c(0,1), c("Male","Female")))
tally(~ female + gender, margin=FALSE, data=newds)

##       gender
## female Male Female
##      0  346      0
##      1    0    107

載入套件，排序資料，顯示1~5筆

library(dplyr)
newds = arrange(ds, cesd, i1)
newds[1:5, c("cesd", "i1", "id")]

##   cesd i1  id
## 1    1  3 233
## 2    3  1 139
## 3    3 13 418
## 4    4  4 251
## 5    4  9  95

載入套件，選出女性，計算cesd的平均

library(dplyr)
females = filter(ds, female==1)
with(females, mean(cesd))

## [1] 36.9

# an alternative approach
mean(ds$cesd[ds$female==1])

## [1] 36.9

計算female(包含Male, Female)的cesd平均

with(ds, tapply(cesd, female, mean))

##    0    1 
## 31.6 36.9

library(mosaic)
mean(cesd ~ female, data=ds)

##    0    1 
## 31.6 36.9

2020-03-30-Homework-6

Chang, C. Y.

2020-04-09

Comments on R script

不要顯示程式碼、程式碼要執行

環境設定：超過一行的字串用" "、顯示三個位數、控制顯示寬度

讀取檔案、使用dplyr、選取變項

顯示檔案變項名稱、顯示檔案前十個變項的資料

描述前十個變項：最大值、最小值、平均數、中位數、第一、三四分位數

從頭顯示前三筆資料

對newds做解釋與命名，並儲存

讀取資料

載入foreign套件，讀取檔案

讀取newds檔案前十筆資料

讀取newds檔案大於56的資料

載入dplyr套件，挑選大於56的資料，並顯示id

排序前四筆資料、最小值在第199個位置

載入mosaic套件，分組計數f1g的缺失值，對f1g做描述性統計

編碼cesditems

計算na個數

將cesditems設定為0，計算扣除後的總數，設定imputemeancesd的計算方式

列出上述資料表格

不希望跑出訊息，載入套件

新增&定義abstinent、moderate、highrisk

不要顯示程式碼，載入套件

不要顯示程式碼，卸除套件

載入套件，選擇變項，檢視365~370筆資料

載入套件，選擇飲酒程度moderate的女性

不希望跑出訊息，載入套件，呈現飲酒程度次數分配表

呈現飲酒程度與性別的次數分配表

創立“Male”,“Female”變項，計算個數

載入套件，排序資料，顯示1~5筆

載入套件，選出女性，計算cesd的平均

計算female(包含Male, Female)的cesd平均