Wednesday, September 10, 2014

R 语言与简单的回归分析

回归模型是计量里最基础也最常见的模型之一。究其原因,我想是因为在实际问题中我们并不知道总体分布如何,而且只有一组数据,那么试着对数据作回归分析将会是一个不错的选择。
一、简单线性回归
           简单的线性回归涉及到两个变量:一个是解释变量,通常称为x;另一个是被解释变量,通常称为y。回归会用常见的最小二乘算法拟合线性模型:
yi = β0 + β1xi +εi
其中β0和β1是回归系数,εi表示误差。
在R中,你可以通过函数lm()去计算他。Lm()用法如下:
lm(formula, data, subset, weights, na.action,
  method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
  singular.ok = TRUE, contrasts = NULL, offset, ...)

      参数是formula模型公式,例如y ~ x。公式中波浪号(~)左侧的是响应变量,右侧是预测变量。函数会估计回归系数β0和β1,分别以截距(intercept)和x的系数表示。
      有三种方式可以实现最小二乘法的简单线性回归,假设数据wage1(可以通过names函数查看数据框各项名称)
(1)lm(wage1$wage ~ wage1$educ + wage1$exper)
(2)lm (wage ~ educ + exper, data= wage1)
(3)attach(wage1)
    lm(wage~educ+exper)#不要忘记处理完后用detach()解出关联
         我们以数据wage1为例,可以看到工资与教育水平的线性关系:
运行下列代码:
library(foreign)
A<-read.dta("D:/R/data/WAGE1.dta")#导入数据
lm(wage~educ,data=A)
>lm(wage~educ,data=A)
Call:
lm(formula = wage~ educ, data = A)
Coefficients:
(Intercept)         educ 
-0.9049      0.5414
           当然得到这些数据是不够的,我们必须要有足够的证据去证明我们所做的回归的合理性。那么如何获取回归的信息呢?
          尝试运行以下代码:
result<-lm(wage~educ,data=A)
summary(result)
我们可以得到以下结果:
Call:
lm(formula = wage~ educ, data = A)
Residuals:
    Min     1Q     Median     3Q     Max
-5.3396   -2.1501    -0.9674     1.1921    16.6085
Coefficients:
            Estimate   Std.Error    t value  Pr(>|t|)   
(Intercept)   -0.90485   0.68497   -1.321   0.187   
educ         0.54136    0.05325 10.167   <2e-16 ***
---
Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1
Residual standarderror: 3.378 on 524 degrees of freedom
MultipleR-squared: 0.1648,     AdjustedR-squared: 0.1632
F-statistic: 103.4on 1 and 524 DF,   p-value: < 2.2e-16
              解读上述结果,我们不难看出,单从判决系数R-squared上看,回归结果是不理想的,但是,从p值来看,我们还是可以得到回归系数是很显著地(注意,这里的P<0.05就可以认为拒绝回归系数为0,即回归变量与被解释变量无关的原择假设,选择备择假设)所以说我们的回归的效果不好但还是可以接受的。当然,这一点也可以通过做散点图给我们直观的印象:
             但是影响薪酬的因素不只是education,可能还有其他的,比如工作经验,工作任期。为了更好地解释影响薪酬的因素,我们就必须用到多元线性回归。
             如果不需要那么多的信息,比如我们只需要F统计量,那么运行anova(result)即可得到方差分析表,更适合阅读。
            另外,再介绍一下函数predict,用法如下:
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,
        interval = c("none", "confidence", "prediction"),
        level = 0.95, type = c("response", "terms"),
        terms = NULL, na.action = na.pass,
        pred.var = res.var/weights, weights = 1, ...)
      说明一下,newdata的数据结构是一个数据框。
           这里还值得一提的参数时interval,他有三个选项:none代表不作区间预测,仅给出响应变量的估计;confidence是给出E(Y|X=x)的置信区间;prediction是给出真实Y的置信区间,运行代码你就会发现两者的差别。出于稳健性考虑,给出自变量取值时,预测Y值通常会采用prediction参数,但是对于X=x时Y的均值的预测就应该用confidence参数(因为回归模型是Y=βX+e,所以自变量相同,响应变量也未必一样)
二、多元线性回归
               还是使用lm函数。在公式的右侧指定多个预测变量,用加号(+)连接:
> lm(y ~ u + v+ w)
                显然,多元线性回归是简单的线性回归的扩展。可以有多个预测变量,还是用OLS计算多项式的系数。三变量的回归等同于这个线性模型:
