數(shù)據(jù)建模與預(yù)測(cè)期末復(fù)習(xí)

萬瑾琳

信息管理與電子商務(wù)系

辦公室：7＃305

R語言基礎(chǔ)

安裝包

install.packages("name")

加載包

Library(name)

變量賦值

a=1; a<-1

向量:

賦值

a=c(1,2); a<-c("xxxx","yyyy")

min(x), max(x), range(x)分別表示向量x的最小分量胶坠、最大分量和向量x的范圍

與min(x), max(x)有關(guān)的函數(shù)是which. min(), which. max(),表示在第幾個(gè)分量求到其值

等差數(shù)列

x<-1:30

等差數(shù)列運(yùn)算大于加減乘除運(yùn)算

等間隔數(shù)列

seq(from,to,by = 間隔)

檢測(cè)是否為空值并替換

s<-is.na(data)

data[is.na(data)]<-num

字符型向量

分隔用的字符可用sep參數(shù)指定趟紊，如：

paste("result.", 1:4, sep="")

[1] “result.1” “result.2” “result.3” “result.4”

image-20211213192203407

多維數(shù)組和矩陣

生成數(shù)組矩陣

dim(a)<-c(3,4)
column.names <- c("col1","col2","col3") #定義列名
row.names <- c("row1","row2","row3") #定義行名
dimnames = list(column.names,row.names)
#矩陣行數(shù)
nrow(X)
#矩陣列數(shù)
ncol(X)

生成數(shù)組或矩陣－用matrix()函數(shù)構(gòu)造矩陣

A<- matrix(1:15,nrow=3,ncol=5,byrow=TRUE)

矩陣轉(zhuǎn)置

T()函數(shù)

如果矩陣A和B具有相同的維數(shù)癣猾，則A＊B表示矩陣中對(duì)應(yīng)的元素的乘積彼绷，A％＊％B表示通常意義下的兩個(gè)矩陣的成積（要求矩陣A的列數(shù)等于矩陣B的行數(shù)

生成對(duì)角陣和矩陣取對(duì)角陣運(yùn)算

函數(shù)diag()依賴于它的變量：

當(dāng)v是一個(gè)向量時(shí)，diag(v)表示以v的元素為對(duì)角線元素的對(duì)角陣.當(dāng)M是一個(gè)矩陣時(shí)斋配，則diag(M)表示的是取M對(duì)角線上的元素的向量.

解線性方程組和求矩陣的逆矩陣

#解方程組Ax＝b的解x
A<-t(array(c(1:8,10),dim=c(3,3)))
b<-c(1,1,1)
x<-solve(A,b)

#求逆矩陣
b<-solve(A)

求矩陣的特征值與特征向量

#求矩陣特征根與特征向量
ev<-eigen(A)
ev

ev存放著對(duì)稱矩陣Sm特征值和特征向量孔飒，是由列表形式給出的，其中evvectors是Sm的特征向量構(gòu)成的矩陣.

image-20211213194050021

矩陣合并

#求矩陣的合并 函數(shù)cbind()把其向量縱向合并
x1<-c(1,2,3,4,5,6,7)
x2<-c(2,3,4,6,7,9,10)
cbind(x1,x2)

#函數(shù)rbind()把其向量橫向合并
x1<-c(1,2,3,4,5,6,7)
x2<-c(2,3,4,6,7,9,10)
rbind(x1,x2)

image-20211213194445173

數(shù)據(jù)框

第二章 P63

構(gòu)建數(shù)據(jù)框

data.frame()

X<-data.frame(
name=c("amy","james","jeffrey","john"),
sex=c("F","M","F","M"), 
age=c(13,12,13,12),
height=c(56.5,57.3,62.5,59.0),
weight=c(84.0,83.0,84.8,99.5))

image-20211213194945928

數(shù)據(jù)的調(diào)用?

從剪切板調(diào)用

