线性线性混合效应模型及R语言实现

1.A very basic tutorial for performing linear mixed effects analyses（入门极品）

Tutorial 1: http://www.bodowinter.com/tutorial/bw_LME_tutorial1.pdf

Tutorial 2: http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf

nlme (Non-Linear Mixed Effects), lme4 (Linear Mixed Effects) and asreml (average spatial real)

2.其它资料（Google搜索即有）

Mixed-Eects Models in R （较好）

An Appendix to An R Companion to Applied Regression, Second Edition

John Fox & Sanford Weisberg

Linear Mixed-Effects Models Using R（一本教材，进阶选用）

A Step-by-Step Approach

Andrzej Ga?ecki ? Tomasz Burzykowski

3.R中的线性混合模型介绍(简单了解不同的包）

http://blog.sciencenet.cn/blog-2577109-949820.html

https://www.r-bloggers.com/linear-mixed-models-in-r/

3.语法备忘

三种模型：
AOD固定斜率，DAY随机截距：LMM.model = lmer(PM25 ~ AOD + (1|Day) , data=LMMexcdata)
AOD随机斜率，DAY固定截距：LMM.model3 = lmer(PM25 ~ AOD + (0 + AOD|Day) , data=LMMexcdata)
AOD随机斜率，DAY随机截距：LMM.model2 = lmer(PM25 ~ AOD + (1 + AOD|Day) , data=LMMexcdata)

师姐发来的查看斜率和截距的程序：

sslopef=as.numeric(as.matrix(lme4::fixef(fm1)[2]))

sloper=as.numeric(LMM.model(lme4::ranef(fm1)$Day[,2]))

intersf= as.numeric(LMM.model(lme4::fixef(fm1)[1]))

intersr=as.numeric(LMM.model(lme4::ranef(fm1)$Day[,1]))

lme4::fixef(yourmodle)读取固定截距和斜率

lme4::fixef(yourmodle)读取随机截距和斜率

文档中示例（A very basic tutorial for performing linear mixed effects analyses）：

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

#load data into R

politeness=read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")

#================view the dataset======================

#show the head of politeness

head(politeness)

#show the tail of politeness

tail(politeness)

# other commands you may use

summary(politeness)

str(politeness)

colnames(politeness)

#================check for missing values======================

which(!complete.cases(politeness))

#======look at the relationship between politeness and pitch by means of a boxplot==========

boxplot(frequency ~ attitude*gender,

col=c("white","lightgray"),politeness)

#================A random intercept model======================

library(Matrix)

library(lme4)

lmer(frequency ~ attitude, data=politeness)

politeness.model = lmer(frequency ~ attitude +

gender + (1|subject) +

(1|scenario), data=politeness)

summary(politeness.model)

#Statistical significance test

politeness.null = lmer(frequency ~ gender +

(1|subject) + (1|scenario), data=politeness,

REML=FALSE)

politeness.model = lmer(frequency ~ attitude +

gender + (1|subject) + (1|scenario),

data=politeness, REML=FALSE)

anova(politeness.null,politeness.model)

#look at the coefficients of the model by subject and by item

coef(politeness.model)

#================A random slope model======================

politeness.model = lmer(frequency ~ attitude +

gender + (1+attitude|subject) +

(1+attitude|scenario),

data=politeness,

REML=FALSE)

coef(politeness.model)

#Statistical significance test (try to obtain a p-value)

#construct the null model

politeness.null = lmer(frequency ~ gender +

(1+attitude|subject) + (1+attitude|scenario),

data=politeness, REML=FALSE)

#do the likelihood ratio test

anova(politeness.null,politeness.model)

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

结果分析

1.使用数据：subject（F1,F2,F3,M3,M4,M7)，gender（M，F），scenario 语境（1～7），attitude（inf，pol），frequency（e.g.200Hz）

2.随机截距模型：

2.1 代码

politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=politeness, REML=FALSE)

其中：

Fix effect: attitude, gender（固定斜率）

Random intercepts: subject, scenario （随机截距）

＊we’ve told the model with “(1|subject)” and “(1|scenario)” to take by-subject and by-item variability into account.

2.2 输出系数分析：（注：这个系数是总系数，固定系数和随机系数已加和）

－－－－－－－－－－－－－－－－－－－－

PS. 若想分别查看固定系数和随机系数，使用以下命令

#输出固定效应系数

print(lme4::fixef(yourmodel))

#输出随机效应系数

print(lme4::ranef(yourmodel))

－－－－－－－－－－－－－－－－－－

对于随机截距模型：

Random intercepts:

scenario1～7：各自具有不同的随机截距

subject（F1,F2,F3,M3,M4,M7)：各自具有不同的随机截距

Fix effect:

对于所有项都是相同的：attitudepol ：attitudepol的pol斜率

genderM：geder的M斜率

3.随机截距模型＋随机截距模型

3.1 代码

politeness.model = lmer(frequency ~ attitude + gender + (1+attitude|subject) + (1+attitude|scenario), data=politeness, REML=FALSE)

其中：

> Fix effect:

gender为固定斜率，

> Radom effect:

让attitude随机斜率，

subject，scenario为随机截距

＊The notation “(1+attitude|subject)” means that you tell the model to expect differing baseline-levels of frequency (the intercept, represented by 1)

3.2 输出系数分析：（注：这个系数是总系数，固定系数和随机系数已加和）

> Fix effect:

genderM：geder的M斜率,固定斜率

> Random intercepts:

scenario1～7：各自具有不同的随机截距

subject（F1,F2,F3,M3,M4,M7)：各自具有不同的随机截距

attitudepol ：attitudepol的pol斜率,随机斜率

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

自己写的程序：

（读取nc文件并合并，做混合效应模型＋系数验证）

#========Set the working directory (on a Mac) and load the data========

##setwd("/Users/Gewanru/Documents/mix effect/")

#=======Read Data from a folder(循环读入多个txt文件）==========

##a = list.files("input")

##dir = paste("./input/",a,sep="")

##n = length(dir)

##LMMexcdata = read.table(file = dir[1],header=TRUE,dec = ".")

##for (i in 2:n){

## new.data = read.table(file = dir[i], header=TRUE, dec = ".")

## LMMexcdata = rbind(LMMexcdata,new.data)

##}

#read a single table

##LMMexcdata <- read.table(file = "merge.txt", header = TRUE, dec = ".")

#========view the data set（查看数据结构）=============

#see the variable names of the data

##names(LMMexcdata)

#show the head of the data

##head(LMMexcdata)

#show the data structure

##str(LMMexcdata)

#check for missing values（misiing value:"NA")

##which(!complete.cases(LMMexcdata))

#show the tail of the data

##tail(LMMexcdata)

# other commands may use

##summary(LMMexcdata)

##colnames(LMMexcdata)

#==========Read nc files（读取气象nc文件，若不用则用上一个读取方式）===============

setwd("/Users/Gewanru/Documents/mix effect/2015/Resultsnc")

library(ncdf4)

ncname <- "0101"

ncfname <- paste(ncname, ".nc", sep = "")

dname <- "PM25"

#===open a NetCDF file==

ncin <- nc_open(ncfname)

print(ncin) #Obtain some basic information of the nc file

#===Read the data from a variable in the file==

PM25data <- ncvar_get(ncin)

nc_close(ncin)

dim(PM25data)#show the the size of the array

print(PM25)

#=======Drawing some plots=============

#Scatter plot(散点图)

##x <- LMMexcdata$AOD #define x lab

##y <- LMMexcdata$PM25 #define y lab

##plot(x, y, main = "Title", xlab = "PM2.5", ylab = "AOD") # "" : mian title & Title of x(y)lab

#Box plot(箱线图)

#====================Correlation（求相关）=========================

#(Here we need to use all of the data)

#correlation

x <- LMMexcdata[,c("PM25","AOD")]

y <- LMMexcdata[,c("PM25","AOD")]

cor(x,y)

#correlation & Significance test（问题1:怎么让输出的数据位数更多一点？）

library(mnormt)

library(psych)

x <- LMMexcdata[,c("PM25","AOD")]

y <- LMMexcdata[,c("PM25","AOD")]

corr.test(x,y,use = "complete")

#===================== Linear Modle （线性回归模型）=====================

#The simplest modle(without temporal and spatial variations)

LMsimplest.model = lm(PM25 ~ AOD , data=LMMexcdata)

summary(LMsimplest.model)

coef(LMsimplest.model)

#Get the coefficient of the Linear Modle

#c1 <- LMsimplest.model$coefficient

#===========Linear Mixed effects model （混合效应模型）===================

#1.(Without site effect)

library(Matrix)

library(lme4)

LMM.model = lmer(PM25 ~ AOD + (1 + AOD|Day) , data=LMMexcdata)

summary(LMM.model)

coef(LMM.model)

print(lme4::fixef(LMM.model))

print(lme4::ranef(LMM.model))

#2.(With site effect)

##LMMsite.model = lmer(PM25 ~ AOD + (1 + AOD|Day) + (1|Site) , data=LMMexcdata)

##summary(LMMsite.model)

##coef(LMMsite.model)

#输出固定效应系数

##print(lme4::fixef(LMMsite.model))

#输出随机效应系数

##print(lme4::ranef(LMMsite.model))

#============Significance test of the model（混合效应模型显著性检验）============