yi = β0 + β1ui +β2vi + β3wi + εi
              在R中,简单线性回归和多元线性回归都是用lm函数。只要在模型公式的右侧增加变量即可。输出中会有拟合的模型的系数:
>result1<-lm(wage~educ+exper+tenure,data=A)
>summary(result1)
Call:
lm(formula = wage~ educ + exper + tenure, data = A)
Residuals:
    Min     1Q    Median      3Q    Max
-866.29    -249.23   -51.07   189.62   2190.01
Coefficients:
            Estimate    Std.Error    t value   Pr(>|t|)   
(Intercept)    -276.240   106.702   -2.589   0.009778 **
educ          74.415      6.287  11.836   <2e-16 ***
exper         14.892      3.253  4.578   5.33e-06 ***
tenure         8.257      2.498  3.306   0.000983 ***
---
Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1
Residual standarderror: 374.3 on 931 degrees of freedom
MultipleR-squared: 0.1459,     AdjustedR-squared: 0.1431
F-statistic:    53 on 3 and 931 DF,   p-value: < 2.2e-16
           我们将数据稍作平稳化处理,将wage换成log(wage),再来看看。
>plot(wage~educ,data=A)
>A$logwage<-log(A$wage)
>result1<-lm(logwage~educ+exper+tenure,data=A)
>summary(result1)
Call:
lm(formula =logwage ~ educ + exper + tenure, data = A)
Residuals:
     Min      1Q    Median      3Q      Max
-2.05802     -0.29645   -0.03265  0.28788   1.42809
Coefficients:
            Estimate   Std. Error   t value   Pr(>|t|)   
(Intercept)   0.284360  0.104190   2.729   0.00656**
educ        0.092029   0.007330 12.555  < 2e-16 ***
exper       0.004121   0.001723  2.391  0.01714 * 
tenure      0.022067  0.003094   7.133 3.29e-12 ***
---
Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1
Residual standarderror: 0.4409 on 522 degrees of freedom
MultipleR-squared: 0.316,      AdjustedR-squared: 0.3121
F-statistic: 80.39on 3 and 522 DF,   p-value: < 2.2e-16
             看得出,平稳化后的数据线性性是更加好的。
            下面我们来提取回归分析的各项统计数据:
           一些统计量和参数都被存储在lm或者summary中
output <-summary (result1)
SSR<- deviance(result1)#残差平方和;(另一种方法:RSquared <- output$r.squared)
LL<-logLik(result1) #对数似然统计量
DegreesOfFreedom<-result1$df #自由度
Yhat<- result1$fitted.values#拟合值向量
Resid<- result1$residuals
s<-output$sigma #误差标准差的估计值(假设同方差)
CovMatrix <-s^2*output$cov #系数的方差-协方差矩阵(与vcov(result1)同)
ANOVA<-anova(result1)#F统计量
confidentinterval<-confint(result1)#回归系数的置信区间,level=0.95
effects(result1)#计算正交效应向量(Vector of orthogonal effects )
predict函数用法与一元完全相同
三、检查结果
      检查回归结果是一件复杂而痛苦地事情,需要检验的东西也很多,当然有不少事是应该在数据进行回归分析之前就该处理的,比如检查复共线性;也有处理中需要考虑的,比如模型的选择,数据的变换;也有事后需要做的,比如残差正态性检验;还有需要关注与特别处理的数据,比如离群点,杠杆点。
     这里我们只提最简单与最常见的事后处理的基本分析。
    通过图形我们可以以一种十分直观的办法检测我们的拟合效果:
plot(result1)
     通过扩展包car中的函数来检测异常值与主要影响因子。代码如下:
library(car)
outlierTest(result1)
influence.measures(result1)
      在R中,线性回归计算变得无比简单,一个lm函数(或glm函数)基本上就摆平了OLS的一切。但拟合数据还仅仅是万里长征第一步。最终决定成败的是拟合的模型是否能真正地派上用场。这样对结果的检测与分析就显得尤为重要。

R语言 apply函数家族详解





apply {base}
通过对数组或者矩阵的一个维度使用函数生成值得列表或者数组、向量。
apply(X, MARGIN, FUN, ...)
X 阵列,包括矩阵
MARGIN  1表示矩阵行,2表示矩阵列,也可以是c(1,2)
例:
>xxx<-matrix(1:20,ncol=4)
>apply(xxx,1,mean)
[1]  8.5  9.5 10.5 11.5 12.5
>apply(xxx,2,mean)
[1]  3  8 13 18
>xxx
     [,1] [,2] [,3] [,4]