#從剪貼板讀取數(shù)據(jù)
#打開excel數(shù)據(jù)坏瞄，用ctrl+C復(fù)制數(shù)據(jù)
#然后調(diào)用read.delim函數(shù)，從剪貼板讀取數(shù)據(jù)

b <-read.delim("clipboard")

從文本調(diào)用

read.table(file, header = FALSE, sep = "")

常用參數(shù)說明：

file 文件路徑

header 第一行是否為列名稱

sep 字段之間的分隔符

quote字符串的標(biāo)記甩卓，默認(rèn)為"

dec小數(shù)點(diǎn)的符號(hào)鸠匀，默認(rèn)為.

讀取excel/csv

X<-read.xlsx("F:\\OneDrive - chenxiangxi\\桌面\\2.xlsx",sheetIndex = 1) #安裝xlsx
X<-read.csv("textdata.csv")

儲(chǔ)存文件

#存儲(chǔ)到csv
write.csv(mydata,"G:/xxx/Inventory.csv")

#存儲(chǔ)數(shù)據(jù)到txt文件,制表符分隔
write.table(mydata,""G:/xxx/Inventory.txt", sep="\t")

#存儲(chǔ)數(shù)據(jù)到xlsx文件
write.xlsx(mydata, "G:/xxx/Inventor..xlsx", sheetName="Sheet1")

多元數(shù)據(jù)直觀表示

均值

mean(x,trim=0,na.rm=FALSE)

x是對(duì)象，trim是在計(jì)算均值前去掉x兩端觀察值的比例猛频，默認(rèn)值為0狮崩，即包括全部數(shù)據(jù)。當(dāng)na. rm＝TRUE時(shí)鹿寻，允許數(shù)據(jù)中有缺失數(shù)據(jù)。

求矩陣平均值

用apply()函數(shù)

apply(x, 1, mean)  #計(jì)算矩陣各行的均值
apply(x, 2, mean)  #計(jì)算矩陣各列的均值

方差

var(x)計(jì)算樣本方差

sd(x)計(jì)算樣本標(biāo)準(zhǔn)差

變異系數(shù)CV

當(dāng)CV>15%時(shí)诽凌，說明數(shù)據(jù)離散程度過高毡熏，需要注意

cv<-100*sd(x)/mean(x)

偏度

偏度系數(shù)是刻畫數(shù)據(jù)的對(duì)稱性指標(biāo)。關(guān)于均值對(duì)稱的數(shù)據(jù)其偏度系數(shù)為0侣诵，越趨近于0痢法，數(shù)據(jù)越服從正態(tài)分布

右側(cè)更分散的數(shù)據(jù)偏度系數(shù)為正，左側(cè)更分散的數(shù)據(jù)偏度系數(shù)為負(fù)

n<-length(x)
m<-mean(x)
s<-sd(x)
g1<-n/((n-1)*(n-2))*sum((x-m)^3)/s^4

峰度

當(dāng)數(shù)據(jù)的總體分布為正態(tài)分布時(shí)杜顺，峰度系數(shù)近似為0

當(dāng)分布較正態(tài)分布的尾部更分散時(shí)财搁，峰度系數(shù)為正；否則為負(fù)

當(dāng)峰度系數(shù)為正時(shí)躬络，兩側(cè)數(shù)據(jù)較多尖奔，當(dāng)峰度系數(shù)為負(fù)，兩端數(shù)據(jù)較少

g2<-((n*(n+1))/(n-1)*(n-2)*(n-3))*sum((x-m)^4/s^4-(3*(n-1)^2)/((n-2)*(n-3)))

描述統(tǒng)計(jì)函數(shù)

summary()
library(Hmisc)  #加載Hmisc包
describe(mydata)  #調(diào)用函數(shù)

該函數(shù)統(tǒng)計(jì)樣本量穷当、缺失樣本量提茁、樣本唯一值數(shù)、均值馁菜、百分之五位數(shù)茴扁、百分之十位數(shù)、四分之一位數(shù)汪疮、中位數(shù)峭火、四分之三位數(shù)毁习、百分之九十位數(shù)、百分之九十五位數(shù)和最小卖丸、最大五個(gè)數(shù)值蜓洪，分別為：n, nmiss, unique, mean, 5,10,25,50,75,90,95th percentiles, 5 lowest and 5 highest scores

library(psych)      #加載psych安裝包
describe(mydata)    #調(diào)用describe函數(shù)

該函數(shù)統(tǒng)計(jì)變量名、變量序號(hào)坯苹、有效樣本量隆檀、均值、樣本標(biāo)準(zhǔn)差粹湃、中位數(shù)恐仑、眾數(shù)、最小值为鳄、最大值裳仆、偏度、峰度孤钦、標(biāo)準(zhǔn)誤差歧斟，分別為：item name ,item number, nvalid, mean, sd, median, mad, min, max, skew, kurtosis, se

library(pastecs) 
stat.desc（x,basic=TRUE,desc=TRUE,norm=FALSE,p=0.95）

數(shù)據(jù)框或者時(shí)間序列的所有值、空值偏形、缺失值的數(shù)量静袖，以及最小值、最大值俊扭、值域队橙、總和、均值萨惑、中位數(shù)捐康、標(biāo)準(zhǔn)誤、平均數(shù)值信度為95%的置信區(qū)間庸蔼、方差解总、標(biāo)準(zhǔn)差以及變異系數(shù)#nbr.val, nbr.null, nbr.na, min max, range, sum

分組統(tǒng)計(jì)函數(shù)

library(psych) #加載安裝包
describeBy(data, newdata$Gender) #分組統(tǒng)計(jì)

頻數(shù)統(tǒng)計(jì)

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118185.png" alt="image-20211223121610012" style="zoom: 50%;" />

分析各年齡層所占人數(shù)

#綁定數(shù)據(jù)
attach(data)
#一維列聯(lián)表
table(年齡)

分析各年齡層性別比例

#二維列聯(lián)表
table(年齡,性別)
#以年齡、性別排列的結(jié)果頻數(shù)三維列聯(lián)表
ftable(年齡,性別,結(jié)果)

交叉表

頻數(shù)交叉表

mytable <- table(A, B, C) 
ftable(mytable)

調(diào)包實(shí)現(xiàn)

library(gmodels)
CrossTable(A,B)

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118827.png" alt="image-20211223122235029" style="zoom:50%;" /><img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118393.png" alt="image-20211223122305419" style="zoom:50%;" /><img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118058.png" alt="image-20211223122403684" style="zoom:50%;" />

繪制多元分析圖

多元分析（下）P13

#安裝ggplot包
library(ggplot2)
ggplot(data)+geom_type()

繪制條形圖

ggplot(BOD) + geom_bar(aes(x=Time,y=demand),stat="identity")

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118353.png" alt="image-20211225225229238" style="zoom:50%;" />

將x軸由連續(xù)值變成不連續(xù)值

#factor()函數(shù)
BOD$Time<-factor(BOD$Time)
#將BOD中time列變成分類變量姐仅，就不再是連續(xù)值了

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118498.png" alt="image-20211225225342358" style="zoom:50%;" />

簇狀條形圖和堆積面積圖

簇狀條形圖和堆積面積圖包含兩個(gè)定性變量（分類變量）以及一個(gè)定量變量的情況花枫，用fill函數(shù)將不同類別填充為不同的顏色

ggplot(cabbage_exp)+geom_bar(aes(x=Cultivar,y=Weight,fill=date),stat=“identity”)

ggplot(cabbage_exp)+geom_bar(aes(x=Cultivar,y=Weight,fill=date),stat=“identity”, position=“dodge”) #dodge表示并排顯示 

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118248.png" alt="image-20211225230057715" style="zoom:50%;" /><img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118693.png" alt="image-20211225230106487" style="zoom:50%;" />

餅圖

ggplot(mydata,aes(x=factor(1),fill=Loan.Purpose))+geom_bar()+coord_polar(theta = "y")

factor(1) 是虛擬x軸的變量，原本的y軸被映射到極坐標(biāo)的theta萍嬉，因此y沒有實(shí)際的意義

散點(diǎn)圖

ggplot(data, aes(x=wt, y=mpg)) + geom_point(color="blue")

折線圖

多元分析（下）P25

ggplot(data,aes(x=wt,y=mpg)) + geom_line()

分組折線圖

ggplot(mtcars,aes(x=wt,y=mpg,group=factor(gear)))+geom_line(aes(color=factor(gear)))
#將gear變成因子變量

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118815.png" alt="image-20211225230957730" style="zoom:50%;" />

頻數(shù)直方圖

'''
binwidth = binsize 用來設(shè)置直方圖條形的寬度乌昔，組距越寬，直方圖越粗糙壤追。
#fill="pink" 填充的顏色
colour="blue" 邊緣的顏色
'''

mydata <- read.csv("F:/R/data/credit risk.csv")
binsize <- diff(range((mydata$Checking)))/20 #將x取值切分為20組
ggplot(mydata,aes(x=Checking)) + geom_histogram(binwidth = binsize, fill="pink", colour = "blue")

概率密度估計(jì)圖

ggplot(mtcars,aes(x=wt)) + geom_density()

多元線性回歸

讀取數(shù)據(jù)

rdata <- read.csv("F:/Rdata/crime.csv")
library(xlsx)
rdata <- read.xlsx("F:/Rdata/crime.xlsx")

基本描述統(tǒng)計(jì)

#匯總分析
summary(rdata)
library(psych)
a=describe(rdata[,c("crime","poverty","single")])
# 散點(diǎn)圖 
pairs(~crime+poverty+single,data=rdata, main="Simple Scatterplot Matrix")
#or
ggplot(data, aes(x=wt, y=mpg)) + geom_point(color="blue")

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171118750.png" alt="image-20211225233557330" style="zoom:70%;" />

建立回歸模型并估計(jì)參數(shù)

ols<- lm(crime ~ poverty + single, data = rdata))
summary(ols)

