R从入门到精通

areg

###########################################################

# 8 Writing New Models 65

# 8.1 Managing Statistical Model Inputs ....................... 66

###########################################################

8.1 管理统计分析分析模型输入

大多数统计分析模型需要一个解释变量的矩阵和一个因变量矩阵。而不是用户建立的那个矩阵自身，R有一个方便的用户界面来建立响应矩阵和凭空产生的解释变量。用户简单地以因变量~解释变量的格式来指定一个公式，开发者用下列函数来转换公式到适合的矩阵，使mydata成为一个数据框。

> formula <-y ~ x1 + x2 # 用户输入

# 给定上面的公式，程序设计都可以用下列标准命令。

> D <-model.frame(formula, data = mydata) # 明确指出不用子集和列表

> X <-model.matrix(formula, data = D) # 建立X矩阵

> Y <-model.response(D) # 建立Y矩阵

在此

-D是mydata的一个子集，它仅含有公式中指定的变量(y, x1, 和 x2)，并且在子集数据框不执行列表属性。

-X是一个含有1列的矩阵，该列来自D的解释变量x1和x2。

-Y是含有来自D的因变量。

依模型而定，Y可能是一个列向量、矩阵或其它数据结构。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.1 Managing Statistical Model Inputs ....................... 66

# 8.1.1 Describe the Statistical Model ...................... 66

###########################################################

8.1.1 描述统计分析模型

建立X矩阵后，大多数模型的下一步识别与参数对应的响应向量。对一个没有辅助参数的简单响应向量模型，标准R界面是很方面的：给定X，模型的参数是简单的β。

然面，很少有模型属于这种分类。例如，甚至正态回归，有两组参数β和σ^2。为使R的公式格式更灵活，Zelig添加了一组公式来让你描述对模型的输入(对于多组参数)。

当你写下一个统计分析模型后，标识你模型的参数。使参数智能化，第一步是写出你模型中describe.*()function 。如果你的模型称为mymodel，那么describe.mymodel() function没有设参量和返回下列信息的一个列表：

-category：描述因变量的字符串。见13.1部分，变量分类的当前列表。

-parameters：在你模型中设置的参数列表。对每个参数(e.g., theta)，你需要提供下面的信息：

(a)方程equations：给参数提供方程的整数数目。对于参数所用方程，例如，使用c(2,4)是2至4方程式。

(b)允许的标签tagsAllowed：指定逻辑值(TRUE/FALSE)，是否一个给定的参数允许是系统规定的参数。

(c)因变量depVar：指定逻辑值(TRUE/FALSE)是否一个参数要求一个相应的因变量。

(d)解释变量expVar：指定逻辑值(TRUE/FALSE)是否一个参数允许解释变量。

(见13.1示例，用describe.mymodel()来添加输出参量)

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.1 Managing Statistical Model Inputs ....................... 66

# 8.1.2 Single Response Variable Models: Normal Regression Model . . . . . 67

###########################################################

8.1.2 简单的响应变量模型：正态回归模型

Normal(yi | μi,σ^2)公式略

我们假定你试写一个正态回模型，含有标量随机成分的方差参数σ^2>0，系统成分E(Yi)= μi = xiβ。这暗示模型中设有两个参数，和下面的describe.normal.regression()函数。

describe.normal.regression <-function() { <br />
  category <-"continuous" <br />
  mu <-list(equations = 1, # 系统成分<br />
      tagsAllowed = FALSE,<br />
      depVar = TRUE,<br />
      expVar = TRUE)<br />
  sigma2 <-list(equations = 1, # 标量辅助参数<br />
         tagsAllowed = FALSE,<br />
         depVar = FALSE,<br />
         expVar = FALSE)<br />
  pars <-list(mu = mu, sigma2 = sigma2) <br />
  list(category = category, parameters = pars) <br />
}

查找对数似然：(公式略)

在R的编码，翻译成：

ll.normal <-function(par, X, Y, n, terms) { <br />
  beta <-parse.par(par, terms, eqn = "mu") # [1] <br />
  gamma <-parse.par(par, terms, eqn = "sigma2") # [2] <br />
  sigma2 <-exp(gamma) <br />
  -0.5 * (n * log(sigma2) + sum((Y -X %*% beta)^2 / sigma2)) <br />
}

在上面的注释[1]，我们用函数parse.par()抽出beta参数向量(它与系统成分与解释变量xi的μi相关)。不论那有多少协变量，函数parse.par()可以使用terms从par中抽出恰当的参数。在注释[2]，我也用也parse.par()来抽出对应σ^2的标量的辅助参数(转换之后)。

优化这个函数，简单键入：

out <-optim(start.val, ll.normal, control = list(fnscale = -1),

method = "BFGS", hessian = TRUE, X = X, Y = Y, terms = terms)