[1,]    1    6   11   16
[2,]    2    7   12   17
[3,]    3    8   13   18
[4,]    4    9   14   19
[5,]    5   10   15   20

lapply {base}
通过对x的每一个元素运用函数,生成一个与元素个数相同的值列表
lapply(X, FUN, ...)
X表示一个向量或者表达式对象,其余对象将被通过as.list强制转换为list
例:
> x <- list(a = 1:10, beta = exp(-3:3), logic = c(TRUE,FALSE,FALSE,TRUE))
> x
$a
 [1]  1  2  3  4  5  6  7  8  9 10
$beta
[1]  0.04978707  0.13533528  0.36787944  1.00000000  2.71828183  7.38905610
[7] 20.08553692
$logic
[1]  TRUE FALSE FALSE  TRUE
> lapply(x,mean)
$a
[1] 5.5
$beta
[1] 4.535125
$logic
[1] 0.5

sapply {base}
这是一个用户友好版本,是lapply函数的包装版。该函数返回值为向量、矩阵,如果simplify=”array”,且合适的情况下,将会通过simplify2array()函数转换为阵列。sapply(x, f, simplify=FALSE, USE.NAMES=FALSE)返回的值与lapply(x,f)是一致的。
sapply(X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)
X表示一个向量或者表达式对象,其余对象将被通过as.list强制转换为list
simplify 逻辑值或者字符串,如果可以,结果应该被简化为向量、矩阵或者高维数组。必须是命名的,不能是简写。默认值是TRUE,若合适将会返回一个向量或者矩阵。如果simplify=”array”,结果将返回一个阵列。
USE.NAMES  逻辑值,如果为TRUE,且x没有被命名,则对x进行命名。
例:
k<-c("a","b","c")
> sapply(k, paste,USE.NAMES=FALSE,1:5,sep="...")
     [,1]    [,2]    [,3]
[1,] "a...1" "b...1" "c...1"
[2,] "a...2" "b...2" "c...2"
[3,] "a...3" "b...3" "c...3"
[4,] "a...4" "b...4" "c...4"
[5,] "a...5" "b...5" "c...5"
> sapply(k, paste,USE.NAMES=TRUE,1:5,sep="...")
     a       b       c    
[1,] "a...1" "b...1" "c...1"
[2,] "a...2" "b...2" "c...2"
[3,] "a...3" "b...3" "c...3"
[4,] "a...4" "b...4" "c...4"
[5,] "a...5" "b...5" "c...5"
> sapply(k, paste,USE.NAMES=TRUE,1:5,sep="...",simplyfy=TRUE)
     a              b              c          
[1,] "a...1...TRUE" "b...1...TRUE" "c...1...TRUE"
[2,] "a...2...TRUE" "b...2...TRUE" "c...2...TRUE"
[3,] "a...3...TRUE" "b...3...TRUE" "c...3...TRUE"
[4,] "a...4...TRUE" "b...4...TRUE" "c...4...TRUE"
[5,] "a...5...TRUE" "b...5...TRUE" "c...5...TRUE"
> sapply(k, paste,simplify=TRUE,USE.NAMES=TRUE,1:5,sep="...")
     a       b       c    
[1,] "a...1" "b...1" "c...1"
[2,] "a...2" "b...2" "c...2"
[3,] "a...3" "b...3" "c...3"
[4,] "a...4" "b...4" "c...4"
[5,] "a...5" "b...5" "c...5"
> sapply(k, paste,simplify=FALSE,USE.NAMES=TRUE,1:5,sep="...")
$a
[1] "a...1" "a...2" "a...3" "a...4" "a...5"
$b
[1] "b...1" "b...2" "b...3" "b...4" "b...5"
$c
[1] "c...1" "c...2" "c...3" "c...4" "c...5"


vapply {base}
vapply类似于sapply函数,但是它的返回值有预定义类型,所以它使用起来会更加安全,有的时候会更快
在vapply函数中总是会进行简化,vapply会检测FUN的所有值是否与FUN.VALUE兼容,以使他们具有相同的长度和类型。类型顺序:逻辑<</span>整型<</span>实数<</span>复数
vapply(X, FUN, FUN.VALUE, ..., USE.NAMES = TRUE)
X表示一个向量或者表达式对象,其余对象将被通过as.list强制转换为list
simplify 逻辑值或者字符串,如果可以,结果应该被简化为向量、矩阵或者高维数组。必须是命名的,不能是简写。默认值是TRUE,若合适将会返回一个向量或者矩阵。如果simplify=”array”,结果将返回一个阵列。
USE.NAMES  逻辑值,如果为TRUE,且x没有被命名,则对x进行命名。
FUN.VALUE   一个通用型向量,FUN函数返回值得模板
例:
> x<-data.frame(a=rnorm(4,4,4),b=rnorm(4,5,3),c=rnorm(4,5,3))
> vapply(x,mean,c(c=0))
         a          b          c
 1.8329043  6.0442858 -0.1437202
