用模拟来理解混合效应模型之一：random intercept model

混合效应模型（Mixed-effect Model）在之前文章提到过，但感觉仍是雾里看花。此番又研究了一些资料，准备来做一个系列讲讲清楚，也算是自己学习的一个总结。

普通的线性回归只包含两项影响因素，即固定效应（fixed-effect）和噪声（noise）。噪声是我们模型中没有考虑的随机因素。而固定效应是那些可预测因素，而且能完整的划分总体。例如模型中的性别变量，我们清楚只有两种性别，而且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model？因为有些现实的复杂数据是普通线性回归是处理不了的。例如我们对一些人群进行重复测量，此时存在两种随机因素会影响模型，一种是对某个人重复测试而形成的随机噪声，另一种是因为人和人不同而形成的随机效应（random effect）。如果将一个人的测量数据看作一个组，随机因素就包括了组内随机因素（noise）和组间随机因素（random effect）。这种嵌套的随机因素结构违反了普通线性回归的假设条件。

奇异值分解和潜在语义分析

奇异值分解是一种矩阵分解技术，如果将列看作特征，SVD也能看作是一种降维技术。如果将数据先进行标准化之后，再进行SVD操作，其左奇异向量U就是主成分得分，右奇异向量V就是主成分负载。但SVD的优势恰恰在于不进行标准化，当我们处理高维稀疏数据（例如文本数据）时，标准化会使矩阵不再稀疏，影响计算效率。所以一般是用SVD直接对原始数据进行操作。得到U的第一列是平均值的意义，U的后几列才有区别各样本的功能。

在文本挖掘中，潜在语义分析（LSA）假设在词项和文档之间有一个潜在语义层，词汇和语义产生关联，文档也和语义产生关联。这样通过语义概念，可以将词项和文档放在一个语义空间中观察分析。进一步我们可以将语义所为词项的“主成分”放到数据挖掘中进行分析。LSA的一般过程是先构造词项-文档矩阵（其转置为文档-词项矩阵，视需要而定），再计算TFIDF值，然后用SVD来分解这个矩阵，得到的奇异向量就是潜在语义。U的第一列可以认为是词项在文档中的平均出现频率，V的第一列可以认为是文档中包含的词项频率，我们更有兴趣的是后几列。关于SVD和LSA有一份精彩的入门文档，下面的R代码就是对该文档的模仿，但是SVD得到的结果却和文档不同，我很奇怪为什么会这样。（使用的数据）

说英雄谁是英雄之光荣三国

本文的缘起是知乎的这一篇文章，老马这位三国历史专家让大家都知道了有这么一份三国志游戏数据的存在。这份数据很容易从网上找到，而本人只关心三国志四的各人物数据，为什么呢，因为只玩过三国志四嘛。游戏是多年前玩的，如今来玩玩数据吧。下面的问题主要围绕一个焦点，谁是三国志四中文武兼资第一人。

生存分析函数小结

生存分析(survival analysis)适合于处理时间－事件数据。例如中风病人从首次发病到两次复发，其中就涉及到时间和事件。此例中时间就是复发的时间间隔，事件就是是否复发。如果用普通的线性回归对复发时间进行分析，就需要去除那些没有复发的病人样本。如果用Logistic回归对是否复发进行分析，就没有用到时间这个因素。而生存分析同时考虑时间和事情这两个因素，效果会更好些。

在R语言中我们可以使用survival包进行生存分析，其中主要的函数功能罗列如下：

Surv：用于创建生存数据对象
survfit：创建KM生存曲线或是Cox调整生存曲线
survdiff：用于不同组的统计检验
coxph：构建COX回归模型
cox.zph：检验PH假设是否成立
survreg：构建参数模型

培训广告

本月中旬本人将在上海举行一次为期三天的R语言商业培训，有兴趣的朋友可以点这里查看。

电影爱好者的R函数

作为一个伪影迷，经常纠结一些电影该不该下，要不要看。毕竟吾生也有涯而片源无涯。还好可以去豆瓣一类的地方看看大家的评分择优录用。去豆瓣查分需要登录网站搜索再鼠标点点点，如果要查好几部电影就有点费事儿。其实可以用R写个函数，先抓取相应的网页，再筛选返回需要的分值。这样在R里头就可以批量查分了，恩，走起来。

