标签：分析 Affairs pwr level 回归 library sig 线性 table

多元线性回归 回归诊断 方差分析 功效分析 广义线性模型 logistics回归 主成分分析 因子分析 购物篮分析

多元线性回归

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])
#将矩阵转换为数据框
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)
coef(fit)
library(car)
qqPlot(fit, labels=row.names(states), id.method="identify",
       simulate=TRUE, main="Q-Q Plot")


#Mutiple linear regression with a significant interaction term
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)#查看拟合结果
#说明马力和车重的交互项显著
library(effects)
plot(effect("hp:wt", fit,, list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)
#AIC
fit1 <- lm (Murder ~ Population+Illiteracy+Income+Frost,data=states)
fit2 <- lm (Murder ~ Population+Illiteracy,data=states)
AIC(fit1,fit2)
#使用AIC来检测模型,第一个模型包含四个自变量，第二个模型包含两个自变量
#Backward stepwise selection
library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,
          data=states)
stepAIC(fit, direction="backward")


#All subsets regression
library(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income +
                     Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")

回归诊断

opar <- par(no.readonly=TRUE)
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)#生成评价拟合模型的四幅图
par(opar)#对绘图选项进行限制

fit2 <- lm(weight ~ height + I(height^2), data=women)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))#横排两张竖排两张
plot(fit2)
#第一幅图是残差与拟合的图，用来表示因变量与自变量是否呈线性关系
#图中的点是残差分布，曲线是拟合曲线
#第二幅图用来描述正态性，正态分布情况下，应该是一条直线
#第三幅图是位置与尺寸图，用来描述同方差性，如果满足方差不变，图中的点
#随机分布
#第四幅图是用来判断离群点，高杠杆点
par(opar)

方差分析

?aov
#ANOVA :Analysis of Variance
#研究对结果有影响的变量
library(multcomp)
attach(cholesterol)
?cholesterol
table(trt)     
aggregate(response, by=list(trt), FUN=mean) 
#分组统计，可以看出药物E的效果最好
#aggregate(response, by=list(trt), FUN=sd) 
fit <- aov(response ~ trt,data =cholesterol)                                
summary(fit)
fit.lm <- lm(response~trt,data=cholesterol)
plot(fit)
#One-way ANCOVA
data(litter, package="multcomp")
attach(litter)
table(dose) 
aggregate(weight, by=list(dose), FUN=mean)
#分组统计平均数
fit <- aov(weight ~ gesttime + dose) 
#使用aov分析因变量weight与协变量gesttime与自变量dose的关系
summary(fit)

#Two way ANOVA
attach(ToothGrowth)
table(supp,dose)
aggregate(len, by=list(supp,dose), FUN=mean)
aggregate(len, by=list(supp,dose), FUN=sd)
class(ToothGrowth$dose)
dose <- factor(dose)
fit <- aov(len~ supp*dose,data=ToothGrowth)
summary(fit)

install.packages("HH")
library(HH)
interaction.plot(dose, supp, len, type="b", 
                 col=c("red","blue"), pch=c(16, 18),
                 main = "Interaction between Dose and Supplement Type")
#使用HH包中的interaction函数对结果可视化
#One-way MANOVA
library(MASS)
attach(UScereal)
shelf <- factor(shelf)#将shelf列转换为因子
y <- cbind(calories, fat, sugars)
aggregate(y, by=list(shelf), FUN=mean)
cov(y)
fit <- manova(y ~ shelf)
summary(fit)
summary.aov(fit)

功效分析

install.packages("pwr")
library(pwr)
#使用pwr包进行功效分析
#Linear Models
pwr.f2.test(u=3, f2=0.0769, sig.level=0.05, power=0.90)
#0.05显著性水平，0.9的功效
#ANOVA
pwr.anova.test(k=2,f=.25,sig.level=.05,power=.9)
#每一组需要86个样本
#t tests
pwr.t.test(d=.8, sig.level=.05,power=.9, type="two.sample",  
           alternative="two.sided")
pwr.t.test(n=20, d=.5, sig.level=.01, type="two.sample", 
           alternative="two.sided")
#Correlations
pwr.r.test(r=.25, sig.level=.05, power=.90, alternative="greater")

#Tests of proportions
pwr.2p.test(h=ES.h(.65, .6), sig.level=.05, power=.9, 
            alternative="greater")
#Chi-square tests
prob <- matrix(c(.42, .28, .03, .07, .10, .10), byrow=TRUE, nrow=3)
ES.w2(prob)
pwr.chisq.test(w=.1853, df=3 , sig.level=.05, power=.9)

广义线性模型

?glm
#library(glm)
#glm
#getwd()
#setwd("D:\\backup\\Ln\\RIA\\RData")
#使用glm进行广义线性模型分析
data(breslow.dat, package="robust")
names(breslow.dat)
summary(breslow.dat[c(6, 7, 8, 10)])

attach(breslow.dat)
#fit regression
fit <- glm(sumY ~ Base + Age + Trt, data=breslow.dat, family=poisson(link="log"))
summary(fit)

#interpret model parametersR
coef(fit)
exp(coef(fit))

Logistics回归

#使用连续型或类别型变量来预测二值型变量
data(Affairs, package="AER")
#install.packages("AER")
summary(Affairs)
table(Affairs$affairs)
#使用summary和table简单统计分析affair
prop.table(table(Affairs$affairs))#只统计affair这一项
prop.table(table(Affairs$gender))#统计gender项