> k<-function(x)
+ {
+ list(mean(x),sd(x))
+ }
> vapply(x,k,c(c=0))
错误于vapply(x, k, c(c = 0)) : 值的长度必需为1,
 但FUN(X[[1]])结果的长度却是2
> vapply(x,k,c(c=0,b=0))
错误于vapply(x, k, c(c = 0, b = 0)) : 值的种类必需是'double',
 但FUN(X[[1]])结果的种类却是'list'
> vapply(x,k,c(list(c=0,b=0)))
  a        b        c      
c 1.832904 6.044286 -0.1437202
b 1.257834 1.940433 3.649194

tapply {base}
对不规则阵列使用向量,即对一组非空值按照一组确定因子进行相应计算
tapply(X, INDEX, FUN, ..., simplify = TRUE)
x  一个原子向量,典型的是一个向量
INDEX  因子列表,和x长度一样,元素将被通过as.factor强制转换为因子
simplify  若为FALSE,tapply将以列表形式返回阵列。若为TRUE,FUN则直接返回数值
例:
> height <- c(174, 165, 180, 171, 160)
> sex<-c("F","F","M","F","M")
> tapply(height, sex, mean)
  F     M
170   170

eapply {base}
eapply函数通过对environment中命名值进行FUN计算后返回一个列表值,用户可以请求所有使用过的命名对象。
eapply(env, FUN, ..., all.names = FALSE, USE.NAMES = TRUE)
env  将被使用的环境
all.names  逻辑值,指示是否对所有值使用该函数
USE.NAMES  逻辑值,指示返回的列表结果是否包含命名
例:
> require(stats)
>
> env <- new.env(hash = FALSE) # so the order is fixed
> env$a <- 1:10
> env$beta <- exp(-3:3)
> env$logic <- c(TRUE, FALSE, FALSE, TRUE)
> # what have we there?
> utils::ls.str(env)
a :  int [1:10] 1 2 3 4 5 6 7 8 9 10
beta :  num [1:7] 0.0498 0.1353 0.3679 1 2.7183 ...
logic :  logi [1:4] TRUE FALSE FALSE TRUE
>
> # compute the mean for each list element
>        eapply(env, mean)
$logic
[1] 0.5

$beta
[1] 4.535125

$a
[1] 5.5

> unlist(eapply(env, mean, USE.NAMES = FALSE))
[1] 0.500000 4.535125 5.500000
>
> # median and quartiles for each element (making use of "..." passing):
> eapply(env, quantile, probs = 1:3/4)
$logic
25% 50% 75%
0.0 0.5 1.0

$beta
      25%       50%       75%
0.2516074 1.0000000 5.0536690

$a
 25%  50%  75%
3.25 5.50 7.75

> eapply(env, quantile)
$logic
  0%  25%  50%  75% 100%
 0.0  0.0  0.5  1.0  1.0

$beta
         0%         25%         50%         75%        100%
 0.04978707  0.25160736  1.00000000  5.05366896 20.08553692

$a
   0%   25%   50%   75%  100%
 1.00  3.25  5.50  7.75 10.00

mapply {base}
mapply是sapply的多变量版本。将对...中的每个参数运行FUN函数,如有必要,参数将被循环。
mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE)
MoreArgs   FUN函数的其他参数列表
SIMPLIFY   逻辑或者字符串,可以减少结果成为一个向量、矩阵或者更高维阵列,详见sapply的simplify参数
USE.NAMES   逻辑值,如果第一个参数...已被命名,将使用这个字符向量作为名字
例:
> mapply(rep, 1:4, 4:1)
[[1]]
[1] 1 1 1 1

[[2]]
[1] 2 2 2

[[3]]
[1] 3 3

[[4]]
[1] 4

rapply {base}
rapply是lapply的递归版本
rapply(X, FUN, classes = "ANY", deflt = NULL, how = c("unlist", "replace", "list"), ...)
X  一个列表
classes  关于类名的字符向量,或者为any时则匹配任何类
deflt  默认结果,如果使用了how=”replace”,则不能使用
how  字符串匹配三种可能结果