library(RCurl)
library(XML)
movieScore <- function(x) {
    stopifnot(is.character(x))
    # 提交搜索豆瓣表单
    search <- getForm("http://movie.douban.com/subject_search", search_text = x)
    searchweb <- htmlParse(search)
    # 解析搜索结果页面
    resnodes <- getNodeSet(searchweb, "//div[@id='wrapper']//table[1]//a")
    if (is.null(resnodes)) 
        return(NULL) else resurl <- xmlGetAttr(resnodes[[1]], name = "href")
    # 得到影片页面后第二次解析
    resweb <- getURL(resurl, .encoding = "UTF-8")
    content <- htmlParse(resweb, encoding = "UTF-8")
    resnodes <- getNodeSet(content, "//div[@id='interest_sectl']//p[@class='rating_self clearfix']//strong")
    namenodes <- getNodeSet(content, "//div[@id='content']//h1//span")
    # 得到影片评分
    score <- xmlValue(resnodes[[1]])
    name <- xmlValue(namenodes[[1]])
    return(list(name = name, score = score))
}

看看天机这部大烂片多少分。

movieScore("天机")

## $name
## [1] "天机·富春山居图"
## 
## $score
## [1] "2.9"

抓网页比较慢，豆瓣为人民群众着想提供了API，我们也可以使用API来调取分数，函数也比较简单。

library(RCurl)
library(XML)
library(RJSONIO)
movieScoreapi <- function(x) {
    api <- "https://api.douban.com/v2/movie/search?q={"
    url <- paste(api, x, "}", sep = "")
    res <- getURL(url)
    reslist <- fromJSON(res)
    name <- reslist$subjects[[1]]$title
    score <- reslist$subjects[[1]]$rating$average
    return(list(name = name, score = score))
}
movieScoreapi("僵尸世界大战")

## $name
## [1] "僵尸世界大战"
## 
## $score
## [1] 7.5

有了这个查分函数，我们可以在R中批量查阅电影评分了。但是豆瓣对于频繁的访问会有限制，对于没有认证的API使用是每分钟10次，超过就会暂时封IP。对于网页抓取，肖楠在第六次R会议上有个很棒的演讲，有兴趣的同学可以去统计之都看看。

关于可视化的一个slide

前段时间团队内部讨论，做了一个关于数据可视化的slide，放在下面的链接中，供有兴趣的各位看看吧。ubuntu下chrome效果最佳，可尝试按下ESC看全局。

https://dl.dropboxusercontent.com/u/6547489/vis-ppt/index.html

适合初学者的R语言教程

learning RStudio

Beginning_R

Begining-R

实现可重复的统计slides

制作幻灯片是数据分析师的必备技能之一，优秀的slides对外可以忽悠住客户，对内可以震慑住领导。精良的slides要秀外慧中，有逻辑有内容，还要有外形有风格。达到这种标准并不容易。在过去，你可能使用office中的ppt来做传统意义的幻灯片，将图片和代码费力的copy到一张张的slides上去，然后到处找模板。但这种方式已经凹特了。

在HTML5的发展背景下，已经出现了大批以网页形式的slides框架，例如：Google IO 2012\HTML5 slides\HTML5 Rocks\Shower\Deck.js这个名单可以很长。而这里只是其中六种演示框架的介绍。在这些slides框架下，你只需要懂一点点web知识，将图片和数据嵌入到一个html模板中，就可以生成一个动态可交互的slides。

初学D3的感觉

简单来讲有如下三点感受：

• D3很强大
• D3并不容易学
• 学会D3并不等于学会了可视化

一山还有一山高，在翻过ggplot2这座山后，发现还有D3这座珠穆朗玛。ggplot2已经非常好了。能够实现The Grammar of Graphics的精义，有很丰富的对象和灵活的自由度，但这都还不够，因为ggplot2只能够生成静态的图形。如果你只需要写一篇分析报告，那么ggplot2是足够的。但如果你想让数据在网页上飞翔，就需要D3做为翅膀。

