數(shù)據(jù)建模與預(yù)測期末復(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(),表示在第幾個分量求到其值
等差數(shù)列
x<-1:30
等差數(shù)列運(yùn)算大于加減乘除運(yùn)算
等間隔數(shù)列
seq(from,to,by = 間隔)
檢測是否為空值并替換
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”
多維數(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表示矩陣中對應(yīng)的元素的乘積，A％＊％B表示通常意義下的兩個矩陣的成積（要求矩陣A的列數(shù)等于矩陣B的行數(shù)
生成對角陣和矩陣取對角陣運(yùn)算
函數(shù)diag()依賴于它的變量：
當(dāng)v是一個向量時，diag(v)表示以v的元素為對角線元素的對角陣.當(dāng)M是一個矩陣時，則diag(M)表示的是取M對角線上的元素的向量.
解線性方程組和求矩陣的逆矩陣
#解方程組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存放著對稱矩陣Sm特征值和特征向量，是由列表形式給出的，其中evvectors是Sm的特征向量構(gòu)成的矩陣.
矩陣合并
#求矩陣的合并 函數(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)
數(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))
數(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)的符號，默認(rèn)為.
讀取excel/csv
X<-read.xlsx("F:\\OneDrive - chenxiangxi\\桌面\\2.xlsx",sheetIndex = 1) #安裝xlsx
X<-read.csv("textdata.csv")
儲存文件
#存儲到csv
write.csv(mydata,"G:/xxx/Inventory.csv")
#存儲數(shù)據(jù)到txt文件,制表符分隔
write.table(mydata,""G:/xxx/Inventory.txt", sep="\t")
#存儲數(shù)據(jù)到xlsx文件
write.xlsx(mydata, "G:/xxx/Inventor..xlsx", sheetName="Sheet1")
多元數(shù)據(jù)直觀表示
均值
mean(x,trim=0,na.rm=FALSE)
x是對象，trim是在計算均值前去掉x兩端觀察值的比例，默認(rèn)值為0，即包括全部數(shù)據(jù)。當(dāng)na. rm＝TRUE時，允許數(shù)據(jù)中有缺失數(shù)據(jù)。
求矩陣平均值
用apply()函數(shù)
apply(x, 1, mean) #計算矩陣各行的均值
apply(x, 2, mean) #計算矩陣各列的均值
方差
var(x)計算樣本方差
sd(x)計算樣本標(biāo)準(zhǔn)差
變異系數(shù)CV
當(dāng)CV>15%時，說明數(shù)據(jù)離散程度過高，需要注意
cv<-100*sd(x)/mean(x)
偏度
偏度系數(shù)是刻畫數(shù)據(jù)的對稱性指標(biāo)。關(guān)于均值對稱的數(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ù)近似為0
當(dāng)分布較正態(tài)分布的尾部更分散時，峰度系數(shù)為正；否則為負(fù)
當(dāng)峰度系數(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)計函數(shù)
summary()
library(Hmisc) #加載Hmisc包
describe(mydata) #調(diào)用函數(shù)
該函數(shù)統(tǒng)計樣本量、缺失樣本量、樣本唯一值數(shù)、均值、百分之五位數(shù)、百分之十位數(shù)、四分之一位數(shù)、中位數(shù)、四分之三位數(shù)、百分之九十位數(shù)、百分之九十五位數(shù)和最小、最大五個數(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)計變量名、變量序號、有效樣本量、均值、樣本標(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ù)、標(biāo)準(zhǔn)誤、平均數(shù)值信度為95%的置信區(qū)間、方差、標(biāo)準(zhǔn)差以及變異系數(shù)#nbr.val, nbr.null, nbr.na, min max, range, sum
分組統(tǒng)計函數(shù)
library(psych) #加載安裝包
describeBy(data, newdata$Gender) #分組統(tǒng)計
頻數(shù)統(tǒng)計
<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)包實現(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%;" />
簇狀條形圖和堆積面積圖
簇狀條形圖和堆積面積圖包含兩個定性變量（分類變量）以及一個定量變量的情況，用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沒有實際的意義
散點(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")
概率密度估計圖
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)計
#匯總分析
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%;" />
建立回歸模型并估計參數(shù)
ols<- lm(crime ~ poverty + single, data = rdata))
summary(ols)
對模型擬合效果進(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) #計算偏離程度
r <- stdres(ols) #標(biāo)準(zhǔn)化回歸誤差
a <- cbind(rdata, d1, r) #綁定在一起，進(jìn)行分析
a[d1 > 4/51, ] #前4個數(shù)據(jù)
庫克距離：庫克距離用來判斷強(qiáng)影響點(diǎn)是否為Y的異常值點(diǎn)。一般認(rèn)為，當(dāng)D<0.5時認(rèn)為不是異常值點(diǎn)；當(dāng)D>0.5時認(rèn)為是異常值點(diǎn)
在判斷Cook距離大小的時候，通常采用的經(jīng)驗分界點(diǎn)是Cook距離序列的4/n處，其中n是觀測值的個數(shù)。
討論模型是否可以通過robust回歸修正
Huber方法
殘差較小的觀測值被賦予的權(quán)重為1，殘差較大的觀測值的權(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殘差所對應(yīng)觀測值的權(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)常有多個最優(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)計
同上
建立全變量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ù)測
pre1<-predict(logit.step, data.frame(x1=1)) #預(yù)測視力正常司機(jī)Logistic回歸結(jié)果
p1<-exp(pre1)/(1+exp(pre1)) #預(yù)測視力正常司機(jī)發(fā)生事故概率
pre2<-predict(logit.step,data.frame(x1=0)) #預(yù)測視力有問題的司機(jī)Logistic回歸結(jié)果
p2<-exp(pre2)/(1+exp(pre2)) #預(yù)測視力有問題的司機(jī)發(fā)生事故概率
c(p1,p2) #結(jié)果顯示
多類別logitstic回歸
需要r包
library(foreign)
library(nnet)
library(ggplot2)
library(reshape2)
讀取數(shù)據(jù)
同上
描述統(tǒng)計
同上
設(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)
- b13 變量 write 增加一個單位與參加普通課程與學(xué)術(shù)課程的對數(shù)幾率降低 0.058 相關(guān)。
- b23 變量 write 增加一個單位與參加職業(yè)計劃與學(xué)術(shù)計劃的對數(shù)幾率降低有關(guān)。 0.1136 的數(shù)量。
- b12 如果從 ses="low" 變?yōu)?ses="high"，則在普通課程與學(xué)術(shù)課程的對數(shù)幾率將降低 1.163。
- b11 如果從 ses="low" 移動到 ses="middle"，則普通課程與學(xué)術(shù)課程的對數(shù)幾率將降低 0.533，盡管該系數(shù)不顯著。
- b22 如果從 ses="low" 變?yōu)?ses="high"，則參加職業(yè)計劃與學(xué)術(shù)計劃的對數(shù)幾率將降低 0.983。
- b21 如果從 ses="low" 移動到 ses="middle"，職業(yè)計劃與學(xué)術(shù)計劃的對數(shù)幾率將增加 0.291，盡管該系數(shù)不顯著。
計算標(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)計量：
計算p值，即z_score發(fā)生的概率，雙側(cè)檢驗
提取參數(shù)估計值并轉(zhuǎn)換為數(shù)值
exp(coef(test))
ses社會地位變量從1變?yōu)?時，普通項目vs學(xué)術(shù)項目的相對優(yōu)勢比是0.3126
Write變量增加一個單位，普通項目vs學(xué)術(shù)項目的相對優(yōu)勢比是0.9437
查看各樣本所屬分類的概率
head(pp <- fitted(test))
預(yù)測
dses <- data.frame(ses = c("low", "middle","high"),write=mean(ml$write))
predict(test, newdata = dses, "probs")
#預(yù)測：Write不變，檢測社會地位的改變對預(yù)測值的影響
泊松回歸
研究因變量為計數(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)計
同上
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))
計算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)計量來判斷模型的擬合優(yōu)度
with(m1, cbind(res.deviance = deviance, df = df.residual,p = pchisq(deviance, df.residual, lower.tail=FALSE)))
預(yù)測
p$phat <- predict(m1, type="response")
# order by program and then by math
p <- p[with(p, order(prog, math)), ]
主成分分析
導(dǎo)入數(shù)據(jù) 計算相關(guān)矩陣
X=read.table("clipboard",header=T)
cor(X)
求相關(guān)矩陣的特征值和主成分負(fù)荷
PCA = princomp(X,cor=T) #主成分分析
PCA #特征值開根號結(jié)果
summary(PCA)
PCA$loadings #主成分載荷
<img src="https://gitee.com/ayase314/polaris_pic/raw/master/202201171119244.png" alt="image-20211228200952960" style="zoom:80%;" />
最后一行為累計貢獻(xiàn)率，一般取累計貢獻(xiàn)率達(dá)到85%~95%的特征值。
依次可以寫出每個主成分的構(gòu)成形式
確定主成分
screeplot(PCA,type="lines") #繪制碎石圖
按照累積方差貢獻(xiàn)率大于80%原則，選入了兩個主成分，其累積方差貢獻(xiàn)率為80.7%，從碎石圖上也可以看出m取2比較合適。
主成分得分
PCA$scores[,1:2] #主成分得分
主成分分析與多元線性回歸結(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]) #對三個自變量作主成分分析
conomy.pr<-princomp(b, cor=T)
summary(conomy.pr, loadings=TRUE)
將原始數(shù)據(jù)準(zhǔn)換為用主成分變量表達(dá)
上一步可知主成分為comp 1和comp 2，將第1主成分的預(yù)測值和第2主成分的預(yù)測值存放在數(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) #提取主成分對應(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)
計算相關(guān)系數(shù)矩陣
cor(X)
計算特征值、因子載荷及共同度
(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)
如果求出主因子后，各個主因子的典型代表變量不很突出，還需要進(jìn)行因子旋轉(zhuǎn)，通過適當(dāng)?shù)男D(zhuǎn)得到比較滿意的主因子。
進(jìn)行因子旋轉(zhuǎn)，就是要使因子載荷矩陣中因子載荷的絕對值向0和1兩個方向分化，使大的載荷更大，小的載荷更小
(Fa1=factanal(X,3,rot="varimax")) #varimax法旋轉(zhuǎn)因子分析
因子得分
#使用回歸估計法的極大似然法因子分析
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%;" />
或
#使用回歸估計法的主成份法因子分析
Fac1=factpc(X,3,scores="regression")
Fac1$scores