image-20211225233930765

對(duì)模型擬合效果進(jìn)行分析

opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
plot(ols, las = 1)

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119793.png" alt="image-20211225234126308" style="zoom:80%;" />

查看偏離較大的數(shù)據(jù)

P 27

library(MASS)
d1 <- cooks.distance(ols)          #計(jì)算偏離程度           
r <- stdres(ols)                   #標(biāo)準(zhǔn)化回歸誤差
a <- cbind(rdata, d1, r)           #綁定在一起磕道，進(jìn)行分析
a[d1 > 4/51, ]                     #前4個(gè)數(shù)據(jù)

庫克距離：庫克距離用來判斷強(qiáng)影響點(diǎn)是否為Y的異常值點(diǎn)。一般認(rèn)為行冰，當(dāng)D<0.5時(shí)認(rèn)為不是異常值點(diǎn)溺蕉；當(dāng)D>0.5時(shí)認(rèn)為是異常值點(diǎn)

在判斷Cook距離大小的時(shí)候伶丐，通常采用的經(jīng)驗(yàn)分界點(diǎn)是Cook距離序列的4/n處，其中n是觀測(cè)值的個(gè)數(shù)疯特。

討論模型是否可以通過robust回歸修正

Huber方法

殘差較小的觀測(cè)值被賦予的權(quán)重為1哗魂，殘差較大的觀測(cè)值的權(quán)重隨著殘差的增大而遞減

