线性线性混合效应模型及R语⾔实现
1.A very basic tutorial for performing linear mixed effects analys（⼊门极品）
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中的线性混合模型介绍(简单了解不同的包）
3.语法备忘
三种模型：
AOD固定斜率，DAY随机截距：del = lmer(PM25 ~ AOD + (1|Day) , data=LMMexcdata)
AOD随机斜率，DAY固定截距：del3 = lmer(PM25 ~ AOD + (0 + AOD|Day) , data=LMMexcdata) AOD随机斜率，DAY随机截距：del2 = lmer(PM25 ~ AOD + (1 + AOD|Day) , data=LMMexcdata)
师姐发来的查看斜率和截距的程序：
sslopef=as.numeric(as.matrix(lme4::fixef(fm1)[2]))
sloper=as.del(lme4::ranef(fm1)$Day[,2]))
intersf= as.del(lme4::fixef(fm1)[1]))
intersr=as.del(lme4::ranef(fm1)$Day[,1]))
lme4::fixef(yourmodle)读取固定截距和斜率
lme4::fixef(yourmodle)读取随机截距和斜率
⽂档中⽰例（A very basic tutorial for performing linear mixed effects analys）：纯色手机壁纸高清
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
#load data into R
politeness=read.csv("/tutorial/politeness_data.csv")
#================view the datat======================
#show the head of politeness
head(politeness)
#show the tail of politeness
tail(politeness)
# other commands you may u
summary(politeness)
str(politeness)
colnames(politeness)
#================check for missing values======================
which(!complete.cas(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)
gender + (1|subject) +
(1|scenario), data=politeness)
del)
#Statistical significance test
politeness.null = lmer(frequency ~ gender +
(1|subject) + (1|scenario), data=politeness,
REML=FALSE)
gender + (1|subject) + (1|scenario),
data=politeness, REML=FALSE)
anova(politeness.del)
#look at the coefficients of the model by subject and by item
扶组词del)
#================A random slope model======================
gender + (1+attitude|subject) +
(1+attitude|scenario),
data=politeness,
REML=FALSE)
del)
#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.del)
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
结果分析
1.使⽤数据：subject（F1,F2,F3,M3,M4,M7)，gender（M，F），scenario 语境（1～
7），attitude（inf，pol），frequency（e.g.200Hz）
2.随机截距模型：
2.1 代码
其中：
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.2019年五行
民的成语
2.2 输出系数分析：（注：这个系数是总系数，固定系数和随机系数已加和）
－－－－－－－－－－－－－－－－－－－－
PS. 若想分别查看固定系数和随机系数，使⽤以下命令
#输出固定效应系数
print(lme4::fixef(yourmodel))
#输出随机效应系数
print(lme4::ranef(yourmodel))
－－－－－－－－－－－－－－－－－－
线性线性混合效应模型及R语⾔实现 - 影 - 我的博客
对于随机截距模型：
Random intercepts:
scenario1～7：各⾃具有不同的随机截距
subject（F1,F2,F3,M3,M4,M7)：各⾃具有不同的随机截距
Fix effect:
麦冬性味归经
对于所有项都是相同的：attitudepol ：attitudepol的pol斜率
genderM：geder的M斜率
3.随机截距模型＋随机截距模型
3.1 代码
+ (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 baline-levels of frequency (the intercept, reprented 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========
##twd("/Urs/Gewanru/Documents/mix effect/")
#=======Read Data from a folder(循环读⼊多个txt⽂件）==========
##a = list.files("input")
##dir = paste("./input/",a,p="")
##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 = "", header = TRUE, dec = ".")
#========view the data t（查看数据结构）=============
#e 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.cas(LMMexcdata))
#show the tail of the data
##tail(LMMexcdata)
# other commands may u
##summary(LMMexcdata)
##colnames(LMMexcdata)
#==========Read nc files（读取⽓象nc⽂件，若不⽤则⽤上⼀个读取⽅式）=============== twd("/Urs/Gewanru/Documents/mix effect/2015/Resultsnc")
library(ncdf4)
ncname <- "0101"
ncfname <- paste(ncname, ".nc", p = "")
dname <- "PM25"
