Fourth

 读取数据

data=read.csv("willett.csv")

建立模型Linear Mixed-Effects Models

查看数据

head(data)

##   id time cons covar   y
## 1  1    0    1   137 205
## 2  1    1    1   137 217
## 3  1    2    1   137 268
## 4  1    3    1   137 302
## 5  2    0    1   123 219
## 6  2    1    1   123 243

读取相应的package

library(GGally)

library(nlme)
library(MASS)

#计算数据的均值 最大值 最小值 中位数

summary(data)

##        id          time           cons       covar             y
##  Min.   : 1   Min.   :0.00   Min.   :1   Min.   : 81.0   Min.   : 76.0
##  1st Qu.: 9   1st Qu.:0.75   1st Qu.:1   1st Qu.:104.0   1st Qu.:171.5
##  Median :18   Median :1.50   Median :1   Median :115.0   Median :205.0
##  Mean   :18   Mean   :1.50   Mean   :1   Mean   :113.5   Mean   :204.8
##  3rd Qu.:27   3rd Qu.:2.25   3rd Qu.:1   3rd Qu.:123.0   3rd Qu.:231.2
##  Max.   :35   Max.   :3.00   Max.   :1   Max.   :137.0   Max.   :310.0

绘制数据的散点图

ggpairs(data)

## Warning in cor(x, y, method = corMethod, use = corUse): 标准差为零

## Warning in cor(x, y, method = corMethod, use = corUse): 标准差为零

## Warning in cor(x, y, method = corMethod, use = corUse): 标准差为零

## Warning in cor(x, y, method = corMethod, use = corUse): 标准差为零

a

混合效应回归模型

library(nlme)

nlme包的建模公式如下,左边为因变量。右边是变量的交互影响 然后声明随机影响 最后是数据 time | id表示在分区里拟合时间和id的随机变量

m1.nlme = lme(y  ~time  ,
random = ~ time | id,
data = data)

获取模型概要

summary(m1.nlme)

## Linear mixed-effects model fit by REML
##  Data: data
##        AIC      BIC    logLik
##   1278.823 1296.386 -633.4114
##
## Random effects:
##  Formula: ~time | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr
## (Intercept) 34.62341 (Intr)
## time        11.50653 -0.45
## Residual    12.62843
##
## Fixed effects: y ~ time
##                Value Std.Error  DF  t-value p-value
## (Intercept) 164.3743  6.118857 104 26.86356       0
## time         26.9600  2.166602 104 12.44345       0
##  Correlation:
##      (Intr)
## time -0.489
##
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max
## -2.27206204 -0.60035761  0.02631436  0.60730532  1.58439675
##
## Number of Observations: 140
## Number of Groups: 35

intervals()函数返回模型的95%相关系数和随机效应的置信区间

intervals(m1.nlme)

## Approximate 95% confidence intervals
##
##  Fixed effects:
##                 lower     est.     upper
## (Intercept) 152.24036 164.3743 176.50821
## time         22.66355  26.9600  31.25645
## attr(,"label")
## [1] "Fixed effects:"
##
##  Random Effects:
##   Level: id
##                            lower       est.      upper
## sd((Intercept))       26.6890363 34.6234055 44.9165789
## sd(time)               8.5444654 11.5065323 15.4954440
## cor((Intercept),time) -0.7005278 -0.4499471 -0.1005903
##
##  Within-group standard error:
##    lower     est.    upper
## 10.70064 12.62843 14.90354

极大似然测试了自变量的差异性

anova(m1.nlme)

##             numDF denDF   F-value p-value
## (Intercept)     1   104 1428.0506  <.0001
## time            1   104  154.8394  <.0001

由于p值小于0.05,因此可以拒绝原假设不同时间对应的因变量没有显著差异

随机效应图

从下面的图可以看到数据点无规则的分布,残差属于白噪声随机误差。模型拟合较好。

plot(ranef(m1.nlme))

aJPG

# Residuals plots

res_lme=residuals(m1.nlme)
plot(res_lme)

Third

本例使用systemfit包里的Kmenta数据集。