library(MASS)
summary(rr.huber <- rlm(crime ~ poverty + single, data = rdata))

也可以查看加權(quán)情況

hweights <- data.frame(state = rdata$state, resid = rr.huber$resid, weight = rr.huber$w) 
hweights2 <- hweights[order(rr.huber$w), ] 
hweights2[1:15, ]

Bisquare方法

所有的非0殘差所對(duì)應(yīng)觀測(cè)值的權(quán)重都是遞減的

rr.bisquare <- rlm(crime ~ poverty + single, data=rdata, psi = psi.bisquare) 
summary(rr.bisquare)

也可以查看加權(quán)情況

biweights <- data.frame(state = cdata$state, resid = rr.bisquare$resid, weight = rr.bisquare$w) 
biweights2 <- biweights[order(rr.bisquare$w), ]
biweights2[1:15, ]

Huber方法的軟肋在于無法很好的而處理極端離群點(diǎn)，而bisquare方法的軟肋在于回歸結(jié)果不易收斂漓雅，以至于經(jīng)常有多個(gè)最優(yōu)解

逐步回歸

#先建立線性回歸模型
steplm<-lm(crime~pctmetro+pctwhite+pcths+poverty+single,data=rdata)
#逐步回歸
step(steplm, direction = c("both", "backward", "forward"))