LMM.null1 = lmer(PM25 ~ AOD + (0 + AOD|Day) , data=LMMexcdata, REML=FALSE)

LMM.null2 = lmer(PM25 ~ AOD + (1|Day) , data=LMMexcdata, REML=FALSE)

LMM.model = lmer(PM25 ~ AOD + (1 + AOD|Day) , data=LMMexcdata,REML=FALSE)

anova(LMM.null1,LMM.model)

anova(LMM.null2,LMM.model)

#============ Output (如果需要把系数输出成文本，自己完善下代码) ================

#Linear Mixed effects modle

#Get the coefficient of LLM

# coefficient (fixed)

AODslopefix=as.numeric((lme4::fixef(LMM.model)[2]))

interceptfix=as.numeric(as.matrix(lme4::fixef(LMM.model)[1]))

# coefficient (radom)

AODsloperad=as.numeric(as.matrix(lme4::ranef(LMM.model)$Day[,2]))

interceptsrad=as.numeric(as.matrix(lme4::ranef(LMM.model)$Day[,1]))

print(AODslopefix)

print(interceptfix)

print(AODsloperad)

print(interceptsrad)

#==============================================================

线性线性混合效应模型及R语言实现相关推荐

统计语言类教程：R语言贝叶斯统计学、Copula、SEM、极值统计学、混合效应模型、R\Python\matlab机器学习、科研数据可视化、线性回归、分位数回归、GAMS、meta分析、近红外光谱等.
查看原文>>>统计语言类教程:贝叶斯统计学.Copula.SEM.极值统计学.混合效应模型.PyTorch深度学习.科研数据可视化以下给大家整理了一些常用的统计学内容和python ...
vs合并项目_线性混合效应模型 VS 方差分析
线性混合效应模型相比方差分析具有很多优点.为了让大家更直观的了解,我们来做一下两者的对比.基于模拟的方法,生成了以下数据dat_sim: 这种数据格式是长格式,也就是从实验软件导出来的那种数据格式.这 ...
R语言线性混合效应模型（固定效应随机效应）和交互可视化3案例
最近我们被客户要求撰写关于线性混合效应模型的研究报告,包括一些图形和统计输出. 视频:线性混合效应模型(LMM,Linear Mixed Models)和R语言实现案例线性混合效应模型(LMM,Li ...
R语言线性混合效应模型不同类型的比较
在R语言中,可以使用lme4包中的lmer函数进行线性混合效应模型的分析.在分析不同类型的线性混合效应模型时,我们可以在模型公式中通过更改因变量.自变量.随机效应等参数来创建不同的模型. 下面是一个简 ...
基于R语言混合效应模型（mixed model）案例研究
全文链接: http://tecdat.cn/?p=2596 在本文中,我们描述了灵活的竞争风险回归模型.回归模型被指定为转移概率,也就是竞争性风险设置中的累积发生率(点击文末"阅读原文&q ...
R语言建立和可视化混合效应模型mixed effect model
最近我们被客户要求撰写关于混合效应模型的研究报告,包括一些图形和统计输出. 我们已经学习了如何处理混合效应模型.本文的重点是如何建立和可视化混合效应模型的结果. 相关视频:线性混合效应模型(LMM, ...
R语言潜类别混合效应模型(Latent Class Mixed Model ,LCMM)分析老年痴呆年龄数据
最近我们被客户要求撰写关于混合效应模型的研究报告,包括一些图形和统计输出. 背景和定义线性混合模型假设 N 个受试者的群体是同质的,并且在群体水平上由独特的曲线 Xi(t)β 描述.相比之下,潜在类 ...
如何用潜类别混合效应模型(Latent Class Mixed Model ,LCMM)分析老年痴呆年龄数据
全文下载链接:http://tecdat.cn/?p=24647 线性混合模型假设 N 个受试者的群体是同质的,并且在群体水平上由独特的曲线 Xi(t)β 描述(点击文末"阅读原文" ...
生物群落数据分析最常用的统计方法：回归和混合效应模型、多元统计分析技术及结构方程等数量分析方法
原文>>>R语言生物群落数据统计分析应用 R 语言作的开源.自由.免费等特点使其广泛应用于生物群落数据统计分析.生物群落数据多样而复杂,涉及众多统计分析方法.本内容以生物群落数据分析 ...
r语言面板数据回归_工具方法 | “名牌包”:面板、时间序列模型常用R语言包...
计量经济学是数学.统计技术和经济分析的综合,即运用数学.统计方法和相关经济理论,通过计量模型来揭示经济数量关系和规律.R语言包,已经实现了现代计量经济学的很多统计分析功能,下面从面板数据模型和时间序列 ...

线性线性混合效应模型及R语言实现

线性线性混合效应模型及R语言实现相关推荐

最新文章

热门文章