Tuesday, September 9, 2014

一些问题的理解(1)

1.NA与NULL
x <- c(88,NA,12,168,13)
x
[188 NA 12 168 13
mean(x)
[1NA
mean(x,na.rm=T)#na.rm=T是指NA,remove,TRUE大家理解吧.
[170.25
x <- c(88,NULL,12,168,13)
mean(x)
[170.25
NA指缺失值,NULL指不存在的

2.
i <- 2
1:i-1 # 注意这个地方!是指(1:i)-1
[10 1
1:(i-1)
[11

3.which()

which可以很方便的用来寻找一些符合固定条件的下标,比如我想找x=c(1,2,3,4,5,1,8)这里面第2个等于1的位置:
which(x==1)[2]

Factor variables

1. Creating factor variables

Factor variables are categorical variables that can be either numeric or string variables. There are a number of advantages to converting categorical variables to factor variables. Perhaps the most important advantage is that they can be used in statistical modeling where they will be implemented correctly, i.e., they will then be assigned the correct number of degrees of freedom. Factor variables are also very useful in many different types of graphics. Furthermore, storing string variables as factor variables is a more efficient use of memory. To create a factor variable we use the factor function. The only required argument is a vector of values which can be either string or numeric. Optional arguments include the levels argument, which determines the categories of the factor variable, and the default is the sorted list of all the distinct values of the data vector. The labels argument is another optional argument which is a vector of values that will be the labels of the categories in the levels argument. The exclude argument is also optional; it defines which levels will be classified as NA in any output using the factor variable.
First we will generate a vector of numeric data called schtyp. It involves the random number generator so we will set the seed to equal 124 in order to make the results reproducible.
set.seed(124)
schtyp <- sample(0:1, 20, replace = TRUE)
schtyp
##  [1] 0 0 1 0 0 0 1 0 1 0 1 1 1 1 0 0 1 1 1 0
is.factor(schtyp)
## [1] FALSE
is.numeric(schtyp)
## [1] TRUE
Now let's create a factor variable called schtyp.f based on schtyp. The first label, private, will correspond to schtyp=0 and the second label, public, will correspond to schtyp=1 because the order of the labels will follow the numeric order of the data.
schtyp.f <- factor(schtyp, labels = c("private", "public"))
schtyp.f
##  [1] private private public  private private private public  private
##  [9] public  private public  public  public  public  private private
## [17] public  public  public  private
## Levels: private public
is.factor(schtyp.f)
## [1] TRUE
Let's generate a string variable called ses (socio-economic status).
ses <- c("low", "middle", "low", "low", "low", "low", "middle", "low", "middle",
    "middle", "middle", "middle", "middle", "high", "high", "low", "middle",
    "middle", "low", "high")

is.factor(ses)
## [1] FALSE
is.character(ses)
## [1] TRUE
Creating a factor variable ses.f.bad.order based on ses.
ses.f.bad.order <- factor(ses)
is.factor(ses.f.bad.order)
## [1] TRUE
levels(ses.f.bad.order)
## [1] "high"   "low"    "middle"
The problem is that the levels are ordered according to the alphabetical order of the categories of ses. Thus, "high" is the lowest level of ses.f.bad.order, "middle" is the middle level and "low" is the highest level. In order to fix the ordering we need to use the levels argument to indicate the correct ordering of the categories. Let's create a new factor variable called ses.f with the correct order of categories.
ses.f <- factor(ses, levels = c("low", "middle", "high"))
is.factor(ses.f)
## [1] TRUE
levels(ses.f)
## [1] "low"    "middle" "high"

2. Creating ordered factor variables

We can create ordered factor variables by using the function ordered. This function has the same arguments as the factor function. Let's create an ordered factor variable called ses.order based on the variable ses created in the above example.
ses.order <- ordered(ses, levels = c("low", "middle", "high"))
ses
##  [1] "low"    "middle" "low"    "low"    "low"    "low"    "middle"
##  [8] "low"    "middle" "middle" "middle" "middle" "middle" "high"  
## [15] "high"   "low"    "middle" "middle" "low"    "high"
ses.order
##  [1] low    middle low    low    low    low    middle low    middle
## [10] middle middle middle middle high   high   low    middle middle
## [19] low    high  
## Levels: low < middle < high
is.factor(ses.order)
## [1] TRUE

3. Adding and dropping levels in factor variables

Below we will add an element from a new level ("very.high") to ses.f our existing factor variable, ses.f. The number in the square brackets ( [21] ) indicates the number of the element whose label we wish to change.
ses.f[21] <- "very.high"
## Warning: invalid factor level, NA generated
ses.f
##  [1] low    middle low    low    low    low    middle low    middle
## [10] middle middle middle middle high   high   low    middle middle
## [19] low    high     
## Levels: low middle high
We can see that instead of changing from "high" to "very.high", the label was changed from "high" to <NA>. To do this correctly, we need to first add the new level, "very.high", to the factor variable ses.f which we do by using the factor function with the levels argument. Then we can finally add an element to the factor variable from the new level.
ses.f <- factor(ses.f, levels = c(levels(ses.f), "very.high"))
ses.f[21] <- "very.high"
ses.f
##  [1] low       middle    low       low       low       low      
##  [7] middle    low       middle    middle    middle    middle   
## [13] middle    high      high      low       middle    middle   
## [19] low       high      very.high
## Levels: low middle high very.high
levels(ses.f)
## [1] "low"       "middle"    "high"      "very.high"
Dropping a level of a factor variable is a little easier. The simplest way is to first remove all the elements within the level to be removed and then to redeclare the variable to be a factor variable. (The level is not automatically removed if there are no elements in it because we could just by chance have a sample which did not contain elements from a specific level.) Let's illustrate this by removing the level of "very.high" from the ses.f variable.
ses.f.new <- ses.f[ses.f != "very.high"]
ses.f.new
##  [1] low    middle low    low    low    low    middle low    middle
## [10] middle middle middle middle high   high   low    middle middle
## [19] low    high  
## Levels: low middle high very.high
ses.f.new <- factor(ses.f.new)
ses.f.new
##  [1] low    middle low    low    low    low    middle low    middle
## [10] middle middle middle middle high   high   low    middle middle
## [19] low    high  
## Levels: low middle high
levels(ses.f.new)
## [1] "low"    "middle" "high"

4. Examples of the usefulness of factor variables

To illustrate the usefulness of factor variables we are first going to create a data frame with all the variables we have used in the previous examples, plus an additional continuous variable called read which contains the reading scores. We also redefine ses.f to equal the ses.f.new variable which does not have any "very.high" elements.
ses.f <- ses.f.new
read <- c(34, 39, 63, 44, 47, 47, 57, 39, 48, 47, 34, 37, 47, 47, 39, 47,
    47, 50, 28, 60)

# combining all the variables in a data frame
combo <- data.frame(schtyp, schtyp.f, ses, ses.f, read)
Tables are much easier to interpret when using factor variables because they add useful labels to the table and they arrange the factors in a more understandable order.
table(ses, schtyp)
##         schtyp
## ses      0 1
##   high   2 1
##   low    6 2
##   middle 2 7
table(ses.f, schtyp.f)
##         schtyp.f
## ses.f    private public
##   low          6      2
##   middle       2      7
##   high         2      1
Graphics are another area that benefits from the use of factor variables. As in the tables the factor variable will indicate a better ordering of the graphs as well as add useful labels.
library(lattice)
bwplot(schtyp ~ read | ses, data = combo, layout = c(2, 2))
plot of chunk unnamed-chunk-14
bwplot(schtyp.f ~ read | ses.f, data = combo, layout = c(2, 2))
plot of chunk unnamed-chunk-14

Subsetting Data

1. Subsetting variables

To manipulate data frames in R we can use the bracket notation to access the indices for the observations and the variables. It is easiest to think of the data frame as a rectangle of data where the rows are the observations and the columns are the variables. Just like in matrix algebra, the indices for a rectangle of data follow the RxC principle; in other words, the first index is for Rows and the second index is for Columns [R, C]. When we only want to subset variables (or columns) we use the second index and leave the first index blank. Leaving an index blank indicates that you want to keep all the elements in that dimension. In the first example we create the data frame hsb3 containing only the variables idread and write, but all the observations from the original data frame hsb2.small. In order to know which variables correspond to which number in the index we use the names function, which will list the names of the variables in the order in which they appear in the data frame. From this list we see that id is variable 1, read is variable 7 and write is variable 8. We cannot refer to the variables by their names alone until we have attached the data.
hsb2.small <- read.csv("http://www.ats.ucla.edu/stat/data/hsb2_small.csv")

# using the names function to see names of the variables and which column of
# data to which they correspond
names(hsb2.small)
##  [1] "id"      "female"  "race"    "ses"     "schtyp"  "prog"    "read"   
##  [8] "write"   "math"    "science" "socst"
(hsb3 <- hsb2.small[, c(1, 7, 8)])
##     id read write
## 1   70   57    52
## 2  121   68    59
## 3   86   44    33
## 4  141   63    44
## 5  172   47    52
## 6  113   44    52
## 7   50   50    59
## 8   11   34    46
## 9   84   63    57
## 10  48   57    55
## 11  75   60    46
## 12  60   57    65
## 13  95   73    60
## 14 104   54    63
## 15  38   45    57
## 16 115   42    49
## 17  76   47    52
## 18 195   57    57
## 19 114   68    65
## 20  85   55    39
## 21 167   63    49
## 22 143   63    63
## 23  41   50    40
## 24  20   60    52
## 25  12   37    44
If the variables we want are in consecutive columns, we can use the colon notation rather than list them using the c function. In the next example we create the data frame hsb4 containing the first four variables of hsb2.small.
(hsb4 <- hsb2.small[, 1:4])
##     id female race ses
## 1   70      0    4   1
## 2  121      1    4   2
## 3   86      0    4   3
## 4  141      0    4   3
## 5  172      0    4   2
## 6  113      0    4   2
## 7   50      0    3   2
## 8   11      0    1   2
## 9   84      0    4   2
## 10  48      0    3   2
## 11  75      0    4   2
## 12  60      0    4   2
## 13  95      0    4   3
## 14 104      0    4   3
## 15  38      0    3   1
## 16 115      0    4   1
## 17  76      0    4   3
## 18 195      0    4   2
## 19 114      0    4   3
## 20  85      0    4   2
## 21 167      0    4   2
## 22 143      0    4   2
## 23  41      0    3   2
## 24  20      0    1   3
## 25  12      0    1   2

2. Subsetting observations

We subset observations by also using the bracket notation but now we use the first index and leave the second index blank. This indicates that we want all the variables for specific observations. In the first example we create the data frame hsb5, which contains the first 10 observations of hsb2.small.
(hsb5 <- hsb2.small[1:10, ])
##     id female race ses schtyp prog read write math science socst
## 1   70      0    4   1      1    1   57    52   41      47    57
## 2  121      1    4   2      1    3   68    59   53      63    61
## 3   86      0    4   3      1    1   44    33   54      58    31
## 4  141      0    4   3      1    3   63    44   47      53    56
## 5  172      0    4   2      1    2   47    52   57      53    61
## 6  113      0    4   2      1    2   44    52   51      63    61
## 7   50      0    3   2      1    1   50    59   42      53    61
## 8   11      0    1   2      1    2   34    46   45      39    36
## 9   84      0    4   2      1    1   63    57   54      58    51
## 10  48      0    3   2      1    2   57    55   52      50    51
We can also subset observations based on logical tests. In the following example we create the data frame hsb6, which contains only the observations for which ses=1. For a logical equality we need to use the double equal sign notation. We also need to refer to the variable, ses in the data frame hsb2.small, which we do using $.
(hsb6 <- hsb2.small[hsb2.small$ses == 1, ])
##     id female race ses schtyp prog read write math science socst
## 1   70      0    4   1      1    1   57    52   41      47    57
## 15  38      0    3   1      1    2   45    57   50      31    56
## 16 115      0    4   1      1    1   42    49   43      50    56
In the previous example we used a logical test to subset the observations, but we only tested for one variable being equal to a single value. We can also subset using a logical test that will test a single variable being equal to the elements in a list, and we do this by using the %in% function. In the following example we create the data frame hsb7, which contains the observations where id is equal to 11, 12, 20, 48, 86 or 195.
(hsb7 <- hsb2.small[hsb2.small$id %in% c(12, 48, 86, 11, 20, 195), ])
##     id female race ses schtyp prog read write math science socst
## 3   86      0    4   3      1    1   44    33   54      58    31
## 8   11      0    1   2      1    2   34    46   45      39    36
## 10  48      0    3   2      1    2   57    55   52      50    51
## 18 195      0    4   2      2    1   57    57   60      58    56
## 24  20      0    1   3      1    2   60    52   57      61    61
## 25  12      0    1   2      1    3   37    44   45      39    46
It is also possible to combine logical tests. In the following example we create the data frame hsb8, which contains only the observations where ses=3 andfemale=0. Here to avoid having to type hsb2.small multiple times, we use the with function to let R know that it should look for ses and female inside thehsb2.small data frame.
(hsb8 <- hsb2.small[with(hsb2.small, ses == 3 & female == 0), ])
##     id female race ses schtyp prog read write math science socst
## 3   86      0    4   3      1    1   44    33   54      58    31
## 4  141      0    4   3      1    3   63    44   47      53    56
## 13  95      0    4   3      1    2   73    60   71      61    71
## 14 104      0    4   3      1    2   54    63   57      55    46
## 17  76      0    4   3      1    2   47    52   51      50    56
## 19 114      0    4   3      1    2   68    65   62      55    61
## 24  20      0    1   3      1    2   60    52   57      61    61
The subset function with a logical statement will let you subset the data frame by observations. In the following example the write.50 data frame contains only the observations for which the values of the variable write is greater than 50. Note that one convenient feature of the subset function, is R assumes variable names are within the data frame being subset, so there is no need to tell R where to look for write.
(write.50 <- subset(hsb2.small, write > 50))
##     id female race ses schtyp prog read write math science socst
## 1   70      0    4   1      1    1   57    52   41      47    57
## 2  121      1    4   2      1    3   68    59   53      63    61
## 5  172      0    4   2      1    2   47    52   57      53    61
## 6  113      0    4   2      1    2   44    52   51      63    61
## 7   50      0    3   2      1    1   50    59   42      53    61
## 9   84      0    4   2      1    1   63    57   54      58    51
## 10  48      0    3   2      1    2   57    55   52      50    51
## 12  60      0    4   2      1    2   57    65   51      63    61
## 13  95      0    4   3      1    2   73    60   71      61    71
## 14 104      0    4   3      1    2   54    63   57      55    46
## 15  38      0    3   1      1    2   45    57   50      31    56
## 17  76      0    4   3      1    2   47    52   51      50    56
## 18 195      0    4   2      2    1   57    57   60      58    56
## 19 114      0    4   3      1    2   68    65   62      55    61
## 22 143      0    4   2      1    3   63    63   75      72    66
## 24  20      0    1   3      1    2   60    52   57      61    61
There is no limit to how many logical statements may be combined to achieve the subsetting that is desired. The data frame write.1 contains only the observations for which the values of the variable write is greater than 50 and for which the variable read is greater than 60.
(write.1 <- subset(hsb2.small, write > 50 & read > 60))
##     id female race ses schtyp prog read write math science socst
## 2  121      1    4   2      1    3   68    59   53      63    61
## 9   84      0    4   2      1    1   63    57   54      58    51
## 13  95      0    4   3      1    2   73    60   71      61    71
## 19 114      0    4   3      1    2   68    65   62      55    61
## 22 143      0    4   2      1    3   63    63   75      72    66
It is possible to subset both rows and columns using the subset function. The select argument lets you subset variables (columns). The data frame write.2contains only the variables write and read and then only the observations of these two variables where the values of variable write are greater than 50 and the values of variable read are greater than 65.
(write.2 <- subset(hsb2.small, write > 50 & read > 60, select = c(write, read)))
##    write read
## 2     59   68
## 9     57   63
## 13    60   73
## 19    65   68
## 22    63   63
In the data frame write.3 contains only the observations in variables read through science for which the values in the variable science are less than 55.
(write.3 <- subset(hsb2.small, science < 55, select = read:science))
##    read write math science
## 1    57    52   41      47
## 4    63    44   47      53
## 5    47    52   57      53
## 7    50    59   42      53
## 8    34    46   45      39
## 10   57    55   52      50
## 11   60    46   51      53
## 15   45    57   50      31
## 16   42    49   43      50
## 17   47    52   51      50
## 20   55    39   57      53
## 25   37    44   45      39

3. Subsetting both variables and observations

We can subset variables and observations by simply combining the two above methods of subsetting. We accomplish this by subsetting using both indices at the same time. In the following example we create the data frame hsb9 in which we keep only the variables idfemaleraceses and read and only the observations where ses=3. Note again that because we are not using subset, we have to let R know where to find the variable ses by explicitly pointing tohsb2.small.
# using the names function to see names of the variables and which column of
# data to which they correspond
names(hsb2.small)
##  [1] "id"      "female"  "race"    "ses"     "schtyp"  "prog"    "read"   
##  [8] "write"   "math"    "science" "socst"
(hsb9 <- hsb2.small[hsb2.small$ses == 3, c(1:4, 7)])
##     id female race ses read
## 3   86      0    4   3   44
## 4  141      0    4   3   63
## 13  95      0    4   3   73
## 14 104      0    4   3   54
## 17  76      0    4   3   47
## 19 114      0    4   3   68
## 24  20      0    1   3   60