the mode of stepwise search, can be one of "both", "backward", or "forward", with a default of "both". If the scope argument is missing the default for direction is "backward".

logitstic回歸

函數(shù)形式

logistic回歸系數(shù)的意義

p 20 例二全面一點(diǎn)p 31

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119737.png" alt="image-20211228184211718" style="zoom:67%;" />

讀取數(shù)據(jù)

d1=read.table(“clipboard”,header=T)  #讀取例1數(shù)據(jù) 

描述統(tǒng)計(jì)

同上

建立全變量logistic回歸模型

logit.glm<-glm(y~x1+x2+x3,family=binomial,data=d1)  #Logistic回歸模型
summary(logit.glm)  #Logistic回歸模型結(jié)果

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119632.png" alt="image-20211228185650802" style="zoom:80%;" />

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119214.png" alt="image-20211228185755703" style="zoom:70%;" />

AIC為模型擬合效果

得到初步的logistic回歸模型：

逐步篩選變量logistic回歸模型：

logit.step<-step(logit.glm,direction="both")  #逐步篩選法變量選擇
summary(logit.step) #逐步篩選法變量選擇結(jié)果

得到最佳邏輯回歸模型

預(yù)測(cè)

pre1<-predict(logit.step, data.frame(x1=1))  #預(yù)測(cè)視力正常司機(jī)Logistic回歸結(jié)果
p1<-exp(pre1)/(1+exp(pre1))  #預(yù)測(cè)視力正常司機(jī)發(fā)生事故概率
pre2<-predict(logit.step,data.frame(x1=0))  #預(yù)測(cè)視力有問題的司機(jī)Logistic回歸結(jié)果
p2<-exp(pre2)/(1+exp(pre2))  #預(yù)測(cè)視力有問題的司機(jī)發(fā)生事故概率
c(p1,p2)  #結(jié)果顯示

多類別logitstic回歸

需要r包

library(foreign) 
library(nnet) 
library(ggplot2)
library(reshape2)

讀取數(shù)據(jù)

同上

描述統(tǒng)計(jì)

同上

設(shè)定參照類录别，做多分類logistic回歸分析

Library(nnet)
#設(shè)定參照類,新生成一列
ml$prog2 <- relevel(ml$prog, ref = "academic")
#做多分類logistic回歸
test <- multinom(prog2 ~ ses + write, data = ml)

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119142.png" alt="image-20211227222048635" style="zoom:70%;" />

獲取多分類logistic回歸分析結(jié)果

summary(test)

image-20211227223556328