在此

start.val 是一个对par的开始值的向量。用set.start()来对每步中全部系统的参数、辅助的参数建立开始值。

ll.normal 是上面产生的对数似然函数。

"BFGS" 用拟牛顿法(quasi-Newton method)来指定无约束化优化(unconstrained optimization)。

hessian = TRUE 指示令R返回Hessian矩阵(来自你可以计算的方差-协方差矩阵)

X and Y 是在ll.normal()函数中，解释变量和因变量向量的矩阵。

terms 是来自model.frame()命令的元数据结构。

请参考R帮助中option()更多选项。

为使这个步骤概括，我们可以写一个函数来优化程序中关于用户指定数据框、公式和开始值可选项。

normal.regression <-function(formula, data, start.val = NULL, ...) { <br />
  fml <-parse.formula(formula, model = "normal.regression") # [1] <br />
  D <-model.frame(fml, data = data) <br />
  X <-model.matrix(fml, data = D) <br />
  Y <-model.response(D) <br />
  terms <-attr(D, "terms") <br />
  n <-nrow(X) <br />
  start.val <-set.start(start.val, terms)<br />
  res <-optim(start.val, ll.normal, method = "BFGS",<br />
      hessian = TRUE, control = list(fnscale = -1),<br />
      X=X,Y=Y, n=n, terms =terms, ...) #[2]<br />
  fit <-model.end(res, D) # [3] <br />
  fit$n <-n <br />
  class(fit) <-"normal" # [4] <br />
  fit <br />
}

下面的注释对应着上面括号中的数字：

1. parse.formula()命令搜索describe.normal.regression()函数，它改变用户指定的公式成下面的模式：

list(mu = formula,sigma = ~ 1) # formula是用户指定的。

2. ...省略号显示，当调用normal.regression()，如果用户键入任何添加参量，那么那些参量将进入optim()函数。

3. model.end()函数取优化的输出和从D中智能删除数据，将用setx()来运行建立一个对象。

4. 为你的模型输出选择一个类，以便你将能够为你的模型写入恰当的summary()、 param()和qi() 函数。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.1 Managing Statistical Model Inputs ....................... 66

# 8.1.3 Multivariate models: Bivariate Normal example . . . . . 70

###########################################################

8.1.3 多变量模型：二元变量正态的例子

大多数模型都有一个分类成分。对于n个观测值，分类成分变量是全部观测值 i=1,...,n。关于正态回归模型，系统成分是μi (σ^2不是作为相关变异函数的估计)然而，在某些情况下，你的模型可能超过一个分类成分。在二元变量概率情况中，我们有一个因变量 Yi=(Yi1,Yi2) 作为(0,0), (1,0), (0,1), or (1,1) for i =1,...,n的实测。与单响应概率模型相似，随机成分用两个潜在(未实测)的连续变量 (Yi . 1, Yi . 2) 来描述，它所遵循是二元变量正态分布：

## 公式见原文资料，此处略。

从上面公式所示，我们在两个设定参数μi =(μi1,μi2) 和ρ。这意指在describe.bivariate.probit()函数中：

describe.bivariate.probit <-function() {

category <-"dichotomous"

package <-list(name = "mvtnorm", version = "0.7")# 要求PKG的最低版本

mu <-list(equations = 2,tagsAllowed = TRUE,depVar = TRUE, expVar = TRUE),

# 分类成分要求有2个方程

rho <-list(equations = 1,tagsAllowed = FALSE, depVar = FALSE,expVar = TRUE),

# 分类成分可选项(默认时作为辅助参数的估计)

pars <-parameters(mu = mu, rho = rho)

list(category = category, package = package, parameters = pars)

}

由于用户可以选择不同的解释变量给参数μi1 和μi2 (有时还有ρ)，这个模型要求最少有两个公式。例如：

formulae <-list(mu1 = y1 ~ x1 + x2, mu2 =y2 ~x2+x3) # 用户输入

fml <-parse.formula(formulae, model = "bivariate.probit") # [1]

D <-model.frame(fml, data = mydata)

X <-model.matrix(fml, data = D)

Y <-model.response(D)

在注释[1]，parse.formula()去找出describe.bivariate.probit() 函数和分析相应的公式。如果ρ取协变量(变成一个分类成分而不是一

个潜在参数)，那就设了三个解释变量：

formulae <-list(mu1 = y1 ~ x1 + x2,

mu2 =y2 ~x2+x3, rho =~ x4 +x5)

透视程序设计，对于单个和多个方程模型有接近相同的工作框架。parse.formula()行改变从"list"到"multiple"中fml的类和确保model.frame()和model.matrix()采取适当的方法。D、X 和Y 类同于他们上面单个方程的副本：

-D 是mydata的子集，它含有变量y1, y2, x1, x2, 和对子集执行智能删除的x3;