R连接MySQL数据库方法备忘

R语言连接数据库可以利用数据库的存贮能力和R的计算能力，起到取长补短的效果。之前我们也说过了如何在R中使用SQL，很多教材上也提到了连接MySQL的方法，但是在安装上还需要注意一些细节问题。以下的解决方法也是在网上放狗加撞墙实验得出的结论。

R语言玩转资产组合计算

以前都是用MATLAB或是EXCEL给学生讲资产组合的计算问题，实际R语言也可以做一样的事情。资产组合要解决的问题并不复杂，即给定一个可选的资产集合，要从中选择出一个最优组合，使其收益率较大，而风险较小。如果可选资产只有两个的话，问题非常简单，可以通过求导最优化的方式得到结果，即有效资产前沿的解析解。如果资产集合超过两个那么问题要复杂一些。首先需要给定一个要求的期望收益率，再求出在此约束下的方差最小组合，这可以通过二次规划解出。这样得到有效前沿上的一个点。然后更改期望收益率，又可以得到另一个最优解，如此反复即可得到多资产组合的有效前沿。

下面我们先用R语言中quadprog包的二次规划求解函数来计算一个例子，再用fPortfolio包中的函数来完成同样的任务。例子中使用的数据来自于fPortfolio包中的SMALLCAP.RET数据集，使用了其中的四个资产来示例。

灰色模型的R代码

最近帮朋友写了一个灰色模型GM(1,1)的R实现，参考网上现有的matlab代码，比较容易就可以弄出来。下面是具体过程，主函数是GM()，建立的模型是一个S3类，搭配两个自定义的泛型函数print和plot可以得到结果输出和图形。

# 本代码用来计算GM(1,1)灰度模型
# 输入：x 原始数据，adv为外推预测步长
# 输出：actual 原始数据，fit 拟合数据,degree 拟合精度，
#       C 后验差比值，P 小概率误差，predict 外推预测值
 
 
# 主函数
GM <- function(x,adv=0) {
    x0 <- x
    k = length(x0)
    # AGO
    x1 = cumsum(x0)
    # construct matrix B & Y
    B = cbind(-0.5*(x1[-k]+x1[-1]),rep(1,times=k-1))
    Y = x0[-1]
    # compute BtB...
    BtB = t(B)%*%B
    BtB.inv = solve(BtB)
    BtY = t(B)%*%Y
 
    # estimate
    alpha = BtB.inv%*%BtY
 
    # 建立预测模型
 
    predict <- function(k) {
        y = (x0[1] - alpha[2]/alpha[1])*exp(-alpha[1]*k)+alpha[2]/alpha[1]
        return(y)
    }
    pre <- sapply(X=0:(k-1),FUN=predict)
    prediff <- c(pre[1],diff(pre))
    # 计算残差
    error <- round(abs(prediff-x0),digits=6)
    emax <- max(error)
    emin <- min(error)
    # 模型评价
    incidence <- function(x) {
         return((emin + 0.5*emax)/(x+0.5*emax))
    }
    degree <- mean(sapply(error,incidence))
 
    s1 <- sqrt(sum((x0-mean(x0))^2)/5)
    s2 <- sqrt(sum((error-mean(error))^2)/5)
 
    C <- s2/s1
 
    e <- abs(error-mean(error))
    p <- length(e<0.6745)/length(e)
 
    result <- list(actual = x0,
                   fit = prediff,
                   degree = degree,
                   C = C,
                   P = p)
 
    # 外推预测第k+adv项
    if (adv > 0) {
        pre.adv <- predict(k+adv-1)-predict(k+adv-2)
 
        result$predict <- pre.adv
     }
    class(result) <- 'GM1.1'
    return(result)
}
 
