精华内容
下载资源
问答
  • R语言 线性回归分析

    2019-03-04 23:05:11
    回归分析是非常基础的但是也是极其泛用化的统计预测工具之一,而R语言也有自带的function来完成一元和多元回归。 针对一元回归,数据集选用的是R语言自带的women数据集 data(“women”) plot(x=womenweight,y=...

    回归分析是非常基础的但是也是极其泛用化的统计预测工具之一,而R语言也有自带的function来完成一元和多元回归。

    针对一元回归,数据集选用的是R语言自带的women数据集

    data(“women”)
    plot(x=women w e i g h t , y = w o m e n weight,y=women weight,y=womenheight)
    model=lm(height~weight,data = women)
    abline(model)

    这里应该注意,plot中x、y的数据应该与lm中的y~x相对应,不然abline是没法显示在图上的。

    summary(model)
    Call:
    lm(formula = weight ~ height, data = women)
    Residuals:
    Min 1Q Median 3Q Max
    -1.7333 -1.1333 -0.3833 0.7417 3.1167
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
    height 3.45000 0.09114 37.85 1.09e-14 ***
    Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
    Residual standard error: 1.525 on 13 degrees of freedom
    Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
    F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14

    summary可以展示包含最大值最小值、中位数四分位数及截距等多种基础数据,而且还包括基础决定系数和调整后决定系数等多种反映拟合程度的数据

    anova(model)
    Analysis of Variance Table
    Response: weight
    Df Sum Sq Mean Sq F value Pr(>F)
    height 1 3332.7 3332.7 1433 1.091e-14 ***
    Residuals 13 30.2 2.3
    —Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

    anova方差分析能够反映单个参数是否显著,还有其他一些信息。

    而其他信息可以用协方差cor()或者残差residuals()等拓展函数显示出来,也可以画图,这样更直观,比如

    plot(predict(model,type = “response”),residuals(model,type = “deviance”))
    plot(hatvalues(model))

    下面进行预测,这里随便按了两个数

    pre=data.frame(height=c(56,74))
    predict(model,pre,interval = ‘prediction’,level=0.95)
    fit lwr upr
    1 105.6833 101.847 109.5197
    2 167.7833 163.947 171.6197

    R的pre会给出一个区间,上下界及预测值都会给出来,如何取就看个人了。

    而对于多元线性回归,主要部分都与一元线性回归大同小异,只要在lm()中fomula参数位置的~后面多写几个变量就够了。而到底写几个变量,这里可以使用逐步回归来对变量进行筛选。

    mydata<-data.frame(
    x1=c( 7, 1,11,11, 7,11, 3, 1, 2,21, 1,11,10),
    x2=c(26,29,56,31,52,55,71,31,54,47,40,66,68),
    x3=c( 6,15, 8, 8, 6, 9,17,22,18, 4,23, 9, 8),
    x4=c(60,52,20,47,33,22, 6,44,22,26,34,12,12),
    Y =c(78.5,74.3,104.3,87.6,95.9,109.2,102.7,72.5,
    93.1,115.9,83.8,113.3,109.4)
    )
    model<-lm(Y~x1+x2+x3+x4,data=tdata)
    summary(model)
    Call:
    lm(formula = Y ~ x1 + x2 + x3 + x4, data = tdata)
    Residuals:
    Min 1Q Median 3Q Max
    -3.1750 -1.6709 0.2508 1.3783 3.9254
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 62.4054 70.0710 0.891 0.3991
    x1 1.5511 0.7448 2.083 0.0708 .
    x2 0.5102 0.7238 0.705 0.5009
    x3 0.1019 0.7547 0.135 0.8959
    x4 -0.1441 0.7091 -0.203 0.8441
    —Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
    Residual standard error: 2.446 on 8 degrees of freedom
    Multiple R-squared: 0.9824, Adjusted R-squared: 0.9736
    F-statistic: 111.5 on 4 and 8 DF, p-value: 4.756e-07

    从上面的t值判别可以看出,没有任何一项关于Y显著,所以我们需要运用逐步回归来优化这个回归模型。

    mystep=step(model)
    Start: AIC=26.94
    Y ~ x1 + x2 + x3 + x4
    Df Sum of Sq RSS AIC
    - x3 1 0.1091 47.973 24.974
    - - x4 1 0.2470 48.111 25.011
    - - x2 1 2.9725 50.836 25.728
    - 47.864 26.944
    - - x1 1 25.9509 73.815 30.576
    Step: AIC=24.97
    Y ~ x1 + x2 + x4
    Df Sum of Sq RSS AIC
    47.97 24.974
    - x4 1 9.93 57.90 25.420
    - - x2 1 26.79 74.76 28.742
    - - x1 1 820.91 868.88 60.629

    R的step函数判别标准是跟据AIC信息的大小来删除的,所以有时候会出现不太完整的逐步回归。
    比如这个例子

    summary(mystep)
    Call:
    lm(formula = Y ~ x1 + x2 + x4, data = tdata)
    Residuals:
    Min 1Q Median 3Q Max
    -3.0919 -1.8016 0.2562 1.2818 3.8982
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 71.6483 14.1424 5.066 0.000675 ***
    x1 1.4519 0.1170 12.410 5.78e-07 ***
    x2 0.4161 0.1856 2.242 0.051687 .
    x4 -0.2365 0.1733 -1.365 0.205395
    —Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
    Residual standard error: 2.309 on 9 degrees of freedom
    Multiple R-squared: 0.9823, Adjusted R-squared: 0.9764
    F-statistic: 166.8 on 3 and 9 DF, p-value: 3.323e-08

    可以看出变量并没有完全显著,这里可以用drop1()函数来进行后退法(与此相对的也有add1()函数),或者再step()中多写一个“k=2”

    step(model,trace = 1,k=2)
    Start: AIC=26.94
    Y ~ x1 + x2 + x3 + x4
    Df Sum of Sq RSS AIC
    - x3 1 0.1091 47.973 24.974
    - - x4 1 0.2470 48.111 25.011
    - - x2 1 2.9725 50.836 25.728
    - 47.864 26.944
    - - x1 1 25.9509 73.815 30.576
    Step: AIC=24.97
    Y ~ x1 + x2 + x4
    Df Sum of Sq RSS AIC
    47.97 24.974
    - x4 1 9.93 57.90 25.420
    - - x2 1 26.79 74.76 28.742
    - - x1 1 820.91 868.88 60.629
    Call:
    lm(formula = Y ~ x1 + x2 + x4, data = tdata)
    Coefficients:
    (Intercept) x1 x2 x4
    71.6483 1.4519 0.4161 -0.2365

    可以看删除X3 X4是比较好的选择

    data<-lm(Y~x1+x2,data=tdata)
    summary(data)
    Call:
    lm(formula = Y ~ x1 + x2, data = tdata)
    Residuals:
    Min 1Q Median 3Q Max
    -2.893 -1.574 -1.302 1.363 4.048
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 52.57735 2.28617 23.00 5.46e-10 ***
    x1 1.46831 0.12130 12.11 2.69e-07 ***
    x2 0.66225 0.04585 14.44 5.03e-08 ***
    Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1
    Residual standard error: 2.406 on 10 degrees of freedom
    Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744
    F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09

    而最后的summary显示确实是这样,删除3、4后,1、2便都显著了。

    展开全文
  • R语言 线性回归分析实例

    千次阅读 2019-09-22 16:07:39
    y,X1,X2,X3 分别表示第 t 年各项税收收入(亿元),某国生产...(1) 建立线性模型: ① 自己编写函数: > library(openxlsx) > data = read.xlsx("22_data.xlsx",sheet = 1) > x = data[,-c(1,2)] &g...

    y,X1,X2,X3 分别表示第 t 年各项税收收入(亿元),某国生产总值GDP(亿元),财政支出(亿元)和商品零售价格指数(%).

    (1) 建立线性模型

     ① 自己编写函数:

    > library(openxlsx)
    > data = read.xlsx("22_data.xlsx",sheet = 1)
    > x = data[,-c(1,2)]
    > x = cbind(rep(1,17),x)
    > x_mat  = as.matrix(x)
    > y =matrix(data[,2],ncol = 1)
    > res = solve(t(x_mat)%*%x_mat)%*%t(x_mat)%*%y
    > res
                        [,1]
    rep(1, 17) 19412.8597818
    X1             0.2679605
    X2            -0.2874013
    X3          -297.3653736
    

     所以各参数的估计值分别为 

     ② lm函数

    > lm(y~x_mat)
    
    Call:
    lm(formula = y ~ x_mat)
    
    Coefficients:
        (Intercept)  x_matrep(1, 17)          x_matX1  
    19412.859781545               NA      0.267960511  
            x_matX2          x_matX3  
       -0.287401287   -297.365373557  
    

     于是各参数的估计值分别为

    这两个方法的结果是一样的。

     

    (2)要求实验报告中画出矩阵散点图,给出参数的点估计、区间估计、t检验值、判定系数和模型F检验的方差分析表

    绘制矩阵散点图。

    library(graphics)
    pairs(data[,-1]pch = 21,bg = c('red','green3','blue'))
    # pch参数是控制点的形状,bg是控制点的颜色
    

     

    下面代码给出参数的点估计,t检验值,判定系数

    > summary(lm(y~x_mat+1))
    
    Call:
    lm(formula = y ~ x_mat + 1)  #调用
    
    Residuals:      #残差统计量,残差第一四分位数(1Q)和第三分位数(3Q)有大约相同的幅度,意味着有较对称的钟形分布
        Min      1Q  Median      3Q     Max  
    -4397.9 -1102.4   153.8  1184.4  2934.6 
    
    Coefficients: (1 not defined because of singularities)        
                      Estimate Std. Error t value Pr(>|t|)    
    (Intercept)      1.941e+04  3.524e+04   0.551    0.591    
    x_matrep(1, 17)         NA         NA      NA       NA    
    x_matX1          2.680e-01  4.466e-02   6.000 4.45e-05 ***
    x_matX2         -2.874e-01  1.668e-01  -1.723    0.109    
    x_matX3         -2.974e+02  3.688e+02  -0.806    0.435    
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    #标记为Estimate的列包含由最小二乘法计算出来的估计回归系数。
    #标记为Std.Error的列是估计的回归系数的标准误差。
    #从理论上说,如果一个变量的系数是0,那么该变量将毫无贡献。然而,这里显示的系数只是估计,它们不会正好为0.
    #因此,我们不禁会问:从统计的角度而言,真正的系数为0的可能性有多大?这是t统计量和P值的目的,在汇总中被标记为t value和Pr(>|t|) #P值估计系数不显著的可能性,有较大P值的变量是可以从模型中移除的候选变量

    Residual standard error: 2013 on 13 degrees of freedom Multiple R-squared: 0.9982, Adjusted R-squared: 0.9977 F-statistic: 2348 on 3 and 13 DF, p-value: < 2.2e-16

    #Residual standard error 表示残差的标准差,F-statistic 表示F的统计量

      区间估计?方差分析表?

     

    (3)保留模型中线性关系显著的预测变量确定最后的模型,并利用R软件中的"predict"语句预测2017年的税收收入

     根据回归分析结果,只有变量X1具有显著性。所以模型中仅保留变量X1。

    构造模型:

    x_mat = cbind(rep(1,17),data[,3]) 
    y = data[,2]
    res = lm(y~x_mat)
    res
    > res
    
    Call:
    lm(formula = y ~ x_mat)
    
    Coefficients:
    (Intercept)       x_mat1       x_mat2  
     -6213.0189           NA       0.1915 
    

     

     该模型为:Y = -6213.0189 + 0.1915 X1

    接下来预测2017年的税收收入,先根据数据data对 t 和 y 之间的关系进行回归分析

    t = data[,1]
    y = data[,2]
    res = lm(y~t)
    res
    
    > res
    
    Call:
    lm(formula = y ~ t)
    
    Coefficients:
    (Intercept)            t  
      -16428607         8213  

     所以 t 与 y 的关系为:y  = -16428607 + 8213 t

     预测 2017 年的税收收入:

    > newdata = data.frame(t = 2017)
    > pre = predict(res,newdata,interval = "prediction",level = 0.95)
    > pre
           fit      lwr      upr
    1 136337.8 116018.1 156657.4
    

      

     

    转载于:https://www.cnblogs.com/jiaxinwei/p/11523752.html

    展开全文
  • R语言线性回归分析

    千次阅读 2019-10-29 15:01:39
    线性回归


    基本假设:正态性、独立性、线性、同方差性
    常用R包:car(carData)、MASS、leap…

    1. 原始数据的分析

    plot(x,y); abline(lm(y~x))
    library(carData);library(car)
    scatterplot(x,y)  #生成拟合的二元关系图(散点+箱线)
    scatterplotMatrix()  #生成散点图矩阵
    
    cor()  #计算相关系数
    

    2. 回归模型的拟合(参数估计和检验)

    # 数据准备
    n <- ; p <- 
    Y <- as.matrix(dataset[,1])
    X1 <- as.matrix(dataset[,2:p])
    X <- cbind(1,X1)
    XX <- t(X)%*%X
    
    # 参数估计
    Beta <- solve(XX)%*%t(X)%*%Y
    e <- Y-X%*%Beta #<or> resid <- resid(lm) #残差的计算
    SSe <- t(e)%*%e
    sigma2 <- SSe/(n-p) || sigma2 <- sd(e)^2*(n-1)/(n-p)  #σ标准差的估计
    covBeta <- sigma2*solve(XX)  #Beta的协方差阵
    r <- matrix(nrow=p,ncol=p)
    for (i in 1:p)  {
    	for (j in 1:p)
    		r[i,j] <- covBeta[i,j]/(sqrt(covBeta[i,i])*sqrt(covBeta[j,j]))
    		}   #Beta的相关系数阵
    
    # 模型检验
    Beta1 <- beta[-1]
    Xc <- (X-matrix(1,nrow=18)%*%colMeans(X))[,-1]
    RSS_I <- t(Beta1)%*%t(Xc)%*%Y
    RSS_H <- n*mean(Y)^2
    SS_He <- t(Y)%*%Y-RSS_H
    RSS <- t(Beta)%*%t(X)%*%Y  ||  RSS <- RSS_H + RSS_I
    SSe <- t(Y)%*%Y-RSS
    F <- (RSS-RSS_H)/(p-1)/(SSe/(n-p))
    pf(F,p-1,n-p,lower.tail = T)
    
    # 系数检验
    t <- matrix(1,ncol=p)
    for(j in 1:4) t[j]<-beta[j]/sqrt(sigma2*solve(t(X)%*%X)[j,j])
    p-value <- 2*pt(abs(t),n-p,lower.tail = FALSE)
    
    # 复相关系数
    TSS<-t(Y-mean(Y))%*%(Y-mean(Y))
    R2 <- RSS_I  / TSS
    Adjusted_R2 <- 1-(1-R2)*(n-1)/(n-p)
    
    myfit <- lm(y~x1+x2+x3,dataframe)
    
    summary() 展示拟合模型的详细结果
    coefficients()  列出拟合模型的模型参数(截距项和斜率)
    confint()  列出模型参数的置信区间(95%fitted()  列出拟合模型的预测值
    anova()  生成拟合模型的方差分析表
    vcov()  列出模型参数的协方差矩阵
    AIC()  输出赤池信息统计量
    plot()  生成评价拟合模型的诊断图
    predict()  用拟合模型对新的数据集预测响应变量值
    

    3. 变量选择

    变量选择准则

    ## 平均残差平方和准则 (RMSq)
    RMSq <- SSeq/(n-q)
    ## Cp准则
    Cp <- SSeq/sigma2-(n-2*q)
    ## 赤池信息准则 (AIC/BIC)
    AIC <- (-2)*lnL+2*q
    BIC <- (-2)*lnL+q*log(n)
    
    # (下面代码假设有3个自变量)
    e <- list()
    SSeq<-matrix(1,nrow=7); RMS<-matrix(1,nrow=7); Cp<-matrix(1,nrow=7); 
    lnL<-matrix(1,nrow=7); AIC<-matrix(1,nrow=7); BIC<-matrix(1,nrow=7);
    q<-matrix(c(2,2,2,3,3,3,4),nrow=1)
    e[[7]] <- resid(lmlist[[7]])
    for(i in 1:7)
    {
      e[[i]] <- resid(lmlist[[i]])
      SSeq[i]<- t(e[[i]])%*%e[[i]]
      RMS[i]<-SSeq[i]/(n-q[i])
      Cp[i] <-SSeq[i]/(sd(e[[7]])^2*(n-1)/(n-p))-n+2*q[i]
      lnL[i]<-(-n/2)*(log(SSeq[i])+log(2*pi/n)+1)
      AIC[i]<-(-2)*lnL[i]+2*q[i]   ||  AIC(lmlist[[i]])   # 两者的结果有一点差异
      BIC[i]<-(-2)*lnL[i]+q[i]*log(n)  ||  BIC(lmlist[[i]])   # 两者的结果有一点差异
    } 
    COM<-data.frame(row.names=c('x1','x2','x3','x1 x2','x1 x3','x2 x3','x1 x2 x3'),RMS,Cp,AIC)
    COM
    

    逐步回归

    library(MASS)
    exps <- stepAIC(lm,scope= ,direction=<"both"/"forward"/"backward">...)
    summary(exps)
    

    全子集回归

    library(leaps)
    leaps <- regsubsets(y~.,data= ,nbest=<n>)   # nbest表示展示n个最佳的子集回归模型
    plot(leaps,scale="adjr2") #绘图展示子模型和调整R平方的结果
    coef(leaps,6)
    
    sleaps <- summary(leaps)
    names(sleaps)  # "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj" 
    plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
    i <- which.max(reg.summary$adjr2)
    points(i,reg.summary$adjr2[i], col="red",cex=2,pch=20)
    
    res <- data.frame(sleaps$outmat,调整R平方=sleaps$adjr2)
    res   #表格展示子模型和调整R平方的结果
    
    library(car)
    subsets(leaps,statistic="cp",main="Cp Plot for All subsets Regression")   #基于Cp统计量的不同子模型的比较
    abline(1,1,lty=2,col="red")
    

    向前/向后回归

    regfit.fwd=regsubsets(y~.,data=,nvmax=,method="forward")#向前逐步选择
    summary(regfit.fwd)
    regfit.bwd=regsubsets(y~.,data=,nvmax=,method="backward") #向后逐步选择
    summary(regfit.bwd)
    coef(regfit.fwd,7)
    coef(regfit.bwd,7)
    

    模型选择

    set.seed(1)
    train=sample(c(TRUE,FALSE), nrow(dataset),rep=TRUE) #rep是replace
    test=(!train)
    regfit.best=regsubsets(y~.,data=dataset[train,],nvmax=19)
    test.mat=model.matrix(y~.,data=dataset[test,])#此函数用于生成回归设计矩阵
    val.errors=rep(NA,19)
    for(i in 1:19){
      coefi=coef(regfit.best,id=i)
      pred=test.mat[,names(coefi)]%*%coefi
      val.errors[i]=mean((Hitters$Salary[test]-pred)^2)
    }
    val.errors
    which.min(val.errors)
    coef(regfit.best,10)
    
    #以上程序有点冗余,可以编写预测函数
    predict.regsubsets=function(object,newdata,id,...){
      form=as.formula(object$call[[2]])
      mat=model.matrix(form,newdata)
      coefi=coef(object,id=id)
      xvars=names(coefi)
      mat[,xvars]%*%coefi
    }
    regfit.best=regsubsets(y~.,data=dataset,nvmax=19)
    coef(regfit.best,10)
    
    # k折交叉验证
    k=10
    set.seed(1)
    folds=sample(1:k,nrow(dataset),replace=TRUE)
    cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
    for(j in 1:k){
     best.fit=regsubsets(y~.,data=dataset[folds!=j,],nvmax=19)
      for(i in 1:19){
        pred=predict.regsubsets(best.fit,dataset[folds==j,],id=i)#这里可以简化为predict
        cv.errors[j,i]=mean( (dataset$y[folds==j]-pred)^2)
      }
    }
    mean.cv.errors=apply(cv.errors,2,mean)
    mean.cv.errors
    par(mfrow=c(1,1))
    plot(mean.cv.errors,type='b')
    reg.best=regsubsets(y~.,data=dataset, nvmax=19)
    coef(reg.best,11)
    

    4. 回归诊断

    综合诊断

    # 残差图、QQ图、位置尺度图、cook统计量图
    plot(fit)
    par(mfrow=c(2,2)); plot(fit)  # 生成在一张图中
    
    # 模型假设通过与否的综合检验
    library(gvlma)
    gvmodel <- gvlma(fit)
    summary(gvmodel)
    

    对误差项假设的诊断

    ## 残差的计算
    resid() or residuals()   # 拟合模型的残差值(SSe)
    resid()/sd(resid())  #标准化残差
    rstudent()  # 学生化残差(ri)
    
    ## 正态性诊断
    # ①正态QQ图
    library(car)
    qqPlot()
    # ②学生化残差柱状图
    residplot <- function(fit,nbreaks=10)
    {
      z<-rstudent(fit)
      hist(z,breaks=nbreaks,freq=FALSE,xlab="Studentized Residual",main="Distribution of Errors")
      rug(jitter(z),col="brown")
      curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=2)
      lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2)
      legend("topright",legend=c("Normal Curve","Kernel Density Curve"),lty=1:2,col=c("blue","red"),cex=.7)
     }
    residplot(lm)
    
    ## 误差独立性诊断
    # Durbin-Watson检验 (p-value>0.05 → 无自相关性,误差项之间独立)
    durbinWatsonTest(lm)
    
    ## 线性诊断
    # 成分残差图(偏残差图)
    library(car)
    crPlots(lm)
    
    ## 同方差性诊断
    ncvTest(lm) #生成计分检验,零假设为误差方差不变
    spreadLevelPlot(lm) #生成最佳拟合曲线的散点图
    

    异常观测值的影响分析

    ## 异常值点(离群点——模型预测效果不佳的观测点,有很大的残差值)
    # ①QQ图——落在置信区间以外的点
    # ②标准化残差值大于2或小于-2的点
    # ③根据单个最大残差值的显著性判断模型是否存在离群点
    outlierTest()  # p-value<0.05, 判定为离群点
    
    ## 强影响点(对模型参数估计值影响比例失衡的点)
    # ①Cook距离(D统计量 > 4/(n-p-1))
    cutoff <- 4/(nrow(data)-length(fit$coefficients)-2)
    plot(fit,which=4,cook.levels=cutoff)
    abline(h=cutoff,lty=2,col="red")
    # ②变量添加图(提供关于强影响点如何影响模型的信息)
    library(car)
    avPlots(fit,ask=FALSE)
    
    ## 高杠杆值点(与其他预测变量有关的离群点,由许多异常的预测变量值组合起来,与响应变量没有关系)
    # 帽子统计量(hatvalue > 2or3*p/n)
    hat.plot <- function(fit){
    	p <- length(coefficients(fit))
    	n <- length(fitted(fit))
    	plot(hatvalues(fit),main="Index Plot of Hat Values")
    	abline(h=c(2,3)*p/n,col="red",lty=2)
    	idenfify(1:n,hatvalues(fit),names(hatvalues(fit)))
    	}  #identify选项可以在图中交互式标注离群点
    hat.plot(fit)
    
    ## 综合
    library(car)
    influencePlot(fit,main="Influence Plot",sub="Circle size is proportional to Cook's distance")
    # 纵坐标超过2或-2的可能是离群点,横坐标超过0.2或0.3的可能是高杠杆值点,圆圈较大的可能是强影响点
    

    多重共线性诊断

    ## 方差膨胀因子 vif=1/(1-R^2), vif越大,变量之间的线性相关程度越大
    #(一般vif>10表明存在严重多重共线性,剔除该变量;或计vif的均值>1)
    #当样本量不算小,而R2接近1时,可以认为存在严重的多重共线性
    library(car)
    vif(lm) 
    
    ## 条件数 k(lambda(1)/lambda(p-1)) 
    # (<=100:共线性的程度较小 100~1000:存在较强的多重共线性   >1000:存在严重的多重共线性)
    X <- scale(as.matrix(data[2:]));X
    kappa(t(X)%*%X,exact=TRUE)
    
    ## 另:计算矩阵的特征值
    lambda <- eigen(t(X)%*%X)
    lambda$values[1]/lambda$values[p-1]
    

    5. 改进措施

    1.删除观测点
    2.增删变量(删除有多重共线性的变量)
    3.Box-Cox变换(变换Y)

    bc<-boxcox(y~ ,lambda=seq(-2,2,0.01))
    
    lambda<-bc$x[which.max(bc$y)] 
    lambda
    
    y_bc<-(y^lambda-1)/lambda 
    lm_bc<-lm(y_bc~I(x1^2)+x2)
    summary(lm_bc)
    
    par(mfrow=c(2,2))
    plot(lm_bc)
    

    4.Box-Cox变换、Box-Tidwell变换(变换X)

    # Box-Cox变换(适用于模型违反正态假设的情况)
    library(car)
    summary(powerTranform(x))
    # Box-Tidwell变换(适用于模型违反线性假设的情况)
    library(car)
    boxTidwell(y~x1+x2)
    

    5.尝试其他方法(岭回归、稳健回归、非参数回归、非线性回归、时间序列、多层次回归、广义线性回归)

    6. 岭回归

    8. 其他

    有交互项的线性回归

    lm<-lm(y~x1+x2+x1:x2)
    library(effects)
    plot(effect("x1:x2",lm,,list(x1=c()),multiline=TRUE)
    

    深层次分析——交叉验证(刀切法)

    shrinkage <- function(fit,k=10){
    	  require(bootstrap)
    	  x <- fit$model[,2:ncol(fit$model)]
    	  y <- fit$model[,1]
    	  theta.fit <- function(x,y){lsfit(x,y)}
    	  theta.predict <- function(fit,x){as.matrix(cbind(1,x),nrow=18)%*%(as.matrix(fit$coef))}
    	  results <- crossval(t(as.matrix(x,nrow=18)),y,theta.fit,theta.predict,ngroup=k)
    	  r2 <- cor(y,fit$fitted.values)^2
    	  r2cv <- cor(y,results$cv.fit)^2
    	  cat("Original R-square = ",r2,"\n")
    	  cat(k,"Fold Cross-Validated R-square = ",r2cv,"\n")
    	  cat("Change = ",r2-r2cv,"\n")
    }
    shrinkage(lm)
    # 代码有点问题
    

    深层次分析——变量的相对重要性

    # 标准化(先将数据标准化后再做回归分析)
    zdata <- scale(data) 
    
    # 相对权重(对所有可能子模型添加一个预测变量引起的R平方平均增加量的近似值,结果显示各变量对R平方的解释比例)
    relweights <- function(fit,...){
    	  R <- cor(fit$model)
    	  nvar <- ncol(R)
    	  rxx <- R[2:nvar,2:nvar]
    	  rxy <- R[2:nvar,1]
    	  svd <- eigen(rxx)
    	  evec <- svd$vectors
    	  ev <- svd$values
    	  delta <- diag(sqrt(ev))
    	  lambda <- evec %*% delta %*% t(evec)
    	  lambdasq <- lambda^2
    	  beta <- solve(lambda) %*% rxy
    	  rsquare <- colSums(beta^2)
    	  rawwgt <- lambdasq %*% beta^2
    	  import <- (rawwgt / rsquare)*100
    	  import <- as.data.frame(import)
    	  row.names(import) <- names(fit$model[2:nvar])
    	  names(import) <- "Weights"
    	  import <- import[order(import),1,drop=FALSE]
    	  dotchart(import$Weights,labels=row.names(import),xlab="% of R-Square",pch=19,
    	           main="Relative Importance of Predictor Variables",
    	           sub=paste("Total R-Square= ", round(rsquare,digits=3)),...)
    	  return(import)
    }
    relweights(fit)
    
    展开全文
  • 使用R语言对数据分析主成分分析实现多元线性回归。包括源数据和代码。
  • R语言线性回归

    千次阅读 2017-06-01 21:05:46
    线性回归模型线性回归模型的计算

    线性回归模型

    线性回归模型的计算

    #lm()可以完成多元线性回归函数的估计,回归系统与回归方程的检验的工作
    #summary()函数,返回列表内容
    #X1表示体重,x2表示年龄,Y表示对应体重与年龄下的血压
    blood<-data.frame(
        X1=c(76.0, 91.5, 85.5, 82.5, 79.0, 80.5, 74.5, 
             79.0, 85.0, 76.5, 82.0, 95.0, 92.5),
        X2=c(50, 20, 20, 30, 30, 50, 60, 50, 40, 55, 
             40, 40, 20),
        Y= c(120, 141, 124, 126, 117, 125, 123, 125,
             132, 123, 132, 155, 147)
    )
    lm.sol<-lm(Y ~ 1+X1+X2, data=blood)#1表示常数项
    summary(lm.sol)
    call:
    lm(formula = Y ~ 1 + X1 + X2, data = blood)
    
    Residuals:
        Min      1Q  Median      3Q     Max 
    -4.0404 -1.0183  0.4640  0.6908  4.3274 
    
    Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
    (Intercept) -62.96336   16.99976  -3.704 0.004083 ** 
    X1            2.13656    0.17534  12.185 2.53e-07 ***
    X2            0.40022    0.08321   4.810 0.000713 ***
    ---
    Signif. codes:  0***0.001**0.01*0.05.0.1 ‘ ’ 1
    
    Residual standard error: 2.854 on 10 degrees of freedom
    Multiple R-squared:  0.9461,    Adjusted R-squared:  0.9354 
    F-statistic: 87.84 on 2 and 10 DF,  p-value: 4.531e-07

    总结

    1.call:为函数所使用的模型

    2.residuals:残差,列出了最小值,1/4,中位数,3/4,最大值。

    3.coefficients:系数,Estimate表示估计值(截距,x1,x2系数),3.1.std.error表示各个估计值的标准差
    3.2.t value表示t统计值
    3.3.Pr(>|t|)表示对应统计量的p值
    3.4.Signif. codes:显著性检测,p值在那个区间就用对应的星号表

    4.Residal standard error表示残差的标准差,自由度为n-p-1

    5.Multiple R-Squared 表示相关系数的平方,R^2

    6.Adjusted R-Squared 表示修正相关系数平方,这个值小于R^2,其目的是不要轻易做出自变量与因变量相关的判断。消除了R平方对自变量个数的依赖。

    7.F-statistic 表示F统计量,自由度为(p,n-p-1),p-value表示F统计量对应的p值

    Y = -62.96 +2.136*X1 + 0.4002*X2

    预测区间与置信区间

    #在R中使用predict()函数来计算y0的估计值,预测区间和E(y0)的置信区
    间。
    #置信区间估计(confidence interval estimate):利用估计的回归方程,对于自变量 x 的一个给定值 x0 ,求出因变量 y 的平均值的估计区间。
    
    #预测区间估计(prediction interval estimate):利用估计的回归方程,对于自变量 x 的一个给定值 x0 ,求出因变量 y 的一个个别值的估计区间。
    ## 一元回归估计
     x为含碳量,y为合金强度
    x<-c(0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 
         0.17, 0.18, 0.20, 0.21, 0.23)
    y<-c(42.0, 43.5, 45.0, 45.5, 45.0, 47.5, 49.0,
         53.0, 50.0, 55.0, 55.0, 60.0)
    lm.sol<-lm(y ~ 1+x)
    summary(lm.sol)
    
    ## 计算预测值, 并绘图
    new <- data.frame(x = seq(0.10, 0.24, by=0.01))#生成x的横坐标
    pp<-predict(lm.sol, new, interval="prediction")
    #3行15列(fit      lwr      upr)
    pc<-predict(lm.sol, new, interval="confidence")
    #3行15列(fit      lwr      upr)
    par(mai=c(0.8, 0.8, 0.2, 0.2))
    #这句程序?
    matplot(new$x, cbind(pp, pc[,-1]), type="l",
            xlab="X", ylab="Y", lty=c(1,5,5,2,2), 
            col=c("blue", "red", "red", "brown", "brown"), 
            lwd=2)
    #new$x为方向上的15个值,pp与pc的fit值相同,故去除相同一列,cbine表示y方向上的值
    points(x,y, cex=1.4, pch=21, col="red", bg="orange")
    #增加出x,y的点图
    legend(0.1, 63, 
           c("Points", "Fitted", "Prediction", "Confidence"), 
           pch=c(19, NA, NA, NA), lty=c(NA, 1,5,2), 
           col=c("orange", "blue", "red", "brow")
    )
    #增加标签

    回归与预测直线,置信曲线

    展开全文
  • 多元线性回归分析R语言

    万次阅读 多人点赞 2018-12-07 13:35:10
    ▼多元线性回归分析▼ 一、多元线性回归模型 设变量Y与X1,X2,……,Xp之间有线性关系   其中 , 和 是未知参数,p≥2,称上公式为多元线性回归模型。 二、参数估计 我们根据多元线性回归模型,认为误差...
  • 下面以物种多样性为例子展示了如何在R语言中进行相关分析和线性回归分析。 怎么做测试 相关和线性回归示例 Data = read.table(textConnection(Input),header=TRUE) 数据简单图 ...
  • 关于R语言多水平线性回归分析

    千次阅读 2020-09-11 20:24:05
    我再使用R语言分析嵌套模型的时候,遇到了一下问题,希望可以有大神指导一下。 数据是学校,班级,学生组成的嵌套模型。在模型构建的开始,我先确定模型中应该使用几水平。首先我先把三个水平均作为随机项放进模型中...
  • 作者 Lionel Hertzog译者 钱亦欣在一簇散点中拟合一条回归线(即线性回归)是数据分析的基本方法之一。有时,线性模型能很好地拟合数据,但在某些(很多)情形下,变量间的关系未必是线性的。这时,一般有三类方法解决这...
  • PAGE / NUMPAGES fire.txt 数据出自何晓群--应用回归分析 x y 3.4 26.2 1.8 17.8 4.6 31.3 2.3 23.1 3.1 27.5 5.5 36 0.7 14.1 3.0 22.3 2.6 19.6 4.3 31.3 2.1 24 1.1 17.3 6.1 43.2 4.8 36.4 3.8 26.1 #数据准备...
  • R语言多元线性回归分析

    万次阅读 2016-08-27 22:29:34
    我们这里采用逐步回归的方法来优化一下 y( 0.9 , 1.1 , 4.8 , 3.2 , 7.8 , 2.7 , 1.6 , 12.5 , 1.0 , 2.6 , 0.3 , 4.0 , 0.8 , 3.5 , 10.2 , 3.0 , 0.2 , 0.4 , 1.0 , 6.8 , 11.6 , 1.6 , 1.2 , 7.2 , 3.2 ) ...
  • #线性模型中有关函数#基本函数 a<-lm(模型公式,数据源) #anova(a)计算方差分析表#coef(a)提取模型系数#devinace(a)计算残差平方和#formula(a)提取模型公式#plot(a)绘制模型诊断图#predict...#多元线性回归分析 ...
  • R语言简单分析实际数据,进行简单线性回归和非线性回归
  • R语言与大数据编程实战》 学习笔记
  • 多元线性回归的适用条件: (1)自变量对应变量的变化具有显著影响 (2)自变量与应变量间的线性相关必须是真实的,而非形式上的 (3)自变量之间需有一定的互斥性 (4)应具有完整的统计数据 训练数据:csv格式,...
  • 分段回归( piecewise regression ),顾名思义,回归式是“分段”拟合的。... 分段回归最简单最常见的类型就是分段线性回归( piecewise linear regression ),即各分段内的局部回归均为线性回归。..
  • R语言与多元线性回归分析计算实例

    万次阅读 多人点赞 2019-11-16 12:37:33
    计算结果通过线性回归系数检验和回归方程检验,由此得到销售量与价格差与广告费之间的关系为: 模型的进一步分析 为进一步分析回归模型,我们画出y与x1和y与x2散点图。从散点图上可以看出,对于y与x1,用...

空空如也

空空如也

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

r语言线性回归分析