areg
显示结果
1. 用summary(s.out)显示你模拟数量的概述。你可以如下指定有效位数:
> print(summary(s.out), digits = 2)
2. 另一方法,用plot(s.out)制出结果图。
3. 你也可以用names(s.out)来查看在这个对象中的名字和元素的描述,用.运算符来提取特定的结果。对于大多数模型,它们是:s.out$qi$pr(表示提取预测值),s.out$qi$ev (提取期望值),s.out$qi$fd (提取期望值中的首要差异)。对于逻辑、概率、多项式逻辑和那些有序逻辑、有序概率模型、"感兴趣的统计量" ,也包括s.out$qi$rr (风险率)。(# 原文:For the logit, probit, multinomial
logit, ordinal logit, and ordinal probit models, quantities of interest also include s.out$qi$rr (the risk ratio).)
<未完待续>
fx2006ms
跟着areg老师的讲解看了几段,明天继续看!
谢谢
areg
[quote]引用第21楼fx2006ms于2006-12-01 00:11发表的“”:
跟着areg老师的讲解看了几段,明天继续看!
谢谢[/quote]
看了起步阶段的部分,对R有了初步的了解,最好回到原作者的原文去跟原著学习,是比较好的方式。
其实R并不难学,贵在坚持!
areg
#######################################################
# 4 Statistical Commands 29
# 4.2 SupportedModels................................. 37
########################################################
4.2 Zelig支持的模型
我们在这里列出Zelig中全部可用的模型,按因变量的性质来组织进行预测、解释或描述。
# 在文本中,不少字符不能用,只好用汉字等代替。此部分统计术语太多,本人了解有限,敬请先或只参考原文来学习Zelig。
1. 无边界连续因变量的可取任意实数值,范围是(负无穷,正无穷)。大多数这些模型中设的是一个连续性因变量,而Bayesian因子分析多个连续性因变量。
(a) "ls" 线性最小二乘法(见12.16)计算最小残差平方和的系数。这是计算线性回归系数常用的方法,返回无偏差估计量β和方差(对指定模型设有条件)
(b) "normal" 正态(见12.20)模型对正态随机成分和有序线性成分计算极大似然估计。系数标识为ls,但是极大似然估计是对方差的一致性而不是偏离程度。
(c) "normal.bayse" 正态Bayesian回归模型(见12.21)与极大似然Gaussian回归相似,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
(d)"tobit" tobit回归模型(见12.30)检查左侧观测值的正态分布。
(e)"tobit.bayes" Bayesian tobit分布(见12.31)检查左侧和/或右侧观测值的正态分布。
(f)"factor.bayes" Bayesian因子分析模型(见12.7)估计多个连续型因变量作为潜在解释变量的函数。
2. 二歧因变量由两个离散值组成,通常是(0, 1)。
(a) "logit" Logistic回归(见12.13)指定Pr(Y=1)给一个由一组解释变量组成的线性函数(n 逆)逻辑转换。
(b) "relogit" 稀有事件逻辑回归选项(见12.29)估计,作为相同的逻辑模型,但它是对于稀有事件(当输出结果之一比其它的更常见)偏差的校正。它也是对校正抽样设计选择基础(对照)而进行提前校正的选项。
(c) "logit.bayes" 逻辑Bayesian回归(见12.14)与极大似然逻辑回归相信,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
(d) "probit" 概率回归(见12.37)指定Pr(Y=1)给一个由一组解释变量组成的线性函数(n 逆)CDF正态转换。
(e) "probit.bayes" Bayesian概率回归(见12.28)与极大似然概率回归相似,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
(f) "blogit" 二元变量逻辑模型(见12.1),模型中Pr(Yi1 = y1,Yi2 = y2)对应于(y1,y2) = (0, 0), (0, 1), (1, 0), (1, 1),根据二元正态分布密度。
(g) "bprobit" 二元变量概率模型(见12.2),模型中Pr(Yi1 = y1,Yi2 = y2)对应于(y1,y2) = (0, 0), (0, 1), (1, 0), (1, 1),根据二元正态分布密度。
(h) "irt1d" 一维项响应模型(见12.11),取值为多个二歧因变量和把模型中的值作为一个潜在(未观察到的)解释变量的函数。
(i) "irtkd" k维项响应模型(见12.12),取值为多个二歧因变量和把模型中的值作为k个潜在(未观察到的)解释变量的函数。
3. 序号常用于模型中排序和离散因变量。变量的输出值(如kill, punch, tap, bump)是有序的,但是任意两个连续分类间的距离不能精确确定。每个因变量可以认为是线性的,是一个连续的、未观察到的因变量观测值,执行的机制是返回一个有序的选择。
(a) "ologit" 秩序逻辑模型(见12.22)指定未观察变量的随机成分为一个标准逻辑分布。
(b) "oprobit" 秩序概率分布(见12.23)指定未观察变量的随机成分为一个标准正态分布。
(c) "oprobit.bayes" Bayesian秩序概率模型(见12.24)与秩序概率回归模型相似,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
(d) "factor.ord" Bayesian秩序因子分析(见12.9),做观测值模型,用有序因变量作为潜在解释变量的函数。
4. 多项因变量是无序的,为离散分类响应。例如,你可以在各个牌子的橙子汁间做单个人选择,或在一次选举中各个候选人间做单个选择。
(a) "mlogit" 多项逻辑模型(见12.17),根据多项随机成分和逻辑分类成分指定分类响应分布。
(b) "mlogit.bayes" Bayesian多项逻辑回归(见12.18)与极大似然多项逻辑回归相似,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
5.计数因变量是非负整数值,如总统否决的次数或撞到光子探测器的光子数目。
(a) "poisson" 泊松模型(见12.25)在给定观测期限的解释变量的指数函数,指定事件的期望值。泊松随机组分的特性就是λ = E(Yi|Xi) = V(Yi|Xi)
(b) "poisson.bayes" Bayesian Poisson回归(见12.26)与极大似然泊松回归相似,但是它进行的是经由来自exact posterior或priors进行的有效小样本推断。
(c) "negbin" 负二项式模型(见12.19)与泊松分布有相同的组分,但是允许事件计数来结束离散,如V(Yi|Xi) > E(Yi|Xi)
6. 有界连续因变量它的连续性限定有一定范围内,通常是(0,正无穷)。此外某些模型(exponential, lognormal, and Weibull)要求其值大于某些检查点,以至于这些因变量有的具有全部观测值,而其它一些变量仅有部分观测值(检查点)。
(a) "gamma" Gama模型(见12.10) 要求是正数,连续因变量为全部观测值(没有检查点)。
(b) "exp" 指数模型(见12.6) 右检查点因变量假定风险函数是超时常数。对于某此变量,主题可能是脱离现实的假设,或多或
少很可能再长一点对于解释变量将失败。
(c) "weibull" Weibull模型(见12.32)右检查点通过包括添加尺度参数,放宽风险常数的假设,α:如果α>1,失败的风险增加更长主题的存活;如果α<1,失败的风险降低更长主题的存活。然而zelig()默认估计的α,你可以可选α为任何大于0的值。固定α=1的结果是变成一个指数模型。
(d) "lognorm" 逻辑正态模型(见12.15) 对右检查点,在因变量指定的风险函数不是单调的,在观测周期部分增加了风险,其它部分降低了风险。
7. 混合因变量包括的模型,取值为超过一个变量,而这些变变来自两个或多个上面所分的类。(他们可以是不一致的类型。)
(a) Bayesian混合因子分析模型与Bayesianl因子分析模型、有序因子分析模型相反,以两种因变量作为潜在解释变量的函数为模型。
8. 生态推断模型估计未观测内部单元值,用观测值的行和列边界给出一个列联表。
(a) ei.hier 分层ei模型(见12.4)用于2×2列联表层面的估计。
(b) ei.dynamic Quinn’s 动态 Bayesian ei模型(见12.3),用临时依赖于全表来对2×2列联表估计一个动态Bayesian模型。
(c) ei.RxC R×C ei模型(见12.5)为超过2行或2列的列联表估计一个多项控制ei模型。
<未完待续>
areg
#######################################################
# 4 Statistical Commands 29
# 4.3 ReplicationProcedures .............................. 40
# 4.3.1 SavingReplicationMaterials....................... 40
# 4.3.2 ReplicatingAnalyses ........................... 41
########################################################
4.3 重复操作
大部分任何统计分析证明在给定同样数据下你的操作,任何人都能重复你的结果。此外,许多杂志要求引用和"数据集备份"的传播以让其它人能重复你的结果(see King, 1995)。如果你想为你自己的记录建立材料备份,或提供数据给其它人或出版物作为参考,Zelig使这个过程非常容易。
4.3.1 保存材料备份
让mydata作为你最后一个数据集,z.out是你的zlig()输出,s.out是sim()输出。把全部这些保存在一个文件中,键入:
> save(mydata, z.out, s.out, file = "replication.RData")
这就在你的工作目录下建立一个文件replication.RDdata。你可以用zip或gzip工具压缩这个文件。如果你运行几个具体结果,所有这些估计都可以保存为一个.RData文件。甚至你仅从这些模型中建立"感兴趣的统计量",你仍然可以保存全部具体结果在一个文件中。例如:
> save(mydata, z.out1, z.out2, s.out, file = "replication.RData")
尽管.Rdata格式能包括数据集,也能包括输出对象,但是它不是保存大数据集最有效的方式。用未压缩格式,ASCII文本文件保存数据所占空间比用.RData少。(当压缩时,文本格式的数据仍然比.RData格式数据小)。因此,如果你的观测值超过10万个,你希望把数据集与Zelig输出对象分开保存。为此,使用write.table()命令。例如,如果mydata是你工作空间的一个数据框,使用write.table(mydata, file = "mydata.tab") 来保存它为以tab分隔符的ASCII文本文件。你也可以指定其它分隔符;见help.zelig("write.table")中的可选项。
4.3.2 重复分析
如果数据集和分析全都保存在一个.RData文件中,并位于你的工作目录中,你可以简单健入:
> load("replication.RData") # 载入备份文件
> z.rep <-repl(z.out) # 仅重复模型
> s.rep <-repl(s.out) # 重复模型和"感兴趣的统计量"。
默认情况,repl()用相同的选项以来建立原始输出对象。因此,如果原始对象s.out使用bootstrapping进行245次模型,s.rep对象将相似地进行245次bootstraped模拟。此外,当重复"感兴趣的统计量"来再使用而不是重新恢复模拟参数,你可用prev可选项。健入help.zelig("repl")来查看repl()可选项的完整列表。
如果数据是保存为文本文件,用read.table()来载入数据,接着进行重复分析:
> dat <-read.table("mydata.tab", header = TRUE) # dat是在z.out中使用的相同的名
> load("replication.RData")
> z.rep <-repl(z.out)
> s.rep <-repl(s.out)
如果你载入数据有问题,请参考3.2.2部分。
最后,你可以用identical()命令来验证这个重复的相关输出的各步与原始zelig()输入是否一致(脚注4)。例如:
> identical(z.out$coef, z.rep$coef) # 检查系数如果参数是再模拟或再抽样,模拟"感兴趣的统计量" 将与原始数量有变化。如果你想用identical()来证实"感兴趣的统计量"是一致的,你可以用
# 再启用(保存)原始sim()输出中的模拟参数(和保存)
> s.rep <-repl(s.out, prev = s.out$par)
# 检查期望值的一致性。你可以为每个qi如下做。
> identical(s.out$qi$ev, s.rep$qi$ev)
脚注4 identical()命令检查数字值的小数位(16进制)的最大数目的一致性,也检查两个对象有相同的类(numeric, character, integer, logical, or factor)。更多信息参考help(identical)。
<未完待续,后面是第5章 图形命令>
areg
################################################################
# 5 Graphing Commands 43
# 5.1 DrawingPlots ................................... 43
# 5.2 Adding Points, Lines, and Legends to Existing Plots . . . 45
# 5.3 SavingGraphstoFiles .............................. 45
# 5.4 Examples ..................................... 47
# 5.4.1 DescriptivePlots:Box-plots ....................... 47
# 5.4.2 DensityPlots:AHistogram ....................... 48
# 5.4.3 AdvancedExamples ........................... 49
################################################################
第5章 图形命令
R和此类Zelig能生成特别美观的图形。已有许多内部函数,包括散点图、线图、直方图、条形图、饼图、趋势图、等高线图和各种各样的三维图形。如果你期望,你能练习一个绝对控制来生成所需的正确图形。Zelig为每个模型的一组模拟观测值来包含几默认图形。对于全屏查看这些图形,简单键入plot(s.out),此s.out是来自sim()的输出。依所所模型,plot()将返回不同的图形。
如果你想建立你自己的图形,这部分回顾了新建和保存二维图形的大多数基本步骤。R制图材料分成两步:
1.你必需调出输出设备(在5.3部分讨论),选择图形类型、画出图形的范围、画出轴线和给出制图的数据。在此步你也可以确定轴的标签、图的标题和为数据图形的颜色。第1步的概述在下面的5.1部分。
2。可选项,你可以对已有图形添加点、线、文本或图标符号。这些命令在5.2部分概述。
<未完待续>
areg
################################################################
# 5 Graphing Commands 43
# 5.1 DrawingPlots ................................... 43
################################################################
5.1 画图
最一般的制图命令是plot(),它自动识R对象的类型,你尝试制图和选择最好的图形类型。最普通的是用plot()返回图形,如下:
1. 如果X是一个长度为n的变量,plot(X),返回一个(xi,i) for i =1,...n,散点图。如果X没有排序,这将生成一个散乱的图形。用plot(sort(X))来按升序排列制出 (xi,i)各点的值。
2. 对于两个数值变量X和Y,长度都为n,plot(X,Y)绘出每个点(xi,yi) for i =1,...n的散点图。另一方法,如果Z是一个含有两个向量的对象,plot(Z)也建立散点图。
图形的可选参量包括:
(a) main 参量main给图形添加标题,xlab和ylab分别给x和y轴添加。例如:
plot(x, y, main = "My Lovely Plot", xlab = "Explanatory Variable",ylab = "Dependent Variable")
(b) type 参量type控制你要求的图形类型。默认值是plot(x, y, type = "p"),但是你可以在下面的类型中选择:
"p" 点
"l" 线 # 注意是"L"的小写,不是数字1
"b" 点和线
"c" 画线但不包括点
"h" 直方图
"s" 阶梯函数(阶梯图?)
"n" 一个空白图形区域(用轴指定)
(c) 如果你选择type = "p",默认R图形是小圆圈(小空心圆点)。你可以指定参量pch来改变点的类型。例如,plot(x, y, type="p", pch=19)建立实心圆的散点图。参量pch的其它可选项:
19 实心圆圈(碟子)
20 更小的实心圆圈
21 正方形
22 菱形
23 向上的三角形
24 向下的三角形
此外,你也可以用参量pch来指定自己的符号,例如pch = "*" 或 pch = "."。
(d) 如果你选择type = "l",默认是R图形画实线。用参量lty来设置线的类型。例如,plot(x, y, type = "l", lty = "dashed")画的图形是一条虚线。其它可选项是dotted, dotdash, longdash和twodash。
(e) 点、线或条形区颜色设置。例如, plot(x, y, type = "b", pch = 20, lty = "dotted", col = "violet")制成用虚线连接的小实心圆圈,两者都是紫罗兰色。(轴和标签保留黑色)。用colors()查看提供的完整颜色列表。
(f) xlim和ylim设置x轴和y轴的界限。例如,plot(x, y, xlim = c(0, 25), ylim = c(-15, 5)) 设定x轴的范围是[0,25],y轴的范围是[-15,5]。关于添加图形选项,参考help(par)。
<未完待续>
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步间的某些地方,你将需要把输入的数据转成你用来写下模型的数学标记。而不是重复整个编码区块,使用函数来合理化用户将需要运行你模型的命令数目。
输入命令会变得很复杂,用较少的命令来执行更多的步骤。这些结构的输入实际内容很多。如果你的函数有复杂难懂的语名,它将很难以使用,难以解释,难以用文献说明。如果你的函数容易使用,有一个直观语句,那么它很容易解释和文献说明,它将使你的程序更容易接近用户。
<未完待续>