# 增加对应gm1.1类的泛型函数 print & plot
print.GM1.1 <- function(mod){
    cat('the result of GM(1,1)\n')
    cat('Actual Data:', '\n',mod$actual ,'\n')
    cat('Fit Data:', '\n',round(mod$fit,2) ,'\n')
    cat('Degree:', round(mod$degree,3) ,'\n')
    cat('C:', round(mod$C,3) ,'\n')
    cat('P:', round(mod$P,3) ,'\n')
    if (!is.null(mod$predict)) {
        cat('Predict Data:', round(mod$predict,2), '\n')
    }
}
 
plot.GM1.1 <- function(mod,adv=5){
    prex <- numeric(adv)
    x <- mod$actual
    for (k in 1:adv){
        prex[k] <- GM(x,k)$predict    
    }
 
    value = c(x,prex)
 
    res <- data.frame(index = 1:length(value),
                      value = value,
                      type = factor(c(rep(1,length(x)),rep(2,length(prex)))))
    library(ggplot2)
    p <- ggplot(res,aes(x=index,y= value))
    p + geom_point(aes(color=type),size=3)+ 
        geom_path(linetype=2) +
        theme_bw()
}
 
 
# 原始数据
x = c(26.7,31.5,32.8,34.1,35.8,37.5)
 
# 预测第7项
res <- GM(x,1)
print(res)
plot(res,3)

数据挖掘的三行俳句

最近才看到Tom Khabaza写的一篇很有份量的文章，阐述了数据挖掘的九大法则，在最后他以俳句方式进行了总结，可谓是字字珠玑。原文很长，只将俳句和各法则的纲要翻译放在这里。

First the business goal // Is all in data mining // This defines the field
Second is knowledge // Of business at every stage // This is the centre
Third prepare data // The form permits the question // Shape the problem space
Fourth is no free lunch / /By searching find the model // And the goal may change
Fifth is Watkins law / /There will always be patterns // Traces in data
Sixth is perception // Data mining amplifies // Thus we find insight
Seventh prediction // Generalised from data // New facts locally
Eighth is the value // All value comes from business // Not technical things
Ninth and last is change // No pattern lasts forever // We do not stand still
Miners know this truth / /How can business be better? / /If we see what is
Can the process change? // Technology can change it / /None has done so yet

1、商业目标是数据挖掘的起点
2、业务知识贯穿着数据挖掘的每个阶段
3、数据挖掘中的大部分工作是数据准备
4、对于给定的问题，只有不断尝试才能得到适合的模型
5、模式总会存在
6、数据挖掘能增强对业务的理解
7、通过归纳，预测模型增加了信息量
8、数据挖掘的价值并不取决于模型的准确或稳定
9、所有的模式都会改变

R语言资源

以前人的烦恼是没有书可读，现在人的烦恼是书太多了。关于R语言的书已经出版很多了，博主大约读过其中的四十多本，但是书在精，而不在多，学在透，而不在速。把有限的时间放到无限的书海中，这不是阅读的真意。本着造福学习者的角度，博主精选出十二本R书。什么是好书的标准？我以为是：有案例，有代码，有习题，有讲解，逻辑清楚，排版精良，体系完备，互有补充，内容千锤百炼，值得反复揣摩。书单均为英文版，都可以从网上找到。当然这份书单的选择是有主观偏见的。

XML和XPath使用方法备忘

如果把XML看作传统的关系数据库，那么XPath就是SQL。R语言中的XML包可用来解析处理XML或是HTML数据。在之前的文章中，我们了解到readHTMLTable函数，如果页面中的数据是一个规整的表格，用它是很方便的，但如果页面中是一些非结构化的数据，就要用到XML包中的其它函数了。其中最重要两个函数是xmlTreeParse()和getNodeSet()，前者负责抓取页面数据并形成树状结构，后者对抓取的数据根据XPath语法来选取特定的节点集合。下面用一个简单的例子来看一下这两个函数配合使用的效果。

如何学习数据科学

本文翻译自一篇博客文章，作者是一名软件工程师，他描述了在五年时间内学习数据科学的经历和心得，他的学习途径包括了自学（书籍、博客、小项目），课程学习，教学讨论，会议交流和工作实践。

一、入门

1）自学（2 - 4个月）