#create binary outcome variable
Affairs$ynaffair[Affairs$affairs > 0] <- 1
#将年度出轨次数中大于0的数据赋值为1
Affairs$ynaffair[Affairs$affairs == 0] <- 0
#将年度出轨次数中大于0的数据赋值为0
head(Affairs)
Affairs$ynaffair <- factor(Affairs$ynaffair, 
                           levels=c(0,1),
                           labels=c("No","Yes"))
#将affair转换为因子，其中只取年度出轨数据，数据的水平为0和1，其中
#0是no，1是yes
table(Affairs$ynaffair)

#fit full model
attach(Affairs)
fit <- glm(ynaffair ~ gender + age + yearsmarried + children + 
                  religiousness + education + occupation +rating,
                data=Affairs,family=binomial())
#Affairs$ynaffair
#对年度出轨数据进行拟合，ynaffair为因变量，后面的为自变量
summary(fit)
#结果可以看到P值后面的*，*表示这个变量对affair影响是否显著
#fit reduced model
fit1 <- glm(ynaffair ~ age + yearsmarried + religiousness + 
                     rating, data=Affairs, family=binomial())
summary(fit1)
#在新模型中，每一个因素的影响都非常显著
#compare models
anova(fit, fit1, test="Chisq")
#使用卡方检验，证明了两次检验结果差别不大，那么我们就可以精简变量了
#interpret coefficients
coef(fit1)
#使用coef算出回归系数，其他预测变量不变时，1单位预测变量的变化可能引起
#响应变化对数优势比的变化
exp(coef(fit1))#由于对其使用了对数优势比，那么我们这里将其指数运算
#回到了正常模式
#结果可以看到，婚龄增加一年婚外情优势比增加1.11，年龄增加一年
#婚外情优势比增加0.97
#calculate probability of extramariatal affair by marital ratings
testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
                       age = mean(Affairs$age),
                       yearsmarried = mean(Affairs$yearsmarried),
                       religiousness = mean(Affairs$religiousness))
#使用测试数据集预测相应概率
testdata$prob <- predict(fit1, newdata=testdata, type="response")
testdata
#calculate probabilites of extramariatal affair by age
testdata <- data.frame(rating = mean(Affairs$rating),
                       age = seq(17, 57, 10), 
                       yearsmarried = mean(Affairs$yearsmarried),
                       religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit1, newdata=testdata, type="response")
testdata

主成分分析

#Principal components analysis of US Judge Ratings
library(psych)
#主成分分析，将大量相关变量转换为一组很少的不相关变量
head(USJudgeRatings)
fa.parallel(USJudgeRatings,fa="pc",n.iter = 100)
#作pc分析，循环100次，绘制碎石图
#碎石图用来确定使用几个因子比较恰当
pc <- principal(USJudgeRatings, nfactors=1)
pc
#使用principal进行分析，nfactors指定主成分数
#Principal components analysis Score
pc <- principal(USJudgeRatings,nfactors = 1,scores = TRUE)
pc$scores#获得每个变量的得分

#Principal components analysis Harman23.cor data
fa.parallel(Harman23.cor$cov, n.obs=302, fa="pc", n.iter=100,
            show.legend=FALSE, main="Scree plot with parallel analysis")
#继续绘制碎石图，n.obs表示样本大小
#Principal components analysis of body measurements
library(psych)
PC <- principal(Harman23.cor$cov, nfactors=2, rotate="none")
PC#nfactors=2表示有两个主成分
#Principal components analysis with varimax rotation
rc <- principal(Harman23.cor$cov, nfactors=2, rotate="varimax")
rc#主成分的旋转

因子分析

#因子分析法，本质上用来降维
options(digits=2)
library(psych)
covariances <- ability.cov$cov
#convert covariances to correlations
correlations <- cov2cor(covariances)
correlations

#determine number of factors to extract
fa.parallel(correlations, n.obs=112, fa="both", n.iter=100,
            main="Scree plots with parallel analysis")
#判断提取因子数

#Principal axis factoring without rotation
fa <- fa(correlations, nfactors=2, rotate="none", fm="pa")
fa


#Factor extraction with orthogonal rotation
fa.varimax <- fa(correlations, nfactors=2, rotate="varimax", fm="pa")
fa.varimax


#Listing Factor extraction with oblique rotation
fa.promax <- fa(correlations, nfactors=2, rotate="promax", fm="pa")
fa.promax

#plot factor solution
factor.plot(fa.promax, labels=rownames(fa.promax$loadings))
fa.diagram(fa.promax, simple=FALSE)

#factor scores
fa <- fa(correlations,nfactors=2,rotate="none",fm="pa",score=TRUE)
fa.promax$weights

购物篮分析

install.packages("arules")
library(arules)
data(Groceries)
Groceries#内置数据集
inspect(Groceries)#查看数据集内容
fit <- apriori(Groceries,parameter = list(support=0.01,confidence=0.5))
#使用apriori进行建模，最小支持度support为0.01，最小置信度confidence为0.5
summary(fit)
inspect(fit)	

标签：分析,Affairs,pwr,level,回归,library,sig,线性,table
来源： https://www.cnblogs.com/kwq717/p/14696797.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。