library(systemfit)
data( "Kmenta" )
eqDemand <- consump ~ price + income
eqSupply <- consump ~ price + farmPrice + trend
system <- list( demand = eqDemand, supply = eqSupply )
fitols <- systemfit( system, data=Kmenta )
summary(fitols )

##
## systemfit results
## method: OLS
##
##         N DF     SSR detRCov   OLS-R2 McElroy-R2
## system 40 33 155.883 4.43485 0.709298   0.557559
##
##         N DF     SSR     MSE    RMSE       R2   Adj R2
## demand 20 17 63.3317 3.72539 1.93013 0.763789 0.735999
## supply 20 16 92.5511 5.78444 2.40509 0.654807 0.590084
##
## The covariance matrix of the residuals
##         demand  supply
## demand 3.72539 4.13696
## supply 4.13696 5.78444
##
## The correlations of the residuals
##          demand   supply
## demand 1.000000 0.891179
## supply 0.891179 1.000000
##
##
## OLS estimates for 'demand' (equation 1)
## Model Formula: consump ~ price + income
##
##               Estimate Std. Error  t value   Pr(>|t|)
## (Intercept) 99.8954229  7.5193621 13.28509 2.0906e-10 ***
## price       -0.3162988  0.0906774 -3.48818  0.0028153 **
## income       0.3346356  0.0454218  7.36729 1.0999e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.930127 on 17 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 17
## SSR: 63.33165 MSE: 3.725391 Root MSE: 1.930127
## Multiple R-Squared: 0.763789 Adjusted R-Squared: 0.735999
##
##
## OLS estimates for 'supply' (equation 2)
## Model Formula: consump ~ price + farmPrice + trend
##
##               Estimate Std. Error t value   Pr(>|t|)
## (Intercept) 58.2754312 11.4629099 5.08383 0.00011056 ***
## price        0.1603666  0.0948839 1.69013 0.11038810
## farmPrice    0.2481333  0.0461879 5.37226 6.2274e-05 ***
## trend        0.2483023  0.0975178 2.54623 0.02156713 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.405087 on 16 degrees of freedom
## Number of observations: 20 Degrees of Freedom: 16
## SSR: 92.551058 MSE: 5.784441 Root MSE: 2.405087
## Multiple R-Squared: 0.654807 Adjusted R-Squared: 0.590084

可以看到,上述供需结构方程的拟合结果。两个方程的残差相关性达到0.891179。price和income对consump在第一个方程的pvalue非常小,所以影响十分显著,同样的对farmPrice,trend对consump的p第2个方程pvalue也很小,影响也十分显著,但price却影响不显著p值大于0.05.

Second

数据来源R里的内置数据集mtcars.

head(mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

multcomp::glht函数的示例:

library(multcomp)
mtcars$gear <- factor(mtcars$gear)
mtcars$cyl<- factor(mtcars$cyl)
model <- aov(mpg ~ cyl + gear, data = mtcars)
res <- multcomp::glht(model, linfct = mcp(gear = "Tukey"))
plot(print(confint(res)))

##
##   Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = mpg ~ cyl + gear, data = mtcars)
##
## Quantile = 2.4785
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
##            Estimate lwr     upr
## 4 - 3 == 0  1.3241  -3.4535  6.1017
## 5 - 3 == 0  1.5002  -3.0970  6.0974
## 5 - 4 == 0  0.1761  -4.5122  4.8644

a

本示例展示了汽车油耗mpg与两个影响因素:气缸数cyl和前进档位gear的关系。先利用aov函数建立双因素方差分析模型,然后利用multcomp::glht函数对gear的各个因子水平进行Tukey比较,gear的因子水平有3,4,5这三个水平,输出的结果显示gear的4档比3档耗油多的估计值为1.3241,95%的置信区间为[-3.4553 6.1035].同样的5档与3档,5档与4档之间差距的估计值与置信区间也可以很清楚的看到。最后的图形将上述结果直观的展示出来,可以看到这些置信区间都包含了0值,这也意味着前进档位各个因子水平之间对mpg的影响差别不是很显著。