areg
################################################################
# 5 Graphing Commands 43
# 5.2 Adding Points, Lines, and Legends to Existing Plots . . . 45
################################################################
5.2 对已有图形添加点、线和图例说明
一旦你建立了一个图形，你可以添加点、线、文本或图例说明。安排这些元素的每个放置地方，R在定义图形面积的x轴和y轴术语条目时使其与之匹配，没有匹配定义的是图形窗口或设备条目。例如，如果你的图形中已有x值范围[0，100]的x轴，y轴的范围是[50,75]，那么你可在添加在(55,55)处添加一个点。
(a) points() 画一组或多组点。用pch对已有图形添加点。例如， points(P, Q, pch = ".", col = "forest green")，绘制每个(pi,qi)为很小的绿点。
(b) lines() 用线段连接指定的点。可以使用col和lty参量。例如，lines(X, Y, col = "blue", lty = "dotted")从每个设定的点(xi,yi)到下一个点画一条蓝色的虚线。另一方法，线也可对指定匹配的(x,y)采取命令输出。例如，density(Z)建立一个x向量和一个y向量，plot(density(Z))画出kernel密度函数图形。
(c) text() 在指定匹配的(x,y)处添加字符串。例如，text(5, 5, labels = "Key Point")在图形位置(5,5)处添加标签Key Point。你也可以用font可选项来选择字体，使用cex可选项来确定字体的相对于轴标签的大小，用col可选项选择颜色。可选项的完全列表用help(text)来访问。
(d) legend() 在与指定(x,y)相匹配的位置安排图例说明。键入demo(vertci) 来查看legend()示范。
<未完待续>
areg
################################################################
# 5 Graphing Commands 43
# 5.3 Saving Graphs to Files .............................. 45
################################################################
5.3 图形文件保存
默认情况下，R显示在你的屏幕有上有一个显示图形的窗口。把R图形保存为文件(包括在纸质文件中的格式，例如)，以你的图形命令开始：
> ps.options(family = c("Times"), pointsize = 12)
> postscript(file = "mygraph.eps", horizontal = FALSE, paper = "special", width = 6.25, height = 4)
此外ps.options()命令设置字体类型和输出文件的大小，postscript命令允许你指定文件名和其它几个附加的可选项。paper=special允许你指定包围附言区域的宽度和高度，用的是英寸(本例中，6.25英寸长，4英寸高)，horizontal=FALSE语句抑制R的默认图形方向。另一方法是，你可用pdf()而不是postscript()。如果你希望选择postscript的可选项为.pdf输出，你可以用pdf()的可选项来完成。例如：
> pdf(file = "mygraph.pdf", width = 6.25, height = 4, family = "Times", + pointsize = 12)
在每个图形结束，你应该关闭你的输出设备。命令dev.off()停止写入和在你的工作目录中保存.eps或.pdf文件。如果你忘记关闭这个文件，你将写写全部随后系列的图形到相同的文件，覆盖了先前的图形。你也可用dev.off()来关闭屏幕的图形窗口。
在同一个文件中写入多个图，你可用下面的可选项：
(a) 在同一个.pdf文件安排各图形在不同的页面中：
> pdf(file = "mygraph.pdf", width = 6.25, height = 4, family = "Times", pointsize = 12, onefile = TRUE)
(b) 在同一页中的多图，对一个.pdf或.eps文件初始化，接着(在任何图形命令前)键入：
par(mfrow = c(2, 4))
这将建立一个两行四列的图层，你的图形语句将从左向右，从第一行接着第二行来组装图层。
<未完待续>
areg
################################################################
# 5 Graphing Commands 43
# 5.4 Examples ..................................... 47
# 5.4.1 DescriptivePlots:Box-plots ....................... 47
################################################################
5.4 示例
5.4.1 描述性图形：盒形图(箱丝图)
图：收入与受教育的年限关系(略)
使用Zelig包中的rurnout样本数据集，返下面的命令可以生成该图。
> library(Zelig) # 载入Zelig程序包
> data (turnout)
# 载入样本数据，生成一个收入与受教育年限的盒形图。
> boxplot(income ~ educate, data = turnout, col = "grey", pch = ".",
main = "Income as a Function of Years of Education",
xlab = "Education in Years", ylab = "Income in \$10,000s")
################################################################
# 5 Graphing Commands 43
# 5.4 Examples ..................................... 47
# 5.4.2 Density Plots:A Histogram ....................... 48
################################################################
5.4.2 密度图：直方图
直方图以容易的方式来解析“感兴趣的统计量”的密度。
图：年收入频数的直方图
图形编码：
# 载入Zelig程序包
> library(Zelig)
# 载入样本数据集
> data(turnout)
# 调用main plot和可选项
> truehist(turnout$income, col = "wheat1", # Calls the main plot,
with xlab = "Annual Income in $10,000s", main = "Histogram of Income")
# 添加 kernel密度线
> lines(density(turnout$income))
<未完待续>
areg
################################################################
# 5 Graphing Commands 43
# 5.4 Examples ..................................... 47
# 5.4.3 Advanced Examples ........................... 49
################################################################
5.4.3 高级图形示例
上面的例子是简单示例，仅是对R制图潜力的一个表面游览。我们在Zelig 示例脚本中含有更高级的、特殊模型图形，有建立函数来自动生成这些图形，包括：
1. 三元图描述有三个观测值输出的分类因变量的概率预测。你可以有多元逻辑、秩逻辑或秩概率模型中使用这种制图(Katz and King,1999)。样本编码见12.17部分，在R提示符下键入demo(mlogit)来演示示例，参考12.17部分来给这个三元图形添加点。
图：1988 Mexican Presidential Election略
2. ROC Plots 概述二元因变量拟合数据的模型过程(logit, probit, 和 relogit)。ROC图形解析0和1数码来为每种可能的阈限值作出正确的预测，对连续概率(Y=1)可以认为是一个二歧预测。比较接近ROC曲线的是图形的右上角，比较好地拟合了模型规则(King and Zeng, 2002b)。样本编码见3部分，在R提示符下键入demo(roc)来执行这个示例。
3. 垂直的置信区间，几乎可以用于任何模型，虽然允许解释变理之一的值的改变超过给定值的范围，对任何“感兴趣的统计量”置信区间的模拟进行描述(King, Tomz and Wittenberg, 2000)。在R提示符下键入demo(vertci)来执行这个示例，键入help.zelig(plot.ci)来获取帮助手册页。
################################################################
#
# I Zelig用户指南(基础部分介绍完毕)
# II Zelig高级应用(从第6章开始，将稍晚点推出)
#
################################################################
<未完待续>
areg
###########################################################
# 6 R Objects 53
# 6.1 Scalar Values ................................... 53
# 6.2 Data Structures .................................. 54
# 6.2.1 Arrays ................................... 54
# 6.2.2 Lists.................................... 57
# 6.2.3 Data Frames ............................... 58
# 6.2.4 Identifying Objects and Data Structures . . . . 59
###########################################################
###########################################################
# 6 R Objects 53
# 6.1 ScalarValues ................................... 53
###########################################################
Part II Zelig高级应用
第6章 R 对象
在R中，对象可以是一种或多种类型，由标量值的类和保存标量值的数据结构的类组成。使用is()命令来确定一个对象是什么类型。如果你已经熟悉R对象，你可以跳过3.2.2部分的载入数据，或4.1部分描述的Zelig命令。
6.1 标量值
R中有几类标量值，他们组成更大的数据结构。R是以类为基础：某一个操作运算仅运行某一种类型的值，或某一种类型的数据结构。在此我们列出三种类型的标量值供你参考。
1. 数值 数值是大多数数目的默认类型。整数是数值类的子集，可以用作数值。对数值你可以执行方法任意类型或对数值进行逻辑运算，包括：
# 注意pi是内部常数，log()是自然对数函数。
> log(3 * 4 * (2 + pi))
[1] 4.122270
# 基本逻辑运算，包括>, <, >= (大于等于)，<=(小于等于)，== (等于)， !=(不等于)
> 2 > 3
[1] FALSE
# 高级逻辑运算符，包括& (和), && (if and only if控制和), | (或), 还有 || (控制或)
> 3 >= 2 && 100 == 1000/10
[1] TRUE
注意Ing(正无限)，-Inf(负无限)，NA(缺失值)，还有NaN(非数字缺失值)这些特殊数字值大多数学运算不能进行。(然而，逻辑运算可以进行)
2. 逻辑运算的结果是产生 TRUE 或 FALSE 逻辑值。转换逻辑值成数字值，使用as.integer()命令：
> as.integer(TRUE)
[1] 1
> as.integer(FALSE)
[1] 0
3. 字符值是文本字符串，例如：
> text <-"supercalafragilisticxpaladocious"
> text
[1] "supercalafragilisticxpaladocious"
在你的工作空间，把赋值符 <- 右侧的文本字符串赋予对象名。文本字符串主要用于数据框，文本描述的部分。R总是用双引号返回的字符串。
<未完待续>
areg
###########################################################
# 6 R Objects 53
# 6.2 DataStructures .................................. 54
# 6.2.1 Arrays ................................... 54
###########################################################
6.2 数据结构
6.2.1 数组
数组仅由一种类型的标量值组成(例如,一个字符串向量，或一个数字值矩阵)。最常见的是一维数组和二维数组，它们分别为向量和矩阵。
建立数组的方法途径：
1. 建立向量(或一维数组)的常用方法包括：
# 连接数字值成一个向量
> a <-c(3, 7, 9, 11)
# 连接字符串成一个向量
> a <-c("a", "b", "c")
# 建立一个从1到5的整数向量
> a <-1:5
# 建立一个由五个1组成的向量
> a <-rep(1, times = 5)
向量操作：
# 从向量a中提取第10个值
> a[10]
# 在向量a中，用3.14作为第5个值插入(即用新值替换原来的值)
> a[5] <-3.14
# 用2，4和7替换第5至第7的值
与较大的数组不同，不用先建立其它个正正确长度的向量，向量就可以扩展。因此：
> a <-c(4, 6, 8)
> a[5] <-9 # 向量a本来才有三个位置，现在在第5个位置插入9，那么自动在第4个位置插入NA。
2. 因子向量一种特殊类型的向量，允许用户在一个向量中建立j指示符，而不是使用j虚拟向量(如在Stata或SPSS)。R从已有的向量x用factor()命令来建立这种特殊类型的向量，按在x中的离散观测值把x分成不同水平层次。这些值可以是整数值或字符串。例如：
>x<-c(1, 1,1, 1, 1,2,2,2, 2, 9,9,9, 9)
> factor(x)
[1] 1111122 22 99 99
Levels: 1 2 9
默认情况，factor()命令建立无序因子，作为离散值处理而非秩序的水平。在向量中添加参量 ordered = TRUE 给因子排序。
> x <-c("like", "dislike", "hate", "like", "don’t know", "like", "dislike")
> factor(x, levels = c("hate", "dislike", "like", "don’t know"), ordered = TRUE)
[1] like dislike hate like don’t know like dislike
Levels: hate < dislike < like < don’t know
factor()命令根据在参量水平的可选项中的秩序来排定水平秩序。如果你省略水平命令，R将按它们在向量中的出现来确定秩序值。因此，省略水平参量所排序的水平如上面例子中将变成like < dislike < hate < don’t know。如果你在水平列表中省略一个或多个水平，R返回水平值NA作为缺失水平：
> factor(x, levels = c("hate", "dislike", "like"), ordered = TRUE)
[1] like dislike hate like <NA> like dislike
Levels: hate < dislike < like
在秩序逻辑和多元逻辑模型中(见4.2部分)，设置解释变量的值使用setx (见10部分)，数据框中使用因子向量制图(5.1部分)。
3. 从向量(一维数组)构建矩阵(或二维数组)
你可以两种方式来建立矩阵：
(a) 从向量建矩囝：使用命令matrix(vector, nrow = k, ncol = n) 从向量来建立一个 k×n 矩阵，在矩阵的列中从左到右填入。例如：
# 注意当给一个矩阵赋予向量，行或列都没有名字。
> matrix(c(1,2,3,4,5,6), nrow = 2, ncol = 3)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
(b) 从两个或多个长度为k的向量来建立矩阵：使用cbing()垂直组合n向量成格式 k×n 矩阵，或rbind()水平组合n向量成格式为 n×k 矩阵。例如：
> x <-c(11, 12, 13) # 建立由3个值组成的向量x
> y <-c(55, 33, 12) # 建立由3个值组成的向量y
# 建立一个 2x3 矩阵。注意，根据参量值传递给rbind()的秩序，1行名为x，2行名为y。
> rbind(x, y)
[,1] [,2] [,3]
x 11 12 13
y 55 33 12
# 建立一个矩阵，注意列名命名是根据传递值给cbind()的秩序确定。
> cbind(x, y)
x y
[1,] 11 55
[2,] 12 33
[3,] 13 12
R支持许多种函数命令，包括：det()，返回矩阵的行列式； t(), 矩阵转置；solve(), 矩阵反置；%*%, 两个矩阵相乘。此外，dim() 命令返回矩阵的维数。与向量一样，方括号从矩阵中提取指定值和用赋值机制 <- 替换值。例如：
> loo[,3] # 提取loo的第3列
> loo[1,] # 提取loo中的第1行
> loo[1,3] <-13 # 在第1行第3列插入值13
> loo[1,] <-c(2,2,3) # 替换loo中第一行的值
如果你在替换行或列时遇到问题，请检查向量的维数dims()与矩阵的维数 dims()中你想替换的位置是否匹配。
4. n维数组是一个相同维数的堆栈矩阵。例如，通过堆栈z矩阵的每x行和y列，你可以用(x, y, z)维来建立一个三维数组。
> a <-matrix(8, 2, 3) # 建立一个 2x3 矩阵，元素为8
> b <-matrix(9, 2, 3) # 建立一个 2x3 矩阵，元素为9
建立一个 2x3x2 数组，第1水平是[,,1]，组成元素为8，第2水平是[,,2]，组成元素为9。使用方括号来提取值。例如，[1, 2, 2]提取第2水平第一行第二个值。你也可以使用 <- 运算符来替换值。
> array(c(a, b), c(2, 3, 2))
, , 1
[,1] [,2] [,3]
[1,] 8 8 8
[2,] 8 8 8
, , 2
[,1] [,2] [,3]
[1,] 9 9 9
[2,] 9 9 9
> z<-array(c(a, b), c(2, 3, 2))
> z[1, 2, 2]
[1] 9
如果数组是一维向量或二维矩阵，R将以更加有效的方法来处理数组。对数组很有帮助的三个函数：
-is() 返回组成数组的标量值的类型，还有数组的特殊类型(向量, 矩阵, 或更普通的数组)
-dims() 返回数组的容量。
> dims(b)
[1] 33 5
此表示数组是一个二维的(矩阵)，有33行5列。
-单括号 [ ] 显示在数组中的指定值。你想抽出或替换时，用在 [ ]内用逗号显示指定值的索引。
> dims(a)
[1] 14
> a[10] # 抽出向量a中第10个值
> dims(b)
[1] 33 5
> b[1:12, ] # 抽出b中头12个值
> c[1, 2] # 抽出c中第一行第二列的值
> dims(d)
[1] 1000 4 5
> d[ , 3, 1] # 抽出1000个值的向量
<未完待续>
areg
###########################################################
# 6 R Objects 53
# 6.2 DataStructures .................................. 54
# 6.2.2 Lists.................................... 57
###########################################################
6.2.2 列表
与数组不同，数组仅含有一种类型的标量值，列表的数据结构灵活，含有的值其类型和数据结构可以不同。列表的这个灵活性以致可以列出其它列表。例如，列表输出可以含有coef，一个回归系统的向量；方差，方差与协方差矩阵；其它列表条目可以用字符串描述数据。使用
names()函数来查看列表中各元素名和提取一个元素名，用
> names(output)
[1] coefficients variance terms
> output$coefficients
列表中元素是没有名的，用双方括号 [[ ]] 来提取元素：
> L[[4]] # 从列表L中提取第4个元素
> L[[4]] <-b # 用矩阵b替换列表L中的第4个元素，矩阵b是上一节定义的。
与向量一样，列表灵活的数据结构可以扩展，不用先建立元素正确数目的其它列表：
> L <-list() # 建立一个空列表
> L$coefficients <-c(1, 4, 6, 8) # 插入一个向量给列表，并且的把列表中向量命名为coefficients
> L[[4]] <-c(1, 4, 6, 8) # 在列表的第4个位置插入向量。如果列表还没有4个元素，空元素将用NULL值表示。
另一方法，你可以容易地用在你工作空间已有对象来建立一个列表
> L <-list(coefficients = k, variance = v) # 在此k是一个向量，v是一个矩阵。
<未完待续>
areg
###########################################################
# 6 R Objects 53
# 6.2 DataStructures .................................. 54
# 6.2.3 DataFrames ............................... 58
###########################################################
6.2.3 数据框
数据框(数据集)是列表的特殊类型，它的每一个变量限制于有相同数目的观测值。一个数据框可以含有不类型变量(numeric, integer, logical, character, and factor)，只要每一变量有相同数目的观测值。
因此，数据框可以使用矩阵命令和列表命令来处理变量和观测值。
> dat[1:10,] # 提取第1至10观测值和全部相关变量
> dat[dat$grp == 1,] # 提取全部属于group1的全部观测值
> group <-dat$grp # 在工作空间保存变量grp为向量group，不是在数据框。
> var4 <-dat[[4]] # 在工作空间以var4保存第4个变量
数据和数据框更综合的介绍见3.2.2
<未完待续>
areg
###########################################################
# 6 R Objects 53
# 6.2 DataStructures .................................. 54
# 6.2.4 Identifying Objects and Data Structures . . . . 59
###########################################################
6.2.4 识别对象和数据结构
每个数据结构有几个属性来描述它。尽管这些属性对服用户正常情况下是看不见的(如，当键入对象名不在屏幕显示)，有几个帮助函数来显示特定属性：
-对于数组，dims()返回每一维的容量
-对于数组，is()返回标量值类型和数组的特殊类型(vector,matrix,array)。对于更复杂数据结构，is()对那个对象返回默认方法(classes) 。
-对于列表和数据框，names()返回变量名，str()返回变量名和每一成份的简要描述。
对于几乎所有的数据类型，你可以用summary()来获得概述统计。
<未完待续，稍晚一些时间帧出第7章，编程设计语句>
areg
########################################################
# 7 Programming Statements 60
# 7.1 Functions ..................................... 60
# 7.2 If-Statements ................................... 60
# 7.3 For-Loops ..................................... 61
########################################################
########################################################
# 7 Programming Statements 60
# 7.1 Functions ..................................... 60
########################################################
第7章程 程序设计语句
这章介绍主要的程序命令。包括函数、条件语句if-else、for循环和管理统计模型输入的专门步骤。
7.1 函数
函数有内部函数或用户自定义的一组封装命令，用它来运行任意数目的参量。用函数语句为开端和 < 操作符来赋予函数给你工作空间的对象。
在你的工作空间，你可以用函数来对不同对象执行相同的操作。例如：
check <-function(p, q) {
result <- (p-q)/q
result
}
这是一个参量为p和q的简单函数，它计算向量p的第i个元素和向量q的第i个元素的差与q中第i个元素之比，返回一个向量结果。例如，check(p = 10, q = 2)返回4。你可以略于描述，只要你保持参量以正确的秩序：check(10, 2)也返回4。你也可以用其它对象输入给这个函数，如呆again=10，really=2，接着check(p = again, q = really)和check(again, really)返回值也是4。
因为函数是以一个集合来执行命令，你确信在你函数中的每个命令的正常运行，在R提示符下测试函数的每一行。
<未完待续>
areg
########################################################
# 7 Programming Statements 60
# 7.2 If-Statements ................................... 60
########################################################
7.2 If语句
使用if(和可选项，else)来控制R函数的流。例如，让x和y都是标量数字值：
# 如果逻辑语句在()中为真，那么把x改变为NA(缺失值)。如果if语句为假，else语告诉R该如何做。
if (x == y) {
x <-NA
}
else {
x <-x^2
}
作为一个函数，使用 { and } 用if和else语句来定义每个相互联系的一组命令。(如果你在函数内部包含有if语句，你可以用多组嵌套的花括号，即大括号)
<未完待续>
areg
########################################################
# 7 Programming Statements 60
# 7.3 For-Loops ..................................... 61
########################################################
7.3 For循环
用for来执行重复(循环)运算。使用矩阵或向量命令时避免循环通常更快和更巧妙，但是循环在某此时候对赋值是必需的。如果你使用一个循环来给一个数据结构赋值，你必需首先初始化一个空数据结构以来保存你赋予的值。
选择一个数据结构与你的循环将生成的输出类型一致。如果你的循环生成标量，把标量存在一个向量中(向量的中第i个值对应于第i次循环的运行)。如果你的循环生成的输出向量，在一个矩阵中以行(或列)保存，第i行(或列)对应于循环的第i次重复。如果你的输出由矩阵组成，堆栈它们成一个数组。输出的列表(例如回归输出)或在每次重复中改变了输出的维数的列表。初始化这些数据结构，如下：
> x <-vector() # 任意长度的空向量
> x <-list() # 任意长度的空列表
命令 vector() 和 list() 建立一个任意长度的向量或列表，x[5] <-15 自动建立5元素的向量，开头四个值是空值(NA)。相反，命令
matrix() 和 array() 建立的数据结构对它们的最初维数是有严格限制的。
> x <-matrix(nrow = 5, ncol = 2) # 建立一个5行2列的矩阵
> x <-array(dim = c(5,2,3)) # 一个由3堆5行2列矩阵组成的三维数组
如果你试着在(100, 200, 20)位置赋予一个值给上面这两种结构，R将返回一个错误信息(“subscript is out of bounds”，即下标出界)。R不会自动扩展矩阵或数组的维数来接纳添加的值。
例 1: 用逻辑语句建立一个向量
# 初始化一个空数据结构，通过从1到10的每个值的循环，用 i+5 替换在x中满足条件的那个数，把命令用花括号封装。
x <-array()
for (i in 1:10) {
if (is.integer(i/2)) {
x[ i ] <-i + 5
}
}
你可以把 for()用在函数的内部或外部。
例 2: 用手建立一个虚拟变量
你也可以用循环来建立一个虚拟变量的矩阵添加到一个数据框。例如，对于每个州产生的固定效应，让我们假设你已有数据mydata，它含有y, x1, x2, x3, 和 state, state是50个唯一值的字符变量。在此有三种方式来建立虚拟变量：1) 用R的内部命令；2) 用一个循环；或3)用2来循环。
1. R将从不同值的单个变量来凭空生成一个虚拟变量。
> z.out <-zelig(y ~ x1 + x2 + x3 + as.factor(state),
data = mydata, model = "ls")
这个方法返回 k-1 显示为k个州。
2. 另一方法，你可以手动用循环来生成虚拟变量。有两种方式来完成，但都以同样的初始命令开头。用向量命令，先建立一个每个州的索引，初始化矩阵来保存虚拟变量：
idx <-sort(unique(mydata$state))
dummy <-matrix(NA, nrow = nrow(mydata), ncol = length(idx))
现在在两种方法中选择：
(a) 第一种方法的计算效率低，但对不习惯向量操作的用户更直观。第1循环用i作为在索引中让循环通过所有的行，第2循环用j来循环向量索引中全部50个值，它相应为列1虚拟在矩阵中通过50。
for (i in 1:nrow(mydata)) {
for (j in 1:length(idx)) {
if (mydata$state[i,j] == idx[j]) {
dummy[ i,j ] <-1
}
else {
dummy[ i,j ] <-0
}
}
}
接着添加虚拟变量的新矩阵给你的数据框：
names(dummy) <-idx
mydata <-cbind(mydata, dummy)
(b) 当你用向量操作变量更加舒服时，你可以用一个循环来替换上面的双循环操作。
for (j in 1:length(idx)) {
dummy[,j ] <-as.integer(mydata$state == idx[ j ])
}
单循环程序解析在索引中对应于向量mydata$state的每个元素。这建立 n个 TRUE/FALSE 观测值组成的向量，你可以用as.integer()把逻辑
值转换成1和0。赋予结果向理给虚拟中适合的列。组合虚拟矩阵给数据框，如上面来完成操作。
例 3：给子集加权回归
选择zelig()中的可选项划分数据框，接着自动循环指定模型通过每个划分。假定mydata是一个有变量y, x1, x2, x3, 和 state的数据框，
state是50个唯一值的因子变量。让我们说，你将执行一个加权回归，在此每个观测值被关于x1的标准误倒数加权，估计观测的州。换句话
说，首先我们需要为这50个州中每个估计模型，为每个州的state j =1,..., 50计算1 / se(x1j)，接着赋予这些权重给mydata中的每个观测值。
-用zelig()中的可选项分别为每个州估计一个模型：
z.out <-zelig(y ~ x1 + x2 + x3, by = "state", data = mydata, model = "ls")
现在z.out是50个回归输出的列表。
-提取变量x1对于州水平回归中的每个标准误。
# 初始化空数据结构，vcov()建立一个方差矩阵，由于我们有截距，第2对角线值对应于x1。
se <-array()
for (i in 1:50) {
se[ i ] <-sqrt(vcov(z.out[[ i ]])[2,2])
}
-建立加权向量
wts <-1 / se
这个向量wts有50个值，对应着在z.out州水平回归输入的50个集合。
-赋予权重的向量给每个观测值，我们需要每个观测值的州标记与州恰当匹配。比较简洁的方法，假定州是有限数从1到50。
# 初始化空变量
mydata$w <-NA
for (i in 1:50) {
mydata$w[mydata$state == i ] <-wts[ i ]
}
我们用 mydata$state 作为索引(在方括号内)来赋予值给 mydata$w 。因此，每当州观测值等于5，循环赋予向量wts中第5个值给mydata中的变量w。如果我们在mydata中有500个观测值，我们用这种方法来匹配500个观测值的每个给予恰当的wts。
如果州是字符串而不是整数，我们可以用稍微复杂点程序：
mydata$w <-NA
idx <-sort(unique(mydata$state))
for (i in 1:length(idx) {
mydata$w[mydata$state == idx[ i ]] <-wts[ i ]
}
-现在我们可以执行我们的加权回归：
z.wtd <-zelig(y ~ x1 + x2 + x3, weights = w, data = mydata, model = "ls")
<在循环中，避免与论坛斜体字符冲突，前面多加了个空格，未完待续>
areg
###########################################################
# 8 Writing New Models 65
# 8.1 Managing Statistical Model Inputs ....................... 66
# 8.1.1 Describe the Statistical Model ...................... 66
# 8.1.2 Single Response Variable Models: Normal Regression Model . . . . . 67
# 8.1.3 Multivariate models: Bivariate Normal example . . . . . 70
# 8.2 Easy Ways to Manage Matrices ......................... 72
# 8.2.1 TheI ntuitive Layout ........................... 73
# The Computationally-Efficient Layout . . . . 73
# 8.2.3 The Memory-Efficient Layout ...................... 74
# 8.2.4 Interchanging the Three Methods .................... 74
###########################################################
第8章 编写新模型
用Zelig，在R中编写一个新模型是很简单的。(如果你已有一个模型，见第9章，怎么把它包括到Zelig中来)。用户输入时合理化使用工具，写一样新模型并不需要大量的编程设计知识，但要让开发者的焦点在模型的数学方法上。一般来说，写一新统计分析程序或模型应该遵循的步骤是：
1. 写下数学公式的模型。定义你需要的参数，把参数分组成方便的向量或矩阵，无论何时何种可能性(这需要使你的编码清晰)。
2. 写下编码。
3. 检测编码(通常使用Monte Carlo数据，那个让你知道被估计的真实值)和确信它能如预期运行。
4. 写出一些文献，解释你模型和你模型中的函数。
在第1步和第2步间的某些地方，你将需要把输入的数据转成你用来写下模型的数学标记。而不是重复整个编码区块，使用函数来合理化用户将需要运行你模型的命令数目。
输入命令会变得很复杂，用较少的命令来执行更多的步骤。这些结构的输入实际内容很多。如果你的函数有复杂难懂的语名，它将很难以使用，难以解释，难以用文献说明。如果你的函数容易使用，有一个直观语句，那么它很容易解释和文献说明，它将使你的程序更容易接近用户。
<未完待续>
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)矩阵。
<未完待续>