-X 是一个对应于解释变量的矩阵，是下面讨论的三种格式之一(见8.2部分)

-Y 是一个 n×J 矩阵(本例此处J=2),以两列(y1, y2)对应于公式左侧变量的输出。

给定二元概率的概率密度如上所述，那它的似然表达式：见原文

# 以下几个公式及公式的解释说明见原文相应的程序编码见8.2.4。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.2 Easy Ways to Manage Matrices ......................... 72

###########################################################

8.2 管理矩阵的容易途径

大多数统计方法是从每个观测值i =1,...,n，解释变量xi与感兴趣的因变量yi之间有关。让β作为对应X中每列的一个参数，而又 n×k 矩阵有xi行。对于一个简单的方程模型，线性预测是

ηi = xiβ = β0 + β1xi1 + β2xi2 + … + βkxik

因此，η是一组 ηi，i =1,...,n，它通常以 n×1 矩阵来表示。对于两个方程模型，例如二元概率，线性预测变成一个列对应每个因变量(y1i,y2i)的矩阵:

ηi = (ηi1,ηi2) = (xi1β1, xi2β2)

η作为一个n×2矩阵，我们现在对于怎么建立线性预测有几种选择：

1. 一个直观层，就是堆栈解释变量的矩阵，提供显而易见的关于解释变量与相关系数之间关系表示；

2. 一个计算效率层，就是采用向量计算优势的特点；

3. 内存保存层，减少X和β矩阵容量大小。

在这部分使用简单工具描述，你可以选择对你的模型最好的矩阵管理方法。

此外，建立η的方式也影响着参数估计的方式。让我们来说，就是你想两个参数在不同方程中有相同影响。通过以某种方式设置X和β，你可以让用户通过参数来设置系统参数。上面连续的二元概率示例，让模型指定为：

formulae <-list(mu1 = y1 ~ x1 + x2 + tag(x3, "land"),

mu2 = y2 ~ x3 + tag(x4, "land"))

此处tag()是一个特殊函数，变量系统参数对所涉及的方程有相同的效果。因此，对方程中x3的系数mu1是相当于方程中x4的系数mu2，这种在两个方程中，作为land影响是相同的。为了考虑对各方程适合的系统参数，还有X和β的依据结构。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.2 Easy Ways to Manage Matrices ......................... 72

# 8.2.1 The Intuitive Layout ........................... 73

###########################################################

8.2.1 直观层次

# 公式帧出不便，同时在学习介绍中带来不少困难，请参考原文，有问题时，请指明是原文第某页，以便交流。

一个由X和向量β堆栈的矩阵可能是看起来最直觉的结构。让J=2为在二元概率模型中方程的数量，让vt为在两个方程中唯一协变量的总数。

选择model.matrix(..., shape = "stacked") 域为解释变量的(Jn × vt) = (2n × 6) 矩阵。再次，让x1为一个 n × 1 向量来表示变量X1, x2X3等等，那么X=（）# 略，见原文

相应地，β是一个由下面元素组成的向量。 # 组成表达式略，见原文。

此β0j是方程中j = {μ1,μ2}的一个截距术语。由于X是 (2n × 6)和β是(6 × 1)，线性预测结果η也堆栈成(2n × 1)矩阵。虽然难以处理(由于观测的索引是用i和2i来对应每个i =1,...,n,而非仅用i)，但是也容易看到，我们转两个方程为一个大的X矩阵和一个长向量β,这直接类似于我们熟悉的单个不同的方程。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.2 Easy Ways to Manage Matrices ......................... 72

# 8.2.2 The Computationally-Efficient Layout . . . . 73

###########################################################

8.2.2 计算效率层

# 请参考原文，文本帧出公式困难。

选择数组X和向量β可能是最有计算效率的配置：model.matrix(..., shape = "array")生成一个 n × kt × J 数组，此处的J是方程的总体数目，kt是跨越全部方程的参数总数。由于一些参数值在跨越方程时受约束，# 约束公式见原文。如果一个变量不在一定的方程中，它的观测当作0s的向量。以此观点，每个i =1,...,n xi 变成：

# 矩阵略，见原文。

从第一维堆栈这些xi矩阵中的每个，我们得到一个n × kt × J 维的X数组。相应地，β是一个多元素向量。# 元素表达式略。

用 (n × 6 × 2) 维的X数组和 (6 × 1) β 向量，对整个公式进行向量化：

eta <-apply(X, 3, ‘%*%’, beta)

因此线性观测变成一个(n×2)矩阵。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.2 Easy Ways to Manage Matrices ......................... 72

# 8.2.3 The Memory-Efficient Layout ...................... 74

###########################################################

8.2.3 内存效率层

# 请参考原文，文本帧出公式困难。

选择“压缩”X矩阵和β矩阵可能是最有效的内存配置：model.matrix(..., shape = "compact")(默认时)生成一个 n×v 矩阵，此处v是全部方程中唯一变量的数步(这个实例中是5)。让x1成为一个 n × 1 向量表示变量 x1, x2 x2等等，(# X向量和β矩阵略，见原文),βland参数使用两次来完成约束，空单元格的数目通过在β中完成约束而最小化，而不是在X中完成。更进一步，由于X是(n×5)和β是 (5×2)，Xβ=η是 n×2。

<未完待续>

areg

###########################################################

# 8 Writing New Models 65

# 8.2 Easy Ways to Manage Matrices ......................... 72

# 8.2.4 Interchanging the Three Methods .................... 74

###########################################################

8.2.4 三种方法互换

上面的连续型二元概率模型，我们仅需要修改几行编码放到这些不同的框架结构中使用其生效就行。使用的默认可选项(内存效率)，对数似然是：

bivariate.probit <-function(formula, data, start.val = NULL, ...) { <br />
      fml <-parse.formula(formula, model = "bivariate.probit") <br />
      D <-model.frame(fml, data = data) <br />
      X <-model.matrix(fml, data = D, eqn = c("mu1", "mu2")) # [1] <br />
      Xrho <-model.matrix(fml, data = D, eqn = "rho") <br />
      Y <-model.response(D) <br />
      terms <-attr(D, "terms") <br />
      start.val <-set.start(start.val, terms) <br />
      start.val <-put.start(start.val, 1, terms, eqn = "rho") <br />
          <br />
    log.lik <-function(par, X, Y, terms) {<br />
        Beta <-parse.par(par, terms, eqn = c("mu1", "mu2")) # [2]<br />
        gamma <-parse.par(par, terms, eqn = "rho")<br />
        rho <-(exp(Xrho %*% gamma) -1) / (1 + exp(Xrho %*% gamma))<br />
        mu<-X%*% Beta #[3]<br />
        llik <-0<br />
    for (i in 1:nrow(mu)){<br />
      Sigma <-matrix(c(1, rho[i,], rho[i,], 1), 2, 2) <br />
      if (Y[i,1]==1) <br />
        if (Y[i,2]==1) <br />
          llik <-llik + log(pmvnorm(lower = c(0, 0), upper = c(Inf, Inf), <br />
              mean = mu[i,], corr = Sigma)) <br />
        else <br />
          llik <-llik + log(pmvnorm(lower = c(0, -Inf), upper = c(Inf, 0), <br />
          mean = mu[i,], corr = Sigma)) <br />
    else <br />
      if (Y[i,2]==1) <br />
        llik <-llik + log(pmvnorm(lower = c(-Inf, 0), upper = c(0, Inf), <br />
             mean = mu[i,], corr = Sigma)) <br />
      else <br />
        llik <-llik + log(pmvnorm(lower = c(-Inf, -Inf), upper = c(0, 0), <br />
             mean = mu[i,], corr = Sigma)) <br />
      } <br />
<br />
    return(llik)<br />
    }<br />
    res <-optim(start.val, log.lik, method = "BFGS",<br />
         hessian = TRUE, control = list(fnscale = -1), <br />
         X = X, Y = Y, terms = terms, ...)<br />
    fit <-model.end(res, D)<br />
    class(fit) <-"bivariate.probit"<br />
    fit<br />
}

如果你发现默认的方法(内存效率)对执行你的模型并不是最好的方法，你可以使用直观可选项或计算效率可选项，如下改变几行编码就行：

-在上面注释[1]的直观选项处：

X <-model.matrix(fml, data = D, shape = "stacked", eqn = c("mu1", "mu2"))

和在注释[2]

Beta <-parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))

在注释[3]行保留与原始程序相同

-在上面的注释[1]行，计算效率可选项用下面的替换：

X <-model.matrix(fml, data = D, shape = "array", eqn = c("mu1", "mu2"))

注释[2]为

Beta <-parse.par(par, terms, shape = "vector", eqn = c("mu1", "mu2"))

在注释[3]为

mu <-apply(X, 3, ’%*%’, Beta)

甚至你的优化调用C或FORTRAN程序，在你传递这些数据结构给你优化程序之前，你可以用model.matrix() 和 parse.par() 组合来进行数据结构的设置以满足你需要获得的线性预测(或你的等效的模型)

<未完待续，接下来是第9章>

abraham_du

LZ真是辛苦了,谢谢!

areg

#######################################################

# 9 Adding Models and Methods to Zelig 77

# 9.1 Making the Model Compatible with Zelig. . . . . . . . 78

# 9.2 Getting Ready for the GUI............................ 84

# 9.3 Formatting Reference Manual Pages....................... 84

#######################################################

第9章给Zelig添加模型和方法

Zelig是高度模块化的。你可以对Zelig添加方法，如果你想，发行你的程序设计如同发行一个独立的程序包。使你的程序包与Zelig兼容，你申明你的程序包和帮助它完成一个推广使用。

本章假定你的模型是以函数来编写，采用的是用户自定义的公式和数据集(见第8章)，返回一个包括(最低限度是其中一个)估计参数和描述适合用于该模型的术语的输出列表。对于输出列表，你将选择一个类(S3或S4类)，对泛型函数提供恰当的方法，例如summary(), print(), coef() 和 vcov()。

给Zelig添加新模型，你需要提供6个R函数，如图形9.1那样说明(# 略，见原文)。让mymodel以类方法为myclass的新模型。

这些函数如下：

1. zelig2mymodel()翻译zelig()参量为mymodel()的参量。

2. mymodel()估计你的统计分析的程序或操作

3. param.myclass()模拟你模型的参数。另一可选方法，如果你模型的参数是一个对应于观测值为变量-协变量矩阵的向量组成，

你可以写两个简单的函数来置换param.myclall()：

(a) coef.myclass()从你的模型输出中提取系数

(b) vcov.myclass()从你的模型中提取变量-协变量矩阵

4. qi.myclass()计算期望值、模拟预测值和其它你模型中的“感兴趣的统计量”(仅应用于所取解释变量的模型)。

5. plot.zelig.mymodel() 对你模型的“感兴趣的统计量”的模拟图。

6. 写出模型的参考手册页(见9.3部分)

7. 一个描述对你的模型输入的函数(describe.mymodel())，对于使用提供一个图形用户界面(见13.1部分)

8. 一个可选的mymodel.R的样本脚本演示，你的参考手册页的示例中含有模型的注释码。

<未完待续>

areg

#######################################################

# 9 Adding Models and Methods to Zelig 77

# 9.1 Making the Model Compatible with Zelig. . . . . . . . 78

#######################################################

9.1 使模型与Zelig兼容

你可以开发一个模型，编写适合模型的函数，在Zelig框架内检测它，无需Zelig团队的介入。(当然，我们很高兴来回答对于改进的任何问题或建议)。

Zelig的模块性依赖于R的两个程序设计约定：

1. 封装，它从R函数传递参量给其它R函数或外部函数的调用(如C, C++, 或 Fortran函数)

2. 类方法，告诉泛型函数怎么处理给定类的对象。

为R泛型函数指定方法采取的一般格式： method.class()，method是执行的类属过程的名，class是对象的类。你可以定义，例如，以summary.contrib()来概述你的模型输出。注意S4类，泛型函数名不再有method.class()，只要用户调用method()就可以。

使用zelig()执行操作

Zelig对加入新模型的应用有独特方法，让模型编写者不用对任何zelig()函数修改，就可在Zelig的框架中来检测他们的模型。

使用zelig2contrib()封装函数(此处contrib是你的新模型名)，zelig2contrib()重新定义对zelig()的输入来执行对你的函数contrib()所需要的输入。例如，如果你键入

zelig(..., model = "normal.regression")

zelig()搜索zelig2normal.regression()封装程序，在任意环境中进行(在已经载入包或你的工作空间)。如果封装程序存在，那么zelig()就运行这个模型。

如果你有一个先前存在的模型，编写一个zelig2contrib()函数是很容易的。让我们假定你的模型是contrib()，取下列参量：formula、data、 weights和 start。相反，在zelig() 函数中，仅取formula、data、model和by参量。你可以使用...来从zelig()中传递添加的参量给zelig2contrib()，还有 <- NULL来略去你不需要的元素。从8.1.2部分连续性正态回归示例，让formula、model和data为对zelig()的输入，M是了集的数目，而 ... 是在zelig()调用中没有定义的添加参量，但是它们传递给 normal.regression()。

zelig2normal.regression <-function(formula, model, data, M, ...) { <br />
    mf <-match.call(expand.dots = TRUE) # [1] <br />
    mf$model <-mf$M <-NULL # [2] <br />
    mf[[1]] <-as.name("normal.regression") # [3] <br />
    as.call(mf) # [4] <br />
}

上面括号内的数字对应下面的注释：

1. 在zelig2normal.regression()中建立一个参量列表的方式来建立一个调用(一个用来解析的表达式)，包括用

normal.regression()提取的参量，但不是用zelig()。全部封装包必需取同样的标准化参量 (formula, model, data, and M),这些参量可能用在封装函数中来处理zelig()调入到normal.regression()调用中。添加参量给normal.regression()，如此，当start.val从zelig()中

以隐蔽方式传递给 ... 运算符。

2. 从调用mf对象中删除额外信息。在这个封装中，model和M没有被使用。在其它模型中，它们被用来进行处理更多的调用，也是包括在全部封装程序的标准输出中。

3. 对在normal.regression()中用来解析的函数名的调用(现在当前的是zelig2normal.regression)的第一个元素进行再赋值。

4. 返回调用给zelig()，它将为每个多重推理数据集或简单数据所解析调用，每一个子集用by来定义。

如果你使用S4类方法来表示你的模型，比如说mymodel，在zelig.default()内，Zelig的内部函数，建立.ZeligS4()，在两个添加槽(slots)的全局环境中，自动建立一个新S4类调用ZeligS4mymodel。这些包括保存模型名的zelig，如果save.data=TRUE和其它项为空，那还包括保存数据框的zelig.data。这些名字取自原始调用。新的输出继承原始的类mymodel，因此全部的泛型函数与mymodel联系在一起来运行。如果你想查看一个示例，查看该模型的应用，使用VGAM程序包，如多元概率。

使用setx()来执行操作

在setx()情况下，大多数模型将用setx.default()，它依次依敕于普通R函数model.matrix()。对这种程序的运行，你的输出列表必需包括：

-terms 用model.frame()建立或手动建立；

-formula 用户输入的公式对象；

-xlevels 它定义了解释变量的各层；

-contrasts 一个可选元素，它用在解释变量中定义因子变量的的类型。更多信息见help(contrasts)。

如果你的模型输出不能用setx.default()来执行，你必需编定你自己的setx.contrib()函数。例如，适合于多重推理的数据集有来自类方法MI的zelig()输出。特殊的setx.MI()在zelig()输出对象之前预封装程序，传递恰当的参量来setx.default()。

与sim()的兼容性

模拟“感兴趣的统计量”是解释模型结果的一个整体部分。使用这项功能来内嵌入Zelig sim()程序，你需要提供模拟参数的方式(调用param()函数)，对于计算或图示来自模拟参数(调用qi()函数)的“感兴趣的统计量”。

模拟参数

是否你选择使用默认方法或为模拟参数写入一个指定模型方法，这些函数要求同样的三个输入：

-object 估计模型或zelig()输出。

-num 模拟的数量

-bootstrap 选择TRUE或FALSE

来自param()的输出将是下列之一

-如果bootstrap=FALSE(默认)，一个行矩阵对应模拟和列对应模型参数。任意辅助参数将包括在输出矩阵中。

-如果bootstrap=TRUE，含有全部模型参数的向量，包括辅助参数。

有两种方式来模拟参数：

1. 使用param.default()函数来提取从模型中提取参数，如果bootstrapping没有选择，模拟系数使用渐进的正态近似。

(a) coef(): 提取系数。从上面的连续正态回归样本，恰当的coef.normal()函数是容易的：

coef.normal <-function(object)

object$coefficients

(b) vcov(): 提取变量-协变量矩阵。再次从上面的连续性泊松样本：

vcov.normal <-function(object)

object$variance

2. 另一方法，你可以编写你自己的param.contrib()函数。当下列情况中这是恰当的：

(a) 你的模型有辅助参数，如在正态分布情况下有σ。

(b) 你的模型对系数或变量-协变量矩阵执行某些类型的校正系数，而它不能在coef.contrib()或vcov.contrib()函数中执行。

(c) 你的模型对对数似然不依赖于渐近近似。对于Bayesian Markov-chain monte carlo模型，例如，param.contrib()函数

(param.MCMCzelig()在这种情况下)仅仅提取模型参数来模拟适合模型的函数。

连续性正态样本

param.normal <-function(object, num = NULL, bootstrap = FALSE,terms = NULL) {<br />
    if (!bootstrap) { <br />
        par <-mvrnorm(num, mu = coef(object), Sigma = vcov(object)) <br />
        Beta <-parse.par(par, terms = terms, eqn = "mu") <br />
        sigma2 <-exp(parse.par(par, terms = terms, eqn = "sigma2")) <br />
        res <-cbind(Beta, sigma2)<br />
    }<br />
    else {<br />
        par <-coef(object)<br />
        Beta <-parse.par(par, terms = terms, eqn = "mu")<br />
        sigma2 <-exp(parse.par(par, terms = terms, eqn = "sigma2"))<br />
        res <-c(coef, sigma2)<br />
    }<br />
    res<br />
  }

计算“感兴趣的统计量”

对于计算来自模拟参数的“感兴趣的统计量”，所有模型要求指定模型方法。对归类的模型，恰当的qi()函数是qi.contrib()。在没有最小值时，这个函数将计算下列“感兴趣的统计量”：

-ev: 期望值，计算来自对期望值的分析解，作为函数的的系统成分和辅助参数。

-pr: 预测值，从预测值定义的分布提取。对你的函数，如果R没有内部的随机数生成器，你可以从均匀分布随机抽取，使用CDF反置方法来计算预测值。

-fd: 在期望值中的首要差异(均数差)，通过提取预期值来计算从给定x1的期望值来考虑指定的x。

-ate.ev: 使用期望值ev计算平均处理效应。简易方式，y -ev，对每个观测值全局模拟的平均。

-ate.pr: 使用预测值pr计算平均处理效应。简易方式，y -pr, 对每个观测值全局模拟的平均。

qi()函数要求的参量：

-object: zelig输出对象

-par: 模拟参数

-x: 解释变量的矩阵(用setx()建立)

-x1: 对于首要差异(均数差，也是用setx()建立)备择值的可选项矩阵。如果首要差异不合适你的模型，如果x1不是NULL,你应该放入一条warning() or stop()。

-y: 因变量的矩阵或向量的可选项(计算平均处理效应)。如果平均处理效应不适合于你的模型，如果在setx()步骤中的条件预测被选择了，你应当放入一条warning() or stop() 。

从上面的连续性正态回归样本，适合的qi.normal()函数如下：

qi.normal <-function(object, par, x, x1 = NULL, y = NULL) { <br />
      Beta <-parse.par(par, eqn = "mu") # [1] <br />
      sigma2 <-parse.par(par, eqn = "sigma2") # [2] <br />
      ev <-Beta %*% t(x) # [3a] <br />
      pr <-matrix(NA, ncol = ncol(ev), nrow = nrow(ev)) <br />
      for (i in 1:ncol(ev)) <br />
         pr[,i] <-rnorm(length(ev[,i]), mean = ev[,i], # [4] <br />
                sigma = sd(sigma2[i]))<br />
      qi <-list(ev = ev, pr = pr)<br />
      qi.name <-list(ev = "Expected Values: E(Y|X)",<br />
         pr = "Predicted Values: Y|X") <br />
      if (!is.null(x1)){<br />
         ev1 <-par %*% t(x1) # [3b]<br />
         qi$fd <-ev1 -ev<br />
         qi.name$fd <-"First Differences in Expected Values: E(Y|X1)-E(Y|X)"<br />
  } <br />
<br />
if (!is.null(y)) {<br />
      yvar <-matrix(rep(y, nrow(par)), nrow = nrow(par), byrow = TRUE)<br />
      tmp.ev <-yvar -qi$ev<br />
      tmp.pr <-yvar -qi$pr<br />
      qi$ate.ev <-matrix(apply(tmp.ev, 1, mean), nrow = nrow(par))<br />
      qi$ate.pr <-matrix(apply(tmp.pr, 1, mean), nrow = nrow(par))<br />
      qi.name$ate.ev <-"Average Treatment Effect: Y -EV"<br />
      qi.name$ate.pr <-"Average Treatment Effect: Y -PR"<br />
    } <br />
    list(qi=qi, qi.name=qi.name) <br />
  }

上面有五行有编码注释。以下面的四种方式来改变这五行，你可以编写qi()函数来几乎对任何模型都适合：

1. 用你的系统参数的下标名来提取任意系统参数(系统参数在describe.mymodel()定义)。

2. 用他们此处的名字下标名来提取任意辅助参数(在describe.mymodel()定义)

3. 使用反置连接函数和η=Xβ计算期望值。(对于正态模型，这是线性的)。在注释[3a]和[3b]，你需要在这两个位置进行改变。

4. 用函数从你的模型的随机成分中随机抽取来替换rnorm()。

<未完待续>

areg

#######################################################

# 9 Adding Models and Methods to Zelig 77

# 9.2 Getting Ready for the GUI............................ 84

#######################################################

9.2 为GUI作准备

Zelig能执行许多种图形用户界面(GUIs)。GUIs的运行是通过识别接受特定模型的演绎，这些选项对用户仅出现在一些种类的图形界面中。

因而，对于你的模型运行界面的秩序，你必需用术语来描述你的模型，以使GUI能够被理解。使用第8章的指导方针来编写模型，你的模型将会与(最低限度)Virtual Data Center GUI兼容。对于先前存在的模型，对于你的模型跟随在13.1部分的示例，你将需要建立一个describe.*() 函数。

<未完待续>

areg

#######################################################

# 9 Adding Models and Methods to Zelig 77

# 9.3 Formatting Reference Manual Pages....................... 84

#######################################################

9.3 参考手册页的格式

Zelig最主要优点之一就是对模型有完整的文献，与由函数而组织的R文献中面向对象的程序设计的差别。因此，我们要求Zelig的撰稿人提供相似的文献，包括传递给zelig()的句法和参量，还有模型的系统组分和随机组分、“感兴趣的统计量”、输出的值和更多信息(包括参考)。有几个方式提供这些文献：

-如果你已有一个用.Rd帮助格式的PKG文献，help.zelig()将自动搜索除了Zelig帮助文件外的R帮助文件。

-如果你已有一个静态的URLs在线HTML格式的PKG文献(象Zelig或MatchIt)，你需要提供一个PACKAGE.url.tab文件，它是一个两列表，第一列是函数名，而url在第二列。(纵使该文件的扩展是.url.tab,，该文件应该是一个tab-or space-delimited文本文件。)例如：

command http://gking.harvard.edu/zelig/docs/Main_Commands.html

model http://gking.harvard.edu/zelig/docs/Specific_Models.html

如果你想检测来查看是否你的 .url.tab能运行，简单地把它放到你的R library/Zelig/data/ directory中。(你不需要重新安装Zelig来检测你文件的.url.tab。)

-首选方法：你可以提供LATEX 2ε.tex文件。这种文献使用的书籍格式和从下列包中支持命令：graphicx, natbib, amsmath, amssymb, verbatim, epsf 和 html。因为使用\include{}模型页被组合成一体进入这种文献中，在呈送之前，你将确信你的文献编辑。请粘连到下列约定到你的模型页：

1. 全部数学公式将能用equation*、array和eqnarray*，或校准环境来排版。请避免显示方式。(以html看起来有点滑稽)

2. 全部命令或R对象将使用texttt环境。

3. 模型以较大文献的亚组开始，在模型页内的各部分是亚-亚组水平。

4. 对于格式的一致性，请避免使用描述性环境。

每个LATEX模型页将包括下列元素。让投稿人指定新模型。

临时帮助文件

\subsection{{\tt contrib}: Full Name for [type] Dependent Variables}

\label{contrib}

\subsubsection{Syntax}

\subsubsection{Examples}

\begin{enumerate}

\item First Example

\item Second Example

\end{enumerate}

\subsubsection{Model}

\begin{itemize}

\item The observation mechanism, if applicable.

\item The stochastic component.

\item The systematic component.

\end{itemize}

\subsubsection{Quantities of Interest}

\begin{itemize}

\item The expected value of your distribution, including the formula

for the expected value as a function of the systemic component and ancillary paramters.

\item The predicted value drawn from the distribution defined by the corresponding expected value.

\item The first difference in expected values, given when x1 is specified.

\item Other quantities of interest.

\end{itemize}

\subsubsection{Output Values}

\begin{itemize}

\item From the {\tt zelig()} output stored in {\tt z.out}, you may

extract:

\begin{itemize}

\item

\item

\end{itemize}

\item From {\tt summary(z.out)}, you may extract:

\begin{itemize}

\item

\item

\end{itemize}

\item From the {\tt sim()} output stored in {\tt s.out}:

\begin{itemize}

\item

\item

\end{itemize}

\end{itemize}

\subsubsection{Further Information}

\subsubsection{Contributors}

<入门部分介绍完毕，后面各章节是参考手册>

areg

通过多天以来的整理，终于把Zelig的入门部分介绍完毕，如果只是掌握此部分，当然达不到精通，之所以选用“R从入门到精通”这样一个很俗的名，主要是因为我们现在不少人对搜索信息的习惯而选的，实际在以科学的精神来说，根本不提“精通”，在资料介绍时更多提的是某某使用手册。

对于这样一个非常优秀的R学习包，由于本人也才初学，不当之处一定很多，望以此资料作为学习参考的朋友们，当你在学习中发现错误的或不恰当的地方，请你花几分钟时间留贴指出，我尽可能即日更正，同时也非常感谢你提出批评意见。

对于入门部分，我随时会复查，对于不恰当的地方也争取尽快修正。后续部分，不会象入门阶段的知识这样即时贴出，也不会系统地贴出；但是会根据自己学习使用中的体验或各位朋友的反馈信息，针对某个问题贴出探讨。

由于这是一个学习R的贴子，有问题时，欢迎你回贴，大家共同学习！如果没有问题，“说好话或称赞贴之类”就不要贴出，以减少后来者翻阅的时间，同时也便于学习。

如果你在学习和使用R的过程中，觉得某份资料或某个PKG对初学R的人很有帮助，我和大家会非常支持你系统地整理贴出。

最后让我们一起学习R，共同进步！

# 如果你所在的学校或单位，因网络限制，不能访问下列网址，得不到原文资料

http://gking.harvard.edu/zelig/

http://gking.harvard.edu/

访问Zelig的主页不方便，得不到Zelig.PDF文献，请你以短消息或邮件告诉我，当网快时(晚上11点后)，我会尽快给你寄出。

<后面进行的是学习交流的讨论，欢迎你的参与>

biogene

areg，你真是一个大好人！

cny0120

楼主太赞了！

wuguohui

太好了啊，这样学起来就快多了啊

zhhxu1969

Thanks!

yaya_cau

太猛了！楼主！拷贝下来都快绝望啊！