zoukankan      html  css  js  c++  java
  • R in action读书笔记(17)第十二章 重抽样与自助法

    12.4 置换检验点评

    除coin和lmPerm包外,R还提供了其他可做置换检验的包。perm包能实现coin包中的部分功能,因此可作为coin包所得结果的验证。corrperm包提供了有重复测量的相关性的置换检验。

    logregperm包提供了Logistic回归的置换检验。另外一个非常重要的包是glmperm,它涵盖了广义线性模型的置换检验依靠基础的抽样分布理论知识,置换检验提供了另外一个十分强大的可选检验思路。对于上面描述的每一种置换检验,我们完全可以在做统计假设检验时不理会正态分布、t分布、F分布或者卡方分布。当然,置换检验真正发挥功用的地方是处理非正态数据(如分布偏倚很大)、存在离群点、样本很小或无法做参数检验等情况。不过,如果初始样本对感兴趣的总体情况代表性很差,即使是置换检验也无法提高推断效果。置换检验主要用于生成检验零假设的p值,它有助于回答“效应是否存在”这样的问题。

    12.5 自助法

    所谓自助法,即从初始样本重复随机替换抽样,生成一个或一系列待检验统计量的经验分布。

    无需假设一个特定的理论分布,便可生成统计量的置信区间,并能检验统计假设。比如,你想计算一个样本均值95%的置信区间。假设均值的样本分布不是正态分布:

    (1) 从样本中随机选择10个观测,抽样后再放回。有些观测可能会被选择多次,有些可能一

    直都不会被选中。

    (2) 计算并记录样本均值。

    (3) 重复1和2一千次。

    (4) 将1000个样本均值从小到大排序。

    (5) 找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数,它们就限

    定了95%的置信区间。

    12.6 boot 包中的自助法

    boot包扩展了自助法和重抽样的相关用途。可以对一个统计量(如中位数)或一个统计

    量向量(如一列回归系数)使用自助法.

    自助法有三个主要步骤。

    (1) 写一个能返回待研究统计量值的函数。如果只有单个统计量(如中位数),函数应该返回

    一个数值;如果有一列统计量(如一列回归系数),函数应该返回一个向量。

    (2) 为生成R中自助法所需的有效统计量重复数,使用boot()函数对上面所写的函数进行处理。

    (3) 使用boot.ci()函数获取第(2)步生成的统计量的置信区间。

    主要的自助法函数是boot(),它的格式为:bootobject<-boot(data=,statistic=,R=,…)

    data:量、矩阵或者数据框

    statistic:生成k个统计量以供自举的函数(k=1时对单个统计量进行自助抽样)函数需包括indices参数,以便boot()函数用它从每个重复中选择实例

    R:自助抽样的次数

    ...:其他对生成待研究统计量有用的参数,可在函数中传输

    boot()函数调用统计量函数R次,每次都从整数1:nrow(data)中生成一列有放回的随机指

    标,这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算,结果存储在

    bootobject中。boot()函数中返回对象所含的元素

    t0 从原始数据得到的k个统计量的观测值

    t 一个R × k矩阵,每行即k个统计量的自助重复值

    你可以如bootobject$t0和bootobject$t这样来获取这些元素。

    一旦生成了自助样本,可通过print()和plot()来检查结果。如果结果看起来还算合理,

    使用boot.ci()函数获取统计量的置信区间。格式如下:

    Boot.ci(bootobject,couf=type=)

    bootobject boot()函数返回的对象

    conf 预期的置信区间(默认:conf =0.95)

    type 返回的置信区间类型。可能值为norm、basic、stud、perc、bca和all(默认:type =all)

    type参数设定了获取置信区间的方法。perc方法(分位数)展示的是样本均值,bca将根据

    偏差对区间做简单调整

    12.6.1 对单个统计量使用自助法

    > rsq<-function(formula,data,indices){

    + d<-data[indices,]

    + fit<-lm(formula,data=d)

    + return(summary(fit)$r.square)

    + }

    > set.seed(1234)

    > results<-boot(data=mtcars,statistic=rsq,R=1000,formula=mpg~wt+disp)

    > print(results)

    ORDINARY NONPARAMETRICBOOTSTRAP

    Call:

    boot(data = mtcars,statistic = rsq, R = 1000, formula = mpg ~

    wt + disp)

    Bootstrap Statistics :

    original bias std. error

    t1* 0.78093060.0133367 0.05068926

    可以看到,自助的R平方值不呈正态分布。它的95%的置信区间可以通过如下代

    码获得:

    > boot.ci(results,type=c("perc","bca"))

    BOOTSTRAP CONFIDENCE INTERVALCALCULATIONS

    Based on 1000 bootstrap replicates

    CALL :

    boot.ci(boot.out = results, type =c("perc", "bca"))

    Intervals :

    Level Percentile BCa

    95% ( 0.6838, 0.8833 ) ( 0.6344, 0.8549 )

    Calculations and Intervals on OriginalScale

    Some BCa intervals may be unstable

    12.6.2 多个统计量的自助法

    首先,创建一个返回回归系数向量的函数:

    > bs<-function(formula,data,indices){

    + d<-data[indices,]

    + fit<-lm(formula,data=d)

    + return(coef(fit))

    + }

    > results<-boot(data=mtcars,statistic=bs,R=1000,formula=mpg~wt+disp)

    > print(results)

    ORDINARY NONPARAMETRIC BOOTSTRAP

    Call:

    boot(data = mtcars, statistic = bs, R= 1000, formula = mpg ~

    wt + disp)

    Bootstrap Statistics :

    original bias std. error

    t1* 34.96055404 0.1101475088 2.445950503

    t2* -3.35082533 -0.08573079461.114315903

    t3* -0.01772474 0.0003223228 0.008201569

    plot(results,index=2)

    12.7 小结

    本章,我们介绍了一系列基于随机化和重抽样的计算机密集型方法,它们使你无需理论分布

    的知识便能够进行假设检验,获得置信区间。当数据来自未知分布,或者存在严重的离群点,或

    者样本量过小,又或者没有参数方法可以回答你感兴趣的假设问题时,这些方法是非常实用的。

  • 相关阅读:
    Servlet(九):web.xml文件和server.xml文件
    Servlet(八):ServletContext对象和ServletConfig对象
    Servlet(七):session
    Servlet(六):Cookie
    Servlet(五):请求转发和重定向
    Servlet(四):request和response对象
    Servlet(三):生命周期、常用方法、常见错误
    【php】在Windows2003下的IIS配置php5.4
    Spring AOP监控SQL运行
    算法导论—无向图的遍历(BFS+DFS,MATLAB)
  • 原文地址:https://www.cnblogs.com/jpld/p/4478006.html
Copyright © 2011-2022 走看看