自学是起步的关键。两年前，我和几个同事组成了一个研究小组，讨论统计202课程的学习材料。这让我感觉很兴奋，并由此开始数据分析的学习研究。研究小组有5名成员，但最后只有2个人选择去更深入地研究这个领域（数据科学并不适合每一个人）。

• 学习基本的统计知识：统计202课程是非常合适的入门资料
• 学习一种统计工具：作为一个菜鸟，我用了3个月的时间埋头学习R语言，R学起来非常有趣。（为什么要学习R？）
• 解决一些好玩的小问题：好奇心是数据科学的关键。如果你对国家的经济问题，犯罪统计，体育成绩等感兴趣的话，去收集数据并开始回答你的问题吧。
• 学习Unix工具：我选择了O'Reilly出版的数据之魅作为学习材料。
• 学习SQL和脚本语言：我了解的有Java，Ruby和SQL。 Python也在我的名单上。

那些奇葩的R函数

看别人的代码会遇到一些奇葩的函数，一般的教程上很少提到，但却有很好的用处，这类函数基本上分布在base以及utils包中，下面将它们略为归纳一下，以备后用。

1，文件执行：
在用R生成一个PDF文档后，如果想去打开它，你可能会在文件夹里找到再点开。再或者我们想调用系统中的其它程序来做点事情，可能要打开cmd敲点命令。实际上这都可以在R内部完成。举例来说用pandoc转换na.md成docx再打开它。
system('pandoc d:\\rspace\\na.md -o d:\\rspace\\na.docx')
shell.exec('d:\\rspace\\na.docx')

2，网络浏览：
browseURL：浏览某个指定的网页
download.file：下载网络文件到本地

3，文件操作
dir.create：新建一个文件夹
list.dirs：显示目录下的文件夹
list.files：显示目录下的文档
file.create：文档创建
file.exists：判断文档是否存在
file.remove：文档删除
file.rename：重命名
file.append：文档添加
file.copy：文档复制
file.symlink(from, to)
file.show：显示文档内容
file.info：显示文档信息
file.edit：编辑文档
zip： 压缩文件
unzip： 解压缩文件

4，运算进度条
在一个大循环运算时，如果可以看到目前的进度是比较方便的，txtProgressBar和setTxtProgressBar函数可以帮助做到这一点，下面是内置的一个小例子：

testit <- ...="..." function="function" p="p" x="sort(runif(20)),">{
    pb <- p="p" txtprogressbar="txtprogressbar">    for(i in c(0, x, 1)) {Sys.sleep(0.5); setTxtProgressBar(pb, i)}
    Sys.sleep(1)
    close(pb)
}
testit()

用XLConnect包操控Excel表格

作为一个R迷，为什么要去捣鼓XLS文件？其实这种需求场景很多的啦，比如其它部门的同事有批量的Excel文件要处理，或者家里一把手的直接命令。Excel里面已经有不少函数可以处理数据了，包括简单的矩阵运算以及透视表什么的，但归根到底它还是需要鼠标点来点去，伤手腕啊。为了保护右手我们要提倡用代码控制一切需要鼠标的动作。高级的Excel玩家可能会用VBA去做自动处理，更高明的玩家则跳出三界外，从外部来控制单位格数据的输入输出。

R语言中有很多包可以处理表格文档，包括最为通用的RODBC包，XLConnect包也是操控Excel文档的利器，功能很丰富。使用该包的前提是要安装好Java，还要在环境变量里搞好设置。之后就可以安装加载包了。我们可以将表格文档看做是数据输入和输出端，R则是中间的运算单元。二者主要是通过数据框格式和工作表单元格进行交换。下面来看将iris数据框写入和读取的示例（其实是翻译的官方文档）。

浅谈ROC曲线

机器学习中很常见的一个大类就是二元分类器。很多二元分类器会产生一个概率预测值，而非仅仅是0-1预测值。我们可以使用某个临界点（例如0.5），以划分哪些预测为1，哪些预测为0。得到二元预测值后，可以构建一个混淆矩阵来评价二元分类器的预测效果。所有的训练数据都会落入这个矩阵中，而对角线上的数字代表了预测正确的数目，即True Positive+True Nagetive。同时可以相应算出TPR（真正率或称为灵敏度）和TNR（真负率或称为特异度）。我们主观上希望这两个指标越大越好，但可惜二者是一个此消彼涨的关系。除了分类器的训练参数，临界点的选择，也会大大的影响TPR和TNR。有时可以根据具体问题和需要，来选择具体的临界点。