• b13 變量 write 增加一個(gè)單位與參加普通課程與學(xué)術(shù)課程的對(duì)數(shù)幾率降低 0.058 相關(guān)。
• b23 變量 write 增加一個(gè)單位與參加職業(yè)計(jì)劃與學(xué)術(shù)計(jì)劃的對(duì)數(shù)幾率降低有關(guān)邻吞。 0.1136 的數(shù)量组题。
• b12 如果從 ses="low" 變?yōu)?ses="high"，則在普通課程與學(xué)術(shù)課程的對(duì)數(shù)幾率將降低 1.163抱冷。
• b11 如果從 ses="low" 移動(dòng)到 ses="middle"崔列，則普通課程與學(xué)術(shù)課程的對(duì)數(shù)幾率將降低 0.533，盡管該系數(shù)不顯著旺遮。
• b22 如果從 ses="low" 變?yōu)?ses="high"赵讯，則參加職業(yè)計(jì)劃與學(xué)術(shù)計(jì)劃的對(duì)數(shù)幾率將降低 0.983。
• b21 如果從 ses="low" 移動(dòng)到 ses="middle"耿眉，職業(yè)計(jì)劃與學(xué)術(shù)計(jì)劃的對(duì)數(shù)幾率將增加 0.291边翼，盡管該系數(shù)不顯著。

計(jì)算標(biāo)準(zhǔn)化得分 z-score及p-value

z <- summary(test)$coefficients/summary(test)$standard.errors 
z
p <- (1 - pnorm(abs(z), 0, 1)) * 2 
p

基本原理如下：

原假設(shè) H0：βj=0

備擇假設(shè) H1：βj≠0

構(gòu)造統(tǒng)計(jì)量：

計(jì)算p值跷敬，即z_score發(fā)生的概率讯私，雙側(cè)檢驗(yàn)

提取參數(shù)估計(jì)值并轉(zhuǎn)換為數(shù)值

exp(coef(test))

image-20211227225400405

ses社會(huì)地位變量從1變?yōu)?時(shí)，普通項(xiàng)目vs學(xué)術(shù)項(xiàng)目的相對(duì)優(yōu)勢(shì)比是0.3126

Write變量增加一個(gè)單位西傀，普通項(xiàng)目vs學(xué)術(shù)項(xiàng)目的相對(duì)優(yōu)勢(shì)比是0.9437

查看各樣本所屬分類的概率

head(pp <- fitted(test))

預(yù)測(cè)

dses <- data.frame(ses = c("low", "middle","high"),write=mean(ml$write))
predict(test, newdata = dses, "probs")
#預(yù)測(cè)：Write不變，檢測(cè)社會(huì)地位的改變對(duì)預(yù)測(cè)值的影響

泊松回歸

研究因變量為計(jì)數(shù)變量桶癣，自變量為數(shù)值變量或分類變量之間關(guān)系的回歸模型

require(ggplot2) 
require(sandwich)

讀取數(shù)據(jù)

p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")

描述統(tǒng)計(jì)

同上

p <- within(p, { prog <- factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational"))id <- factor(id)  })
summary(p)
with(p, tapply(num_awards, prog, function(x) {sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))}))

畫圖

ggplot(p, aes(num_awards, fill = prog))+geom_histogram(binwidth=.5, position="dodge")

建立回歸模型

summary(m1 <- glm(num_awards ~ prog + math, family="poisson", data=p))

計(jì)算p值及置信區(qū)間

library(sandwich)
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)
r.est

模型的擬合優(yōu)度拥褂，通過卡方統(tǒng)計(jì)量來判斷模型的擬合優(yōu)度

with(m1, cbind(res.deviance = deviance, df = df.residual,p = pchisq(deviance, df.residual, lower.tail=FALSE)))

預(yù)測(cè)

p$phat <- predict(m1, type="response")
# order by program and then by math
p <- p[with(p, order(prog, math)), ]

主成分分析

導(dǎo)入數(shù)據(jù) 計(jì)算相關(guān)矩陣

X=read.table("clipboard",header=T) 
cor(X)

求相關(guān)矩陣的特征值和主成分負(fù)荷

PCA = princomp(X,cor=T)  #主成分分析
PCA                  #特征值開根號(hào)結(jié)果 
summary(PCA)
PCA$loadings  #主成分載荷

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119244.png" alt="image-20211228200952960" style="zoom:80%;" />

