精华内容
下载资源
问答
  • R语言数据降维

    2020-12-31 13:47:19
    今天是2020年的最后一天...目前R语言的几种降维方式。 首先需要配置数据 data<-matrix(rnorm(3000),ncol=6) colnames(data)=paste0('gene',1:6) rownames(data)=c(paste0(rep('C',nrow(data)/2),1:(nrow(data)/2)),

    今天是2020年的最后一天,这一年多灾多难,让我们唏嘘不已。不论2020如何,生活仍要继续,希望即将到来的2021年,可以‘牛’转乾坤。
    这也是2020年的最后一篇博客了,在这里给大家介绍一下。目前R语言的几种降维方式。
    首先需要配置数据

    data<-matrix(rnorm(3000),ncol=6)
    colnames(data)=paste0('gene',1:6)
    rownames(data)=c(paste0(rep('C',nrow(data)/2),1:(nrow(data)/2)),
                     paste0(rep('P',nrow(data)/2),1:(nrow(data)/2)))
    head(data)
    

    在这里插入图片描述1、PCA降维

    data.pca <- prcomp(t(data), scale. = TRUE)
    head(data.pca$rotation)
    

    在这里插入图片描述

    pca_data=data.pca$rotation[,1:2]
    pca_data=data.frame(sample=rownames(pca_data),
                        pca_data,
                        group=c(rep('C',nrow(data)/2),rep('P',nrow(data)/2)))
    head(pca_data)
    

    在这里插入图片描述

    #绘图
    library(ggplot2)
    ggplot(data=pca_data,aes(x=PC1,y=PC2))+
      geom_point(aes(colour=group,shape=group),size=2)+
      theme(text=element_text(size=10),legend.title=element_blank(),
            panel.background = element_rect(fill = "white", colour = "black",size = 0.2), 
            legend.key = element_rect(fill = "white", colour = "white"),
            legend.background = (element_rect(colour= "white",fill = "white")))
    

    在这里插入图片描述2、非度量MDS降维

    dis_data = dist(data,'euclidean',p=2)
    mds_x = cmdscale(dis_data,k=2)
    mds_x = data.frame(mds_x)
    colnames(mds_x)=c('MDS.1','MDS.2')
    head(mds_x)
    

    在这里插入图片描述

    mds_data=data.frame(sample=rownames(mds_x),
                        mds_x,
                        group=c(rep('C',nrow(data)/2),rep('P',nrow(data)/2)))
    head(mds_data)
    

    在这里插入图片描述

    library(ggplot2)
    ggplot(data=mds_data,aes(x=MDS.1,y=MDS.2))+
      geom_point(aes(colour=group,shape=group),size=2)+
      theme(text=element_text(size=10),legend.title=element_blank(),
            panel.background = element_rect(fill = "white", colour = "black",size = 0.2), 
            legend.key = element_rect(fill = "white", colour = "white"),
            legend.background = (element_rect(colour= "white",fill = "white")))
    

    在这里插入图片描述3、umap降维

    library(umap)
    data.umap = umap(data)
    
    umap_data=data.umap$layout
    colnames(umap_data)=c('umap_1','umap_2')
    head(umap_data)
    

    在这里插入图片描述

    umap_data=data.frame(sample=rownames(umap_data),
                        umap_data,
                        group=c(rep('C',nrow(data)/2),rep('P',nrow(data)/2)))
    head(umap_data)
    

    在这里插入图片描述

    library(ggplot2)
    ggplot(data=umap_data,aes(x=umap_1,y=umap_2))+
      geom_point(aes(colour=group,shape=group),size=2)+
      theme(text=element_text(size=10),legend.title=element_blank(),
            panel.background = element_rect(fill = "white", colour = "black",size = 0.2), 
            legend.key = element_rect(fill = "white", colour = "white"),
            legend.background = (element_rect(colour= "white",fill = "white")))
    

    在这里插入图片描述4、tsne降维
    #迭代次数 max_iter
    #困惑度perplexity
    #权衡速度与准确度,越小越精确,越大速度越快 theta

    library(Rtsne)
    #设置随机种子
    set.seed(1500)
    tsne_result <- Rtsne(data,
                         initial_dims = ncol(data),
                         pca = FALSE,
                         dims = 2,
                         check_duplicates = FALSE,
                         perplexity=30,
                         max_iter=2500,
                         theta=0.5)$Y
    colnames(tsne_result)=c('tsne_1','tsne_2')
    head(tsne_result)
    

    在这里插入图片描述

    tsne_data=data.frame(sample=rownames(data),
                         tsne_result,
                         group=c(rep('C',nrow(data)/2),rep('P',nrow(data)/2)))
    head(tsne_data)
    

    在这里插入图片描述

    library(ggplot2)
    ggplot(data=tsne_data,aes(x=tsne_1,y=tsne_2))+
      geom_point(aes(colour=group,shape=group),size=2)+
      theme(text=element_text(size=10),legend.title=element_blank(),
            panel.background = element_rect(fill = "white", colour = "black",size = 0.2), 
            legend.key = element_rect(fill = "white", colour = "white"),
            legend.background = (element_rect(colour= "white",fill = "white")))
    

    在这里插入图片描述

    展开全文
  • R语言数据降维——主成分分析

    千次阅读 多人点赞 2019-07-07 17:27:35
    R语言数据降维——主成分分析 一、项目环境 开发工具:RStudio R:3.5.2 相关包:sqldf,dplyr 二、导入数据 # 这里我们使用的是鸢尾花数据集(iris) data(iris) head(iris) Sepal.Length Sepal.Width ...

    R语言数据降维——主成分分析

    一、项目环境

    • 开发工具:RStudio
    • R:3.5.2
    • 相关包:sqldf,dplyr

    二、导入数据

    # 这里我们使用的是鸢尾花数据集(iris)
    data(iris)
    head(iris)
    
    Sepal.LengthSepal.WidthPetal.LengthPetal.WidthSpecies
    15.13.51.40.2setosa
    24.93.01.40.2setosa
    34.73.21.30.2setosa
    44.63.11.50.2setosa
    55.03.61.40.2setosa
    65.43.91.70.4setosa

    相关数据解释:

    • Sepal.Length:萼片长度

    • Sepal.Width:萼片宽度

    • Petal.Length:花瓣长度

    • Petal.Width:花瓣宽度

    • Species:鸢尾花品种

    三、 数据划分

    library(dplyr)
    library(sqldf)
    
    # 为数据集增加序号列(id)
    iris$id <- c(1:nrow(iris))
    
    # 将鸢尾花数据集中70%的数据划分为训练集
    iris_train <- sample_frac(iris, 0.7, replace = TRUE)
    
    # 使用sql语句将剩下的30%花费为测试集
    iris_test <- sqldf("
                   select *
                   from iris
                   where id not in (
                   select id
                   from iris_train
                   )
                   ")
    
    # 去除序号列(id)
    iris_train <- iris_train[,-6]
    iris_test <- iris_test[,-6]
    

    【注】:这里使用到sqldf包的函数sqldf函数来时间在R语言中使用SQL语句

    四、 进行主成分分析

    # 对鸢尾花数据集的前4列进行主成分分析
    iris_train_pca <- princomp(iris_train[,1:4])
    
    # 绘制碎石图
    screeplot(iris_train_pca, npcs = ncol(iris_train),type="lines")
    
    

    在这里插入图片描述

    【注】:碎石图的分析方法主要是根据纵坐标的值,纵坐标值越大就表示表示该主成分能够解释的方差的比例越大,因此一般情况下我们会选择纵坐标值较大的几个主成分

    五、 提取主成分的信息

    # 第一行是特征值(Standard deviation),
    # 第二列是方差的贡献率(Proportion of Variance)
    # 第三列是累计方差的贡献率(Cumulative Proportion)
    # 方差的贡献率: 标准化后的特征值,全部相加等于100%
    # 累计方差的贡献率:累加后的方差的贡献率
    summary(iris_train_pca)
    

    在这里插入图片描述

    【注】:特征值越大,它所对应的主成分变量包含的信息就越多

    ​ 从输出的结果我们可以得知,在输出的4个主成分中,前两个主成分就包含了原来4个指标98.20%的信息,也就是能够解释98.20%的方差。因此,将前两个作为鸢尾花数据集的主成分。

    六、 查看各个主成分与原始数据的关系

    # loadings代表每一个成分中之前特征系数
    iris_train_pca$loadings
    

    在这里插入图片描述
    ​ loadings显示的是载荷的内容,这个值实际上是主成分对于原始变量Sepal.Length,Sepal.Width ,Petal.Length ,Petal.Width的系数。也是特征值对应的特征向量,它们是线性无关的单位向量。第1列表示主成分一的得分系数,依次类推。据此可以写出由标准化变量所表达的主成分的关系式,即:

    C o m p . 1 = 0.366 × S e p a l . L e n g t h + 0.857 × P e t a l . L e n g t h + 0.351 × P e t a l . W i d t h Comp.1 = 0.366 × Sepal.Length + 0.857 × Petal.Length + 0.351 × Petal.Width Comp.1=0.366×Sepal.Length+0.857×Petal.Length+0.351×Petal.Width
    C o m p . 2 = 0.676 × S e p a l . L e n g t h + 0.709 × S e p a l . W i d t h − 0.194 × P e t a l . L e n g t h Comp.2 = 0.676 × Sepal.Length + 0.709 × Sepal.Width - 0.194 × Petal.Length Comp.2=0.676×Sepal.Length+0.709×Sepal.Width0.194×Petal.Length

    七、 对测试集进行降维处理

    # 使用之间建好的公式对测试集进行降维处理
    # new_test <- as.matrix(iris_test[,1:4])%*%as.matrix(iris_train_pca$loadings[,1:2])
    # 用上面这种方法进行计算可能会出现问题,建议使用下面这种
    new_test <- predict(iris_train_pca, iris_test[,-5])
    # 转化为数据框
    new_test <- as.data.frame(new_test)
    

    【注】:iris_train_pca$loadings[,1:2] 之所以这里取前两个数是因为之前主成分分析确定的。

    ​ 这样就达到了对数据进行降维的作用,同时可以将降维后的数据用与机器学习以降低维度过多,而造成的计算时间过长等问题。

    ​ 需要注意的是,这里的训练集和测试集维度必须完全相同,也就是说如果之前有与预先对训练集进行其他影响维度的操作,那么后续也需要对测试集进行相应的操作,才能保证降维的成功


    很长一段时间没有在csdn中写文章了,事实上后面自己学习过程中的大部分文档都是在语雀中完成的,基本都是自己写自己看。不过后面打算弄个公众号来更新和分享自己写的笔记_(:з」∠)_ 。如果感兴趣的话也可以关注一下。不感兴趣就算了(っ╥╯﹏╰╥c)。
    在这里插入图片描述
    然后写着个的另外一个原因就是,我可能会把一些文档搬运到自己的公众号,为了避免被说抄袭,所以说明一下。这样子_(:з」∠)_。

    展开全文
  • 流行学习,R语言模拟生成Swissroll,Helix, Twinpeaks,圆球等数据,通过pca,lle,isomap,tsne等方法对数据降维并可视化。
  • R语言学习_数据降维

    千次阅读 2019-04-04 21:38:18
    纬度灾难 变量过多(没用的变量) 变量相关(相关的变量) 解决办法 ...降维过程中不影响解的精度 消除多重共线性 数学工具 原变量线性组合得到新变量;方差的重新分配,保留几个方差最大的变量; ...

    纬度灾难
    变量过多(没用的变量)
    变量相关(相关的变量)
    解决办法
    剔除无用变量
    逐步回归
    向前引入法
    向后剔除法
    逐步筛选法
    Step函数
    AIC越小越好 AIC = n ln(SSE) + 2p
    主成分分析
    快速降维技术
    降维过程中不影响解的精度
    消除多重共线性

            数学工具
                原变量线性组合得到新变量;方差的重新分配,保留几个方差最大的变量;
                矩阵对角化
    
            R的函数
                princomp    这个函数是R中的标准PCA函数,可用cor,也可用cov协方差阵来做PCA
                predict
                loadings
                screeplot
        因子分析
            因子分析和主成分分析的区别
                主成分分析从“方差”出发
                因子分析从“相关性”出发
            因子分析的方法
                主成分法
                主因子法
                最大似然法
    
            因子分析的步骤
                观察相关系数矩阵
                提取因子变量
                因子变量命名
                计算因子得分(降维)
    
            R函数
                factanal    (极大似然法做因子分析)
                psych::principal (主成分法)、psync::fa  (主因子法)
                psych::factor.plot、psych::fa.diagram  (可视化)
    

    逐步回归法剔除无用变量代码示例:
    > mtcars
    mpg cyl disp hp drat wt qsec vs am gear carb
    Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4
    Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4
    Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1
    Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1
    Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2
    Valiant 18.1 6 225.0 105 2.76 3.460 20.22 1 0 3 1
    Duster 360 14.3 8 360.0 245 3.21 3.570 15.84 0 0 3 4
    Merc 240D 24.4 4 146.7 62 3.69 3.190 20.00 1 0 4 2
    Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
    Merc 280 19.2 6 167.6 123 3.92 3.440 18.30 1 0 4 4
    Merc 280C 17.8 6 167.6 123 3.92 3.440 18.90 1 0 4 4
    Merc 450SE 16.4 8 275.8 180 3.07 4.070 17.40 0 0 3 3
    Merc 450SL 17.3 8 275.8 180 3.07 3.730 17.60 0 0 3 3
    Merc 450SLC 15.2 8 275.8 180 3.07 3.780 18.00 0 0 3 3
    Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
    Lincoln Continental 10.4 8 460.0 215 3.00 5.424 17.82 0 0 3 4
    Chrysler Imperial 14.7 8 440.0 230 3.23 5.345 17.42 0 0 3 4
    Fiat 128 32.4 4 78.7 66 4.08 2.200 19.47 1 1 4 1
    Honda Civic 30.4 4 75.7 52 4.93 1.615 18.52 1 1 4 2
    Toyota Corolla 33.9 4 71.1 65 4.22 1.835 19.90 1 1 4 1
    Toyota Corona 21.5 4 120.1 97 3.70 2.465 20.01 1 0 3 1
    Dodge Challenger 15.5 8 318.0 150 2.76 3.520 16.87 0 0 3 2
    AMC Javelin 15.2 8 304.0 150 3.15 3.435 17.30 0 0 3 2
    Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
    Pontiac Firebird 19.2 8 400.0 175 3.08 3.845 17.05 0 0 3 2
    Fiat X1-9 27.3 4 79.0 66 4.08 1.935 18.90 1 1 4 1
    Porsche 914-2 26.0 4 120.3 91 4.43 2.140 16.70 0 1 5 2
    Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
    Ford Pantera L 15.8 8 351.0 264 4.22 3.170 14.50 0 1 5 4
    Ferrari Dino 19.7 6 145.0 175 3.62 2.770 15.50 0 1 5 6
    Maserati Bora 15.0 8 301.0 335 3.54 3.570 14.60 0 1 5 8
    Volvo 142E 21.4 4 121.0 109 4.11 2.780 18.60 1 1 4 2
    > carlm = lm(mpg~cyl+disp+hp+drat+wt+qsec,data = mtcars)
    > summary(mtcars)
    mpg cyl disp hp
    Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
    1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
    Median :19.20 Median :6.000 Median :196.3 Median :123.0
    Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
    3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
    Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
    drat wt qsec vs
    Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
    1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
    Median :3.695 Median :3.325 Median :17.71 Median :0.0000
    Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
    3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
    Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
    am gear carb
    Min. :0.0000 Min. :3.000 Min. :1.000
    1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
    Median :0.0000 Median :4.000 Median :2.000
    Mean :0.4062 Mean :3.688 Mean :2.812
    3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
    Max. :1.0000 Max. :5.000 Max. :8.000
    > carlm.step = step(carlm)
    Start: AIC=66.19
    mpg ~ cyl + disp + hp + drat + wt + qsec

           Df Sum of Sq    RSS    AIC
    - qsec  1     3.949 167.43 64.954
    - drat  1     5.209 168.69 65.194
    - cyl   1     6.652 170.13 65.466
    - disp  1     7.870 171.35 65.695
    - hp    1     8.744 172.22 65.857
    <none>              163.48 66.190
    - wt    1    72.580 236.06 75.947
    
    Step:  AIC=64.95
    mpg ~ cyl + disp + hp + drat + wt
    
           Df Sum of Sq    RSS    AIC
    - drat  1     3.018 170.44 63.526
    - disp  1     6.949 174.38 64.255
    <none>              167.43 64.954
    - cyl   1    15.411 182.84 65.772
    - hp    1    21.066 188.49 66.746
    - wt    1    77.476 244.90 75.124
    
    Step:  AIC=63.53
    mpg ~ cyl + disp + hp + wt
    
           Df Sum of Sq    RSS    AIC
    - disp  1     6.176 176.62 62.665
    <none>              170.44 63.526
    - hp    1    18.048 188.49 64.746
    - cyl   1    24.546 194.99 65.831
    - wt    1    90.925 261.37 75.206
    
    Step:  AIC=62.66
    mpg ~ cyl + hp + wt
    
           Df Sum of Sq    RSS    AIC
    <none>              176.62 62.665
    - hp    1    14.551 191.17 63.198
    - cyl   1    18.427 195.05 63.840
    - wt    1   115.354 291.98 76.750
    > summary(carlm.step)
    
    Call:
    lm(formula = mpg ~ cyl + hp + wt, data = mtcars)
    
    Residuals:
        Min      1Q  Median      3Q     Max
    -3.9290 -1.5598 -0.5311  1.1850  5.8986
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept) 38.75179    1.78686  21.687  < 2e-16 ***
    cyl         -0.94162    0.55092  -1.709 0.098480 .
    hp          -0.01804    0.01188  -1.519 0.140015
    wt          -3.16697    0.74058  -4.276 0.000199 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 2.512 on 28 degrees of freedom
    Multiple R-squared:  0.8431,	Adjusted R-squared:  0.8263
    F-statistic: 50.17 on 3 and 28 DF,  p-value: 2.184e-11
    
    > carlm2 = lm(mpg~cyl+wt,data = mtcars)
    > summary(carlm2)
    
    Call:
    lm(formula = mpg ~ cyl + wt, data = mtcars)
    
    Residuals:
        Min      1Q  Median      3Q     Max
    -4.2893 -1.5512 -0.4684  1.5743  6.1004
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)  39.6863     1.7150  23.141  < 2e-16 ***
    cyl          -1.5078     0.4147  -3.636 0.001064 **
    wt           -3.1910     0.7569  -4.216 0.000222 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 2.568 on 29 degrees of freedom
    Multiple R-squared:  0.8302,	Adjusted R-squared:  0.8185
    F-statistic: 70.91 on 2 and 29 DF,  p-value: 6.809e-12
    
    > carlm2 = lm(mpg~cyl*wt,data = mtcars)
    > summary(carlm2)
    
    Call:
    lm(formula = mpg ~ cyl * wt, data = mtcars)
    
    Residuals:
        Min      1Q  Median      3Q     Max
    -4.2288 -1.3495 -0.5042  1.4647  5.2344
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)  54.3068     6.1275   8.863 1.29e-09 ***
    cyl          -3.8032     1.0050  -3.784 0.000747 ***
    wt           -8.6556     2.3201  -3.731 0.000861 ***
    cyl:wt        0.8084     0.3273   2.470 0.019882 *
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 2.368 on 28 degrees of freedom
    Multiple R-squared:  0.8606,	Adjusted R-squared:  0.8457
    F-statistic: 57.62 on 3 and 28 DF,  p-value: 4.231e-12
    

    主成分分析R语言实例:
    # 第一步:将student.csv中的数据读入到程序中
    > # 30名学生的身高,体重,胸围,坐高
    > student = read.csv(‘e:/student.csv’,header = T)
    # 注:header = T表示将students_data.csv中的第一行数据设置为列名,这种情况下,
    student.csv中的第二行到最后一行数据作为data中的有效数据。header = F
    表示不将student.csv中的第一行数据设置为列名,这种情况下,student.csv
    中的第一行到最后一行数据作为data中的有效数据。
    > str(student)
    ‘data.frame’: 30 obs. of 4 variables:
    $ X1: int 148 139 160 149 159 142 153 150 151 139 …
    $ X2: int 41 34 49 36 45 31 43 43 42 31 …
    $ X3: int 72 71 77 67 80 66 76 77 77 68 …
    $ X4: int 78 76 86 79 86 76 83 79 80 74 …
    # 第二步:进行主成分分析
    > student.pr <- princomp(student, cor = TRUE)
    # 注:cor = T的意思是用相关系数进行主成分分析。
    # 第三步:观察主成分分析的详细情况
    > summary(student.pr)
    Importance of components:
    Comp.1 Comp.2 Comp.3 Comp.4
    Standard deviation 1.8817805 0.55980636 0.28179594 0.25711844
    Proportion of Variance 0.8852745 0.07834579 0.01985224 0.01652747
    Cumulative Proportion 0.8852745 0.96362029 0.98347253 1.00000000
    # 说明: 结果中的Comp.1、Comp.2、Comp.3和Comp.4是计算出来的主成分,Standard deviation代表每个主成分的标准差,
    Proportion of Variance代表每个主成分的贡献率,Cumulative Proportion代表各个主成分的累积贡献率。
    每个主成分都不属于X1、X2、X3和X4中的任何一个。第一主成分、第二主成分、第三主成分和第四主成分都是X1、X2、X3和X4的线性组合,
    也就是说最原始数据的成分经过线性变换得到了各个主成分。然而并不是每个主成分的作用都非常关键,因此,我们只选择作用比较关键的几个,
    一般地,选择累积贡献率达到八成的前几个主成分即可(这个实例中我们选择前两个,毕竟第二主成分的贡献率也比较大)。
    接下来,在得到主成分的基础上进行回归也好进行聚类也好,就不再使用原始的X1、X2、X3和X4了,而是使用主成分的数据。
    但现在还没有各个样本的主成分的数据,所以,最后一步就是得到各个样本主成分的数据。

    # 第四步:计算得到各个样本主成分的数据
    > predict(student.pr)
               Comp.1      Comp.2      Comp.3       Comp.4
     [1,] -0.06990950 -0.23813701  0.35509248 -0.266120139
     [2,] -1.59526340 -0.71847399 -0.32813232 -0.118056646
     [3,]  2.84793151  0.38956679  0.09731731 -0.279482487
     [4,] -0.75996988  0.80604335  0.04945722 -0.162949298
     [5,]  2.73966777  0.01718087 -0.36012615  0.358653044
     [6,] -2.10583168  0.32284393 -0.18600422 -0.036456084
     [7,]  1.42105591 -0.06053165 -0.21093321 -0.044223092
     [8,]  0.82583977 -0.78102576  0.27557798  0.057288572
     [9,]  0.93464402 -0.58469242  0.08814136  0.181037746
    [10,] -2.36463820 -0.36532199 -0.08840476  0.045520127
    [11,] -2.83741916  0.34875841 -0.03310423 -0.031146930
    [12,]  2.60851224  0.21278728  0.33398037  0.210157574
    [13,]  2.44253342 -0.16769496  0.46918095 -0.162987830
    [14,] -1.86630669  0.05021384 -0.37720280 -0.358821916
    [15,] -2.81347421 -0.31790107  0.03291329 -0.222035112
    [16,] -0.06392983  0.20718448 -0.04334340  0.703533624
    [17,]  1.55561022 -1.70439674  0.33126406  0.007551879
    [18,] -1.07392251 -0.06763418 -0.02283648  0.048606680
    [19,]  2.52174212  0.97274301 -0.12164633 -0.390667991
    [20,]  2.14072377  0.02217881 -0.37410972  0.129548960
    [21,]  0.79624422  0.16307887 -0.12781270 -0.294140762
    [22,] -0.28708321 -0.35744666  0.03962116  0.080991989
    [23,]  0.25151075  1.25555188  0.55617325  0.109068939
    [24,] -2.05706032  0.78894494  0.26552109  0.388088643
    [25,]  3.08596855 -0.05775318 -0.62110421 -0.218939612
    [26,]  0.16367555  0.04317932 -0.24481850  0.560248997
    [27,] -1.37265053  0.02220972  0.23378320 -0.257399715
    [28,] -2.16097778  0.13733233 -0.35589739  0.093123683
    [29,] -2.40434827 -0.48613137  0.16154441 -0.007914021
    [30,] -0.50287468  0.14734317  0.20590831 -0.122078819
    # 我们只保留Comp.1和Comp.2的数据即可。
    > screeplot(student.pr,type = 'lines')      #碎石图
    
    > consumedata = read.csv('e:/consume.csv')
    > lm1 = lm(Y~.,data = consumedata)
    > summary(lm1)
    
    Call:
    lm(formula = Y ~ ., data = consumedata)
    
    Residuals:
            1         2         3         4         5         6         7         8         9        10
     0.024803  0.079476  0.012381 -0.007025 -0.288345  0.216090 -0.142085  0.158360 -0.135964  0.082310
    
    Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
    (Intercept) -17.66768    5.94360  -2.973  0.03107 *
    X1            0.09006    0.02095   4.298  0.00773 **
    X2           -0.23132    0.07132  -3.243  0.02287 *
    X3            0.01806    0.03907   0.462  0.66328
    X4            0.42075    0.11847   3.552  0.01636 *
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.2037 on 5 degrees of freedom
    Multiple R-squared:  0.9988,	Adjusted R-squared:  0.9978
    F-statistic:  1021 on 4 and 5 DF,  p-value: 1.827e-07
    
    > # 观察自变量之间的相关性
    > library(corrgram)
    
    > corrgram(consumedata[,1:4],lower.panel = panel.conf,upper.panel = panel.pie,text.panel = panel.txt)
    > consumedata.pre = princomp(consumedata[,1:4],cor = TRUE)
    > summary(consumedata.pre)
    Importance of components:
                              Comp.1      Comp.2     Comp.3       Comp.4
    Standard deviation     1.9859037 0.199906992 0.11218966 0.0603085506
    Proportion of Variance 0.9859534 0.009990701 0.00314663 0.0009092803
    Cumulative Proportion  0.9859534 0.995944090 0.99909072 1.0000000000
    > screeplot(consumedata.pre,type = 'lines')
    
    
    > loadings(consumedata.pre)
    
    Loadings:
       Comp.1 Comp.2 Comp.3 Comp.4
    X1  0.502  0.237  0.579  0.598
    X2  0.500 -0.493 -0.610  0.367
    X3  0.498  0.707 -0.368 -0.342
    X4  0.501 -0.449  0.396 -0.626
    
                   Comp.1 Comp.2 Comp.3 Comp.4
    SS loadings      1.00   1.00   1.00   1.00
    Proportion Var   0.25   0.25   0.25   0.25
    Cumulative Var   0.25   0.50   0.75   1.00
    > # 做线性回归
    > consumedata$z1 = predict(consumedata.pre)[,1]
    > consumedata$z2 = predict(consumedata.pre)[,2]
    > consumedata.lm = lm(Y~z1+z2,data = consumedata)
    >
    > summary(consumedata.lm)
    
    Call:
    lm(formula = Y ~ z1 + z2, data = consumedata)
    
    Residuals:
         Min       1Q   Median       3Q      Max
    -0.74323 -0.29223  0.01746  0.30807  0.80849
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept) 14.03000    0.17125  81.927 1.06e-11 ***
    z1           2.06119    0.08623  23.903 5.70e-08 ***
    z2           0.62409    0.85665   0.729     0.49
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.5415 on 7 degrees of freedom
    Multiple R-squared:  0.9879,	Adjusted R-squared:  0.9845
    F-statistic: 285.9 on 2 and 7 DF,  p-value: 1.945e-07
    
    > consumedata.lm2 = lm(Y~z1,data = consumedata)
    > summary(consumedata.lm2)
    
    Call:
    lm(formula = Y ~ z1, data = consumedata)
    
    Residuals:
         Min       1Q   Median       3Q      Max
    -0.72237 -0.20946  0.05154  0.21032  0.81856
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept) 14.03000    0.16615   84.44 4.32e-13 ***
    z1           2.06119    0.08367   24.64 7.87e-09 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.5254 on 8 degrees of freedom
    Multiple R-squared:  0.987,	Adjusted R-squared:  0.9854
    F-statistic: 606.9 on 1 and 8 DF,  p-value: 7.873e-09
    
    
    
    > # install.packages('psych')
    > install.packages('psych')
    > library(psych)
    > R = read.csv('e:/sw.csv')
    > R = R[-1]
    > R
       身高  臂长 上肢长  腿长  体重  颈围  胸围  腰围
    1 1.000 0.846  0.805 0.859 0.473 0.398 0.301 0.382
    2 0.846 1.000  0.881 0.826 0.376 0.326 0.277 0.277
    3 0.805 0.881  1.000 0.801 0.380 0.319 0.237 0.345
    4 0.859 0.826  0.801 1.000 0.436 0.329 0.327 0.365
    5 0.473 0.376  0.380 0.436 1.000 0.762 0.730 0.629
    6 0.398 0.326  0.319 0.329 0.762 1.000 0.583 0.577
    7 0.301 0.277  0.237 0.327 0.730 0.583 1.000 0.539
    8 0.382 0.277  0.345 0.365 0.629 0.577 0.539 1.000
    > # 主成分法
    > pc = principal(r = R,nfactors = 2,rotate = 'varimax')
    > pc
    Principal Components Analysis
    Call: principal(r = R, nfactors = 2, rotate = "varimax")
    Standardized loadings (pattern matrix) based upon correlation matrix
       RC1  RC2   h2    u2 com
    1 0.90 0.27 0.88 0.120 1.2
    2 0.93 0.17 0.90 0.097 1.1
    3 0.92 0.18 0.87 0.129 1.1
    4 0.90 0.24 0.86 0.137 1.1
    5 0.25 0.89 0.85 0.151 1.2
    6 0.18 0.84 0.74 0.264 1.1
    7 0.11 0.84 0.71 0.289 1.0
    8 0.20 0.77 0.63 0.370 1.1
    
                           RC1  RC2
    SS loadings           3.46 2.98
    Proportion Var        0.43 0.37
    Cumulative Var        0.43 0.81
    Proportion Explained  0.54 0.46
    Cumulative Proportion 0.54 1.00
    
    Mean item complexity =  1.1
    Test of the hypothesis that 2 components are sufficient.
    
    The root mean square of the residuals (RMSR) is  0.05
    
    Fit based upon off diagonal values = 0.99
    > # 因子载荷图
    > par(mfrow = c(2,1))
    > factor.plot(pc)
    > # 因子结果图
    > fa.diagram(pc)
    > par(mfrow = c(1,1))
    > # 主因子法
    > fa = fa(r = R,nfactors = 2,rotate = 'varimax')
    > par(mfrow = c(2,1))
    > factor.plot(fa)
    
    > fa.diagram(fa,simple = T)
    > par(mfrow = c(1,1))
    
    展开全文
  • 探索性数据降维分析 本报告主要包含以下内容: 数据介绍 基本原理介绍 结合案例数据进行分析 最后总结 附上代码和参考 数据介绍 本报告所使用的是洛杉矶街区数据,其中包含每个街区的名字、收入中位数、公立学校...

    探索性数据降维分析

    本报告主要包含以下内容:

    1. 数据介绍
    2. 基本原理介绍
    3. 结合案例数据进行分析
    4. 最后总结
    5. 附上代码和参考

    数据介绍

    本报告所使用的是洛杉矶街区数据,其中包含每个街区的名字、收入中位数、公立学校API中位数、种族多样性、年龄中位数、有房家庭占比等14项字段,共有110个观测数据。本报告的主要目的是对这个数据的字段(变量)进行分析,并且探索性地尝试使用主成分分析和因子分析等降维方法来对数据进行降维分析。

    基本原理介绍

    主成分分析

    主成分分析是一种降维方法,通过原始数据一系列的线性变换找到对数组总体变异性(信息量)贡献比较大的主成分,这些线性变换保留了原始变量的大部分信息。其中,总体的变异性是指总体中所包含的信息,保留了原始变量的大部分信息是指保留了大部分的变异(PCA中使用的是方差)信息。

    主成分分析步骤:

    X 1 , … , X p X_1, \dots,X_p X1,,Xp p p p个随机变量。

    1. 将原始数据标准化,得到标准化数据矩阵
      KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲%开始数学环境 X= \lef…

    2. 计算变量X的相关系数矩阵: R = ( r i j ) p × p = X ′ X R = (r_{ij})_{p\times p}=X^{'}X R=(rij)p×p=XX

    3. R R R的特征值 λ 1 ≥ λ 2 ≥ ⋯ ≥ λ p > 0 \lambda_1 \ge \lambda_2 \ge \cdots \ge\lambda_p \gt 0 λ1λ2λp>0及相应的单位特征向量。
      KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ u_1={ \left[ \…

    4. 将特征值排序,计算方差贡献率和累计方差贡献率,找到达到累计贡献率的前几个特征值 λ 1 , λ 2 , ⋯   , λ q \lambda_1, \lambda_2, \cdots, \lambda_q λ1,λ2,,λq所对应的特征向量 u 1 , u 2 , ⋯   , u q u_1, u_2, \cdots,u_q u1,u2,,uq,不妨设我们所需要的累积贡献率为80%,即 ∑ i = 1 q λ i ≥ 0.8 \sum\limits_{i=1}^{q}\lambda_i \ge 0.8 i=1qλi0.8

    5. 确定主成分,用原指标的线性组合来计算各主成分得分:以各主成分对原指标的相关系数为权,将各主成分表示远原指标的线性组合,即:
      C j = u j 1 x 1 + u j 2 x 2 + ⋯ + u j p x p , j = 1 , 2 , ⋯   , 1 C_j = u_{j1}x_1 + u_{j2}x_2 + \cdots + u_{jp}x_p,\quad j=1, 2, \cdots, 1 Cj=uj1x1+uj2x2++ujpxp,j=1,2,,1

    主成分分析的主要理论结果:

    1. 对任意 1 ≤ i ≤ p 1\le i \le p 1ip,第 i i i各主成分的系数向量 a i a_i ai取为 u i u_i ui,因此第 i i i个主成分是 Y i = u i 1 x 1 + u i 2 x 2 + ⋯ + u i p x p Y_i = u_{i1}x_1 + u_{i2}x_2 +\cdots + u_{ip}x_p Yi=ui1x1+ui2x2++uipxp.
    2. 对任意 1 ≤ i ≤ p , V a r ( Y i ) = λ i 1\le i \le p, Var(Y_i) = \lambda_i 1ip,Var(Yi)=λi.
    3. 对任意 1 ≤ i ≠ j ≤ p , C o v ( Y i , Y j ) = 0 1 \le i \ne j \le p, Cov(Y_i, Y_j)=0 1i=jp,Cov(Yi,Yj)=0,即 Y i Y_i Yi Y j Y_j Yj不线性相关。
    4. 数据经过以上的标准化之后,总方差为 ∑ i = 0 p λ i = p \sum\limits_{i=0}^p \lambda_i=p i=0pλi=p,,第 i i i个主成解释的方差比例为 λ i / p \lambda_i/p λi/p,前 q q q个主成分解释总方差的比例为 ∑ i = 1 q λ i / p \sum_{i=1}^q \lambda_i / p i=1qλi/p.
    5. 对任意 1 ≤ i , k ≤ p 1 \le i, k \le p 1i,kp, Y i Y_i Yi X k X_k Xk的相关系数为 C o r r ( Y i , X k ) = e i k λ i Corr(Y_i, X_k)=e_{ik}\sqrt{\lambda_i} Corr(Yi,Xk)=eikλi .

    主成分个数的选择

    1. Kaiser准则:保留那些对应特征值大于所有特征值的平均之的主成分,解释总方差比例大于平均解释比例的主成分,在这里数据标准化周就保留对应的特征值大于1的主成分。
    2. 总方差中被前 q q q个主成分解释的比例达到一定大小**(常用)**。
    3. 保留的主成分在实际应用中具有可解释性。
    4. 使用崖底碎石图绘出特征值与其顺序的关系,找到一个拐点,使得此节点后对应的特征向量都比较小,然后选择拐点之前的一点。

    主成分的含义:

    1. 对第 i i i个主成分,选择主成分中原变量系数绝对值比较大的变量,作为主成分的主要解释。
    2. 根据相关系数来解释,计算第 i i i个主成分与各个原始变量之间的相关系数,根据相关系数比较大的来作为主成分的解释。

    探索性因子分析

    不同主成分分析直接将将一些对总体变异贡献率比较高的特征值所对应的特征向量,来对原始变量进行线性组合,得到新的主成分。因子分析主要是将变量的变异性解释为由两种类型的因子的变异,分别是潜在的、不可观测的公共因子和只与该变量有关的特殊因子。这里因子分析主要是寻找少量公共因子(对应于主成分分析的主成分)来解释一组输入变量的共同的变异性。

    一些符号:

    假设 X = ( x 1 , x 2 , ⋯   , x n ) X = (x_1, x_2, \cdots, x_n) X=(x1,x2,,xn)是一个 p p p维随机向量,它的均值向量为 μ = ( μ 1 , μ 2 , ⋯   , μ p ) T \mu = (\mu_1, \mu_2, \cdots, \mu_p)^T μ=(μ1,μ2,,μp)T,协方差矩阵为 Σ \Sigma Σ,其对角线上的值 σ k 2 \sigma_k^2 σk2表示 X k X_k Xk的方差 ( k = 1 , ⋯   , p ) (k=1, \cdots, p) (k=1,,p),令 F 1 , ⋯   , F q ( q ≤ p ) F_1, \cdots, F_q (q \le p) F1,,Fq(qp)表示q个潜在的公共因子,令 ϵ 1 , ⋯   , ϵ q \epsilon_1, \cdots, \epsilon_q ϵ1,,ϵq表示特殊因子。正交矩阵模型:
    KaTeX parse error: No such environment: split at position 8: \begin{̲s̲p̲l̲i̲t̲}̲ &X_1 - \mu_1 =…
    写成矩阵形式 X − μ = L F + ϵ X - \mu = LF + \epsilon Xμ=LF+ϵ, F = ( F 1 , ⋯   , F q ) F=(F_1, \cdots, F_q) F=(F1,,Fq)是共因子, ϵ = ( ϵ 1 , ⋯   , ϵ q ) \epsilon = (\epsilon_1, \cdots, \epsilon_q) ϵ=(ϵ1,,ϵq)是特殊因子, L L L是载荷矩阵,其中第 k k k行第 i i i列的值 l k i , k = 1 , ⋯   , p , i = 1 , ⋯   , q l_{ki}, k=1, \cdots, p, i=1, \cdots, q lki,k=1,,p,i=1,,q表示 X k X_k Xk在因子 F i F_i Fi上的载荷。

    由于公共因子和特殊因子是不可观测的,所以需要作出一些假定

    1. E ( F i ) = 0 , V a r ( F i ) = 1 , ∀ 1 ≤ i ≤ q E(F_i)=0, Var(F_i)=1, \forall 1 \le i \le q E(Fi)=0,Var(Fi)=1,1iq.

      C o v ( X i , X k ) = 0 , ∀ i , k = 1 , ⋯   , q Cov(X_i, X_k)=0, \forall i, k = 1, \cdots, q Cov(Xi,Xk)=0,i,k=1,,q

    2. E ( ϵ k ) = 0 , V a r ( ϵ k ) = Ψ k E(\epsilon_k)=0, Var(\epsilon_k)=\Psi_k E(ϵk)=0,Var(ϵk)=Ψk

      C o v ( ϵ k , ϵ m ) = 0 , ∀ k , m = 1 , ⋯   , q Cov(\epsilon_k, \epsilon_m)=0, \forall k, m = 1, \cdots, q Cov(ϵk,ϵm)=0,k,m=1,,q

    3. C o v ( F i , ϵ k ) = 0 , ∀ 1 ≤ i ≤ q , 1 ≤ k ≤ p Cov(F_i, \epsilon_k)=0,\forall 1\le i \le q, 1 \le k \le p Cov(Fi,ϵk)=0,1iq,1kp

    主要结论:

    1. V a r ( X k ) = l k 1 2 + l k 2 2 + ⋯ + l k q 2 + Ψ k Var(X_k)=l_{k1}^2 + l_{k2}^2 + \cdots + l_{kq}^2 + \Psi_k Var(Xk)=lk12+lk22++lkq2+Ψk

      l k 1 2 + l k 2 2 + ⋯ + l k q 2 l_{k1}^2 + l_{k2}^2 + \cdots + l_{kq}^2 lk12+lk22++lkq2平方和是方差中能被公共因子解释的部分,称为共性方差。

      Ψ k \Psi_k Ψk是不能被公共因子解释的部分,成为 X k X_k Xk的特殊方差。

    2. C o v ( X k , X m ) = l k 1 l m 1 + l k 2 l m 2 + ⋯ + l k q l m q , ∀ 1 ≤ k ≠ m ≤ p Cov(X_k, X_m)=l_{k1}l_{m1} + l_{k2}l_{m2} + \cdots +l_{kq}l_{mq}, \forall 1 \le k \ne m \le p Cov(Xk,Xm)=lk1lm1+lk2lm2++lkqlmq,1k=mp

    3. C o v ( X k , F i ) = l k i , ∀ k = 1 , ⋯   , p , i = 1 , ⋯   , q Cov(X_k, F_i)= l_ki, \forall k=1, \cdots, p, i =1, \cdots, q Cov(Xk,Fi)=lki,k=1,,p,i=1,,q

    公共因子的解释,同主成分分析,对 F i F_i Fi载荷系数的绝对值较大的输入的变量来解释。

    载荷矩阵估计:

    主要有主成分法(主因子法)和最大似然估计法。

    主成分法将协方差 Σ \Sigma Σ(标准化后为相关系数矩阵 R R R)拆分为 Σ = λ 1 e 1 e 1 T + λ 2 e 2 e 2 T + ⋯ + λ p e p e p T \Sigma=\lambda_1 e_1 e_1^T+\lambda_2 e_2 e_2^T+ \cdots +\lambda_p e_p e_p^T Σ=λ1e1e1T+λ2e2e2T++λpepepT,将式子中前 q q q项归结为主因子解释的部分,后面 p − q p-q pq项归结为特殊因子,得到 L L L Ψ \Psi Ψ的估计 L ~ \tilde{L} L~ Ψ ~ \tilde{\Psi} Ψ~

    最大似然估计法:

    假定公共因子 F F F和特殊因子 ϵ \epsilon ϵ服从正态分布,可以得因子载荷的最大似然估计。

    得到载荷矩阵之后,如果得到的因子对原始变量的可解释性不强,则还需要进行因子旋转,以得到较为明显的实际含义。

    旋转后的载荷矩阵满足以下几个条件:

    (1)对于任意因子 F i , i = 1 , ⋯   , q F_i, i=1, \cdots, q Fi,i=1,,q而言,只有少数输入变量在该因子上的绝对值较大,其他都接近于0.

    (2)对任意输入变量 X 1 , ⋯   , X p X_1,\cdots, X_p X1,,Xp而言,它只在少数因子上的载荷 l k i , k = 1 , ⋯   , p , i = 1 , ⋯   , q l_{ki},k=1, \cdots, p, i=1, \cdots, q lki,k=1,,p,i=1,,q的绝对值较大,在其他因子上的载荷接近于0.

    (3)任何两个因子对应的载荷呈现不同的模式,因为在解释时这两个因子具有不同的含义。

    ​ 比如对于因子 F 1 , F 2 F_1,F_2 F1,F2, 对于第一个变量 X 1 X_1 X1,因子 F 1 F_1 F1的载荷是 l 11 l_{11} l11,因子 F 2 F_2 F2的载荷是 l 12 l_{12} l12
    KaTeX parse error: No such environment: split at position 8: \begin{̲s̲p̲l̲i̲t̲}̲ &X_1 - \mu_1 =…

    多维标度分析

    **介绍:**它是一个在低维空间展现高维数据的可视化方法,使得低维空间中观测点之间的距离与高维空间中观测点之间的距离大致匹配。

    多维度标度根据度量的形式分为度量非度量,度量形式直接采用观测点之间的距离,非度量形式采用观测点之间距离的排序。这里介绍非度量形式。

    假设一共有N个观测,他们之间两两匹配的对数为 M = N ( N − 1 ) / 2 M=N(N-1)/2 M=N(N1)/2,计算高维空间中M对观测数据的距离,对其近似排序可以得到
    d i 1 k 1 ∗ < d i 2 k 2 ∗ < ⋯ < d i M k M ∗ d_{i_1k_1}^* \lt d_{i_2k_2}^* \lt \cdots < d_{i_Mk_M}^* di1k1<di2k2<<diMkM
    设低维空间的维度为q,将N个观测点放置在q维空间中,即每个观测点用一个q维坐标向量代表。令 d i j k j ( q ) d_{i_j k_j}^{(q)} dijkj(q)来表示q维空间中观测点之间的距离。

    给出一些定义:

    应力函数:
    S t r e s s ( q ) = { ∑ j = 1 M ( d i j k j ( q ) − d ^ i j k j ( q ) ) 2 ∑ j = 1 M [ d i j k j ( q ) ] 2 } 1 / 2 Stress(q) = \{\cfrac{\sum_{j=1}^M(d_{i_jk_j}^{(q)} - \hat{d}_{i_jk_j}^{(q)})^2}{\sum_{j=1}^M[d_{i_jk_j}^{(q)}]^2}\}^{1/2} Stress(q)={j=1M[dijkj(q)]2j=1M(dijkj(q)d^ijkj(q))2}1/2
    或者:
    S S t r e s s ( q ) = { ∑ j = 1 M [ ( d i j k j ( q ) ) 2 − ( d ^ i j k j ( q ) ) 2 ] 2 ∑ j = 1 M [ d i j k j ( q ) ] 4 } 1 / 2 SStress(q) = \{\cfrac{\sum_{j=1}^M[(d_{i_jk_j}^{(q)})^2 - (\hat{d}_{i_jk_j}^{(q)})^2] ^2} {\sum_{j=1}^M[d_{i_jk_j}^{(q)}]^4}\}^{1/2} SStress(q)={j=1M[dijkj(q)]4j=1M[(dijkj(q))2(d^ijkj(q))2]2}1/2
    应力函数的值越小,低维空间中观测点之间距离的排序与原来高维空间中观测点之间的距离的排序的一致性越高。

    算法步骤:

    1. 初始化N个观测点在q维空间的坐标向量,并据此计算 d i j k j ( q ) , j = 1 , ⋯   , M d_{i_j k_j}^{(q)}, j=1, \cdots, M dijkj(q),j=1,,M

    2. 在每次循环中:

      (1)固定 d i j k j ( q ) , j = 1 , ⋯   , M d_{i_j k_j}^{(q)}, j=1, \cdots, M dijkj(q),j=1,,M,寻找M个数值 d ^ i j k j ( q ) , j = 1 , ⋯   , M \hat{d}_{i_j k_j}^{(q)}, j=1, \cdots, M d^ijkj(q),j=1,,M,使得他们的排序与 d i j k j ∗ , j = 1 , ⋯   , M d_{i_j k_j}^{*}, j=1, \cdots, M dijkj,j=1,,M的排序完全一致,并且应力函数的值达到最小;

      (2)固定 d ^ i j k j ( q ) , j = 1 , ⋯   , M \hat{d}_{i_j k_j}^{(q)}, j=1, \cdots, M d^ijkj(q),j=1,,M,寻找N个观测点在q维空间的坐标向量,使得应力函数的值达到最小。

    ​ 持续循环直到应力函数的值无法减小为止。

    案例分析

    尝试对数据进行主成分降维

    在进行主成分分析之前,先对数据标准化,这里R语言中的princomp函数就在将原始数据输入进去之后就已经是标准化数据了,所以我们不需要再对数据进行标准化,直接将原始数据输入到princomp函数中去即可。

    在进行主成分分析之后,查看数据的因子载荷,找到对总体变异性贡献比较大的q个主成分。
    在这里插入图片描述

    如上图所示,第一个主成分对总体数据变异的贡献率为0.370107,第二个主成分对总体贡献率为0.1507379,此时的累计贡献率为0.5208449.

    这里我们选择使用Kaiser准则来找到前q个主成分:

    首先画图:
    在这里插入图片描述

    由图可知,第五个主成分是碎石图的一个拐点,在此点之前,特征值都比较大,在此点之后,特征值都比较小,所以根据Kaiser准则,我们选择前4个主成分来作为代表原始变量,进而达到降维的目的。
    在这里插入图片描述

    再对这里的主成分的含义进行解释,由上图所示,对于第一个主成分,它代表是是总体情况,第一个主成分主要由Income, Age, White三个变量来解释;第二个主成分主要有Black, Population, Latitude来解释;第三个主成分主要由Diversit, Asian, Area来解释。依次类推可以对剩下的一个主成分进行解释。

    尝试使用双标图来将各个观测和各个变量绘制的同一张图上。得到
    在这里插入图片描述

    由图可知,第19,16,7观测为异常观测。图中第一个主成分在Latitude和Area等变量上的系数为正,在Latino和Black等变量上的系数为负。

    **最终降维结果:**取成分的个数为4,最后将原始14个变量降维至4个变量。

    尝试使用因子分析来降维

    再接着尝试使用因子分析来降维,再R语言中的factanal函数中的代表因子个数的参数factors设置为8,得到因子分析的结果:
    在这里插入图片描述

    这暂时还不能够直观地看出在探索性因子分析中,选的多少个因子才是合适的,因此我们根据以上因子分析的因子方差来画图:
    在这里插入图片描述

    由上图可知,因子个数取3是一个拐点,结合上面的累计误差知道当取因子个数为2时,累计误差仅为0.373,这两个因子并不能够很好地代表整个总体。所以我们结合Kaiser准则和总方差中的累计贡献率,假设我们要求在80%左右就可以了,这里我们选择因子个数为7,这时7个因子在总方差中的累计贡献率近似达到80%。而且上图也表明,当因子个数超过7时,累计贡献率会从1.081急剧下降到0.188.

    最终选择因子个数为7,得到的载荷矩阵为:
    在这里插入图片描述

    由图可知,第一个因子可以很好地由Lantino, White来解释;第二个因子可以由第Vets来解释;第三个因子可以由Black解释。此外,黑人和白人是不相关的,这符合因子分析假设,各个因子之间是不相关的假设。

    **最终降维结果:**取因子的个数为7,最后将原始14个变量降维至7个变量。

    尝试使用多维标度分析来降维查看数据

    分析的流程:首先,将原始110行,14列的数据转置,得到14列,110行的数据。再求出利用曼哈顿的度量方法来找到数据之间的距离。将数据转换成矩阵形式,使用非度量的方式来分析原始数据,得到应力值为14.12094. 这说明使用2维空间对原来高维空间拟合的效果一般。

    当低维维度设置为2时,进行排序之后的数据进行排序得到:
    在这里插入图片描述

    由图所示,在低维空间中,Income收入和White白种人距离比较接近,Latitude和Area比较接近,Longitude以及Latino距离其他变量都比较远。

    因为当将低维维度设置为2时,效果不好,所以这里我们再尝试使用其他的较低的维度来对原来高维空间进行拟合,发现将低维的维度设置为7应力值才能够降到1以下。效果都不太好,所以最后不考虑使用多维标度分析来进行降维

    总结

    本报告首先介绍了原始数据,再比较详细地分别介绍了主成分分析,因子分析以及多标度分析的原理,并且还结合案例数据来探索进行数据降维。

    在因子分析选择主要因子个数时,仅凭借个人直觉判定应该首先使用碎石图来进行初步地判断,再结合因子对总方差的累计贡献率,最终选择因子个数为7,感觉有些多理论性又不强,是否正确我觉得在后续的学习中还需要进一步地学习去验证。

    模型最后多标度分析并只尝试了使用使用二维空间的数据来对高维空间进行拟合,这里得到的应力函数的值维14.12094,使用课本上的评价方法来说,是很一般的拟合,在后续的探索中我发现,随着维度的增大,对原始数据的拟合的应力值也是相应增大的,但是考虑到篇幅的问题就没有添加上去,后续要需要继续探索学习。

    参考

    [1]高惠璇.应用多元统计分析[M] 北京:北京大学出版社.2008.

    [2]王斌会.多元统计分析及R语言建模[M]广州:暨南大学出版社.2014.

    [3] 如何理解“方差越大信息量就越多”?

    [4] 机器学习实战之PCA

    [5] 一篇深入剖析PCA的好文

    [6] PCA:详细解释主成分分析

    [7] R语言 | 典型相关分析、多维标度法以及综合评价方法及R使用

    [8] R语言 3.14 多维标度法MDS

    代码

    ##加载程序包
    library(dplyr)
    library(ggplot2)
    
    # ---------读入数据---------
    setwd("D:/lagua/CODING/R-learn/R-code/Chap5_DimensionReduction")
    street <- read.csv("LANeighborhoods.csv", header=T, skip=1)
    colnames(street)
    # 去除掉第一列,因为后面需要计算相关系数矩阵
    street = street[, -1]
    summary(street)
    
    
    # ---------主成分分析---------
    help("princomp")
    streetout <- princomp(street,cor = T,scores = T)
    summary(streetout,loadings=T)
    
    streetout <- princomp(scale(street),cor = T,scores = T)
    #streetout$scores记录了每个观测的主成分得分。
    streetout$scores
    ##显示分析结果
    summary(streetout,loadings=T)
    
    ##画崖底碎石图
    screeplot(streetout,type = "lines")
    
    ##画前两个主成分的双标图
    biplot(streetout,choices = 1:2,col="black")
    help(biplot)
    
    
    # ---------因子分析---------
    # help("factanal")
    # 使用极大似然估计,
    # 因子旋转:方差最大
    streetout <- factanal(scale(street),fm="ml",factors = 8,rotation = "varimax")
    streetout # 查看结果
    # 画碎石图,查看方差变化情况,以便选择因子个数
    plot(c(3.037,   2.183,   1.311,   1.265,   1.099,   1.082,   1.081,   0.188),
         type="o",
         xlab="factor",
         ylab="Variance")
    
    # 选择因子个数为7
    streetout <- factanal(scale(street),fm="ml",factors = 4,rotation = "varimax")
    streetout
    plot(c(3.037,   2.183,   1.311,   1.265),
         type="o",
         xlab="factor",
         ylab="Variance")
    
    
    # 显示分析结果
    streetout
    #streetout数据集包含两个公共因子的载荷矩阵Loadings,
    summary(streetout)
    
    
    # ---------多维标度分析---------
    tmpstreet <- t(scale(street))
    help(dist)
    # 求出距离矩阵, 使用L_1范数--曼哈顿距离
    diststreet <- dist(tmpstreet,method = "manhattan") 
    
    
    # ---------度量形式
    # help("cmdscale")
    out <- cmdscale(diststreet) %>% as.data.frame()
    #使用cmdscale函数进行多维标度分析。
    
    ggplot(out,aes(x=V1,y=V2))+
      geom_point()+
      geom_text(aes(y=V2+5,label=row.names(out)))
    
    # ---------非度量形式
    library(MASS)
    D = as.matrix(diststreet) # 转化为矩阵形式
    help("isoMDS")
    MDS = isoMDS(D, k=2) # 非度量形式的多维标度分析
    # MDS
    MDS$stress # 查看应力值
    x = MDS$points[, 1]
    y = MDS$points[, 2]
    # 对代表高维空间数据的低维空间数据画图
    ggplot(out,aes(x=x,y=y))+
      geom_point()+
      geom_text(aes(y=y+5,label=row.names(D)))
    
    
    # ----------查看其他更低维度的应力值
    for (i in 3:9){
      MDS = isoMDS(D, k=i) # 非度量形式的多维标度分析
      # MDS
      print(MDS$stress) # 查看应力值
    }
    
    as.matrix(diststreet) # 转化为矩阵形式
    help("isoMDS")
    MDS = isoMDS(D, k=2) # 非度量形式的多维标度分析
    # MDS
    MDS$stress # 查看应力值
    x = MDS$points[, 1]
    y = MDS$points[, 2]
    # 对代表高维空间数据的低维空间数据画图
    ggplot(out,aes(x=x,y=y))+
      geom_point()+
      geom_text(aes(y=y+5,label=row.names(D)))
    
    
    # ----------查看其他更低维度的应力值
    for (i in 3:9){
      MDS = isoMDS(D, k=i) # 非度量形式的多维标度分析
      # MDS
      print(MDS$stress) # 查看应力值
    }
    
    展开全文
  • 本节书摘来自华章出版社《R语言数据挖掘》一书中的第1章,第1.13节,作者[哈萨克斯坦]贝特·麦克哈贝尔(Bater Makhabel),李洪成 许金炜 段力辉 译,更多章节内容可以访问云栖社区“华章计算机”公众号查看。...
  • R语言变量降维分析

    千次阅读 2019-02-13 09:21:25
    变量降维: (Variable dimension reduction) 涉及因子分析/主成分分析等,通过使用这个工具,可以将多个变量减少,用新的核心变量进行替代,并将新变量用线性关系表示。从而减少变量字段过多造成的数据分析复杂度。...
  • *降维的意义:克服维数灾难,获取本质特征,节省存储空间,去除无用噪声,实现数据可视化 强烈推荐几篇博客: https://www.douban.com/note/469279998/ http://bindog.github.io/blog/2016/06/04/f
  • 数据降维和参数降维的简单理解

    千次阅读 2017-03-22 10:46:45
    数据降维
  • R语言实现PCA降维

    2021-02-20 09:52:52
    R代码实现PCA降维过程 链接: http://blog.sciencenet.cn/blog-3448646-1270918.html 数据集是经典的鸢尾花 做简单的记录,只附代码和结果, 下面展示一些 内联代码片。 具体自行操作,可交流 install.packages(...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 7,036
精华内容 2,814
关键字:

r语言数据降维