2012年的学习、工作和生活

在《英雄志》里面，沉毅木讷的伍定远一直到35岁才跳出公门、踏入江湖，由此获得一系列的奇遇、成长和体验，这番际遇让人感慨。既然2012的玛雅末日没有来，就将过去一年发生的事情在第200篇博文中简单梳理一下。

学习：

最初的博客只是一个粗糙的读书笔记，但仍得到了许多同好的鼓励。赞扬的力量是强大的，于是越发的用心，花了很多时间去读去写去钻研。学习的体会可以总结为四句笨办法：精选资料，反复阅读，归纳笔记，动手操练。2012年大体上钻研了四个方面：数据挖掘、R语言、Python、数据库和SQL。

数据挖掘的算法学了一些皮毛，也知道怎么使用R的一些包，但是还需练习更底层的实现。参加了Coursera的很多公开课，但只有Machine Learning算是拿到最后的成绩，其它的都没有时间跟完进度，所谓贪多嚼不烂。R方面的书已经看得不少了，代码也敲了一些，但其实还是有一些不理解的地方，还需要深入系统的总结，2013年应该会有更进一步的成果。Python断断续续玩了半年，目前还只能用它折腾点小东西，整体感觉很好玩，应该会坚持下去。SQL的话目前工作需要用到，处理速度的确很快，不过它的思维方式和R还真是很有一些不一样，不是太感冒。

五月份和十一月份分别参加了两次R会议，见到了神交已久的统计之都各位高手，也从他们那里学到了很多东西。六月份参加了Supstat办的统计夏令营，这是我第一次参加商业培训，有点小紧张，整体还算圆满顺利。这些会议中我发现很多人都是用LaTeX来做PPT，见贤思齐嘛，回头就买了本LaTeX学习大全，日夜操练。目前正式的成果就是用LaTeX帮老婆作了一份简历，应该可以震慑住一些HR吧。

工作：

写博客的收获不止是知识，还能认识很多朋友，并在职业道路上得到他们的支持和帮助。芒果的S君是最早给我博客好评的，很是鼓舞啊。三月份的时候到上海见到了他本人和同乡L君，聊了很多，也开始了解这一行业。

八月份时候便下定决心要离开所谓的“舒适生活圈”，进入数据江湖冒险。当时还在微博上发了一个贴，原文比较煽情：“早上醒来，如系统自检一样，我知道了我是谁，我从哪来，要到哪去，心情便沉了下去。我希望今后的人生历程变得不一样，所以我想换一份工作。能和数据有关，能用到R语言更好。”

通过微博和朋友们的帮助获得了一些面试机会，因为地域、经验、专业、能力等等原因，蛮多想去的地方都没能如愿。这个过程有困扰、有磨练，也有思考和提高。最后，机会大门还是给我留了一条缝，这样我最终得以从事数据挖掘和R有关的工作，也算是从一档起步了。

生活：

2012年生活中最大的变化就是结婚了，当年达尔文曾用17页纸写下了有关结婚的系统思考。因为结婚会带来收益：例如拥有终身伴侣，避免孤独。结婚也会有成本：放弃了一定程度上的自由。两个人的生活当然会比一个人来得复杂和琐碎，这样看来是否结婚也是可以定量计算的。在我计算很长时间之后，终于让我找到了最后的稳定解，没有留下白卷。

还是用2012读书收获的一句话来结束吧：“要让宿命的人生变得有趣，就尽可能地多走几条路；要让荒谬的人生变得有意义，就不妨细细地体味每一过程的悲欢痛快，结局如何也就无关紧要了。”对我来说，新一年的生活就好象鲁滨逊踏上荒岛，虽然手中只有一把板斧，心中也无所畏惧，因为我将由此获得一系列的奇遇、成长和体验。

2013我们一起进步！