最后一行為累計(jì)貢獻(xiàn)率，一般取累計(jì)貢獻(xiàn)率達(dá)到85%~95%的特征值牙寞。

image-20211228201428843

確定主成分

screeplot(PCA,type="lines") #繪制碎石圖

按照累積方差貢獻(xiàn)率大于80%原則饺鹃，選入了兩個(gè)主成分，其累積方差貢獻(xiàn)率為80.7%间雀，從碎石圖上也可以看出m取2比較合適悔详。

主成分得分

PCA$scores[,1:2]  #主成分得分

image-20211228201658597

主成分分析與多元線性回歸結(jié)合

P 32

導(dǎo)入數(shù)據(jù)

conomy<-data

同上

建立多元回歸模型

lm.sol<-lm(y~x1+x2+x3, data=conomy)
summary(lm.sol)

主成分分析

b<-data.frame(conomy[,1:3])  #對(duì)三個(gè)自變量作主成分分析
conomy.pr<-princomp(b, cor=T)
summary(conomy.pr, loadings=TRUE)

image-20211228205605313

將原始數(shù)據(jù)準(zhǔn)換為用主成分變量表達(dá)

上一步可知主成分為comp 1和comp 2，將第1主成分的預(yù)測(cè)值和第2主成分的預(yù)測(cè)值存放在數(shù)據(jù)框里

pre<-predict(conomy.pr)
conomy$z1<-pre[,1]
conomy$z2<-pre[,2]

作主成分回歸

lm.sol<-lm(y~z1+z2, data=conomy)
summary(lm.sol)

這是因變量與主成分的關(guān)系惹挟，應(yīng)用起來不方便茄螃，需變換成因變量與自變量的關(guān)系

變換公式代碼：

beta<-coef(lm.sol)       #提取回歸系數(shù)
A<-loadings(conomy.pr)     #提取主成分對(duì)應(yīng)的特征向量
x.m<-conomy.pr$center        #數(shù)據(jù)的中心,即各類數(shù)據(jù)均值
x.sd<-conomy.pr$scale       #數(shù)據(jù)的標(biāo)準(zhǔn)差         
coef<-(beta[2]*A[,1]+beta[3]*A[,2])/x.sd #
beta0<-beta[1]-sum(x.m*coef)
c(beta0,coef)

因子分析

P 26

讀入數(shù)據(jù)

X=read.table("clipboard",header=T)

數(shù)據(jù)標(biāo)準(zhǔn)化

data<-data.frame(X[,3:9])
X<-scale(data)

計(jì)算相關(guān)系數(shù)矩陣

cor(X)

計(jì)算特征值、因子載荷及共同度

(FA0=factanal(X,3,rot="none")) #極大似然法因子分析

library(mvstats) 
(Fac=factpc(X,3)) #主成份法因子分析

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119697.png" alt="image-20211228213947051" style="zoom:80%;" />

因子旋轉(zhuǎn)

如果求出主因子后连锯，各個(gè)主因子的典型代表變量不很突出归苍，還需要進(jìn)行因子旋轉(zhuǎn)用狱，通過適當(dāng)?shù)男D(zhuǎn)得到比較滿意的主因子。

進(jìn)行因子旋轉(zhuǎn)拼弃，就是要使因子載荷矩陣中因子載荷的絕對(duì)值向0和1兩個(gè)方向分化夏伊，使大的載荷更大，小的載荷更小

(Fa1=factanal(X,3,rot="varimax")) #varimax法旋轉(zhuǎn)因子分析

image-20211228213917398

因子得分

#使用回歸估計(jì)法的極大似然法因子分析 
Fa1=factanal(X,3,scores="regression") 
Fa1$scores

<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119634.png" alt="image-20211228214348570" style="zoom:67%;" />

或

#使用回歸估計(jì)法的主成份法因子分析 
Fac1=factpc(X,3,scores="regression") 
Fac1$scores