Selecting the Best Linear Mixed Model

合集下载

广义线性混合效应模型及其应用

广义线性混合效应模型及其应用

从 结 果 看 到 , 考 虑 了 不 同 中 心 的 差 异 , A、B 两 种 药 物 的 副 作 用 的 发 生 有 差 异 ( β1 = - 0.9298, P = 0.0354, B 药 更 易 发 生 副 作 用 ) , 不 同 实 验 中 心 的 间 的 方 差 为 1.5809, P = 0.1842, 差异无统计学意义。
设随机效应 ui 的密度函数为: fu( ui, G) , 有边际似然函数:
& Li ! β, γ" = Li ! β, ui " fu ! ui, γ" dui ui
&#% $ ni

ui
fy ! yij ui,
j=1
Xij,
β"fu ! ui,
γ" dui
( 1.5)
γ是 Ui 的方差协方差矩阵, 是 G 的参数估计值。
得到似然函数:
% L! β, γ" = Li ! β, γ" ( 1.6) i
从上式可以看到, 计算似然函数比线性混合效应模型复杂
· 2104 ·
现代预防医学 2007 年第 34 卷第 11 期 Modern Preventive Medicine, 2007, Vol.34, NO.11
得 多 , 需 要 解 决 随 机 效 应 ui 的 高 维 积 分 的 问 题 , 许 多 最 大 化 似然函数的近似的推断方法被提出, 目前积分近似方法主要有 Laplace 近似 ( Liu and Pierce, 1993) , Adaptive Gaussian 积分, 一阶泰勒 ( first- order Taylor) 序列展开近似。 2 实例分析
为了研究 A、B 两种药物的的副作用情况 , 研 究 者 随 机 选 取 了 15 个 中 心 做 临 床 实 验 , 在 每 个 中 心 中 , 随 机 抽 取 一 定 数 量 的 病 人 , 其 中 nA 个 病 人 接 受 A 药 物 , nB 个 病 人 接 受 B 药 物。数据格式见表 1。

广义线性混合效应模型在临床疗效评价中的应用

广义线性混合效应模型在临床疗效评价中的应用

广义线性混合效应模型在临床疗效评价中的应用【摘要】目的:探讨临床疗效评价中分类重复测量资料的广义线性混合效应模型(GLMMs)及的GLIMMIX宏实现。

方法:利用GLIMMIX宏ERROR和LINK语句来指示疗效指标的分布及连接函数,通过REPEATED 和RANDOM语句的TYPE选项选择合适方差协方差结构矩阵来模拟不同时间疗效指标的相关性,采用基于线性的伪似然函数进行模型参数估计。

结果:广义线性混合效应模型允许临床疗效评价指标是指数家族中任意分布,可以通过连接函数将疗效指标的均数向量与模型参数建立线性关系,简化运算过程。

结论:广义线性混合效应模型建模灵活,可为临床疗效评价提供更丰富的信息。

【关键词】广义线性混合效应模型临床疗效评价分类重复测量资料 GLIMMIX宏Apllications of Generalized Linear Mixed Models in Clinical CurativeEffects EvaluationLuo Tiane, et al Abstract Objective :To discuss generalized linear mixed models(GLMMs) of categorical repeated measurement datas in clinical curative effect evaluation, implementing with GLIMMIX macro in soft. Methods: Using the ERROR and LINK sentences of GLIMMX macro to sign the distribution and link function of the index ,adopting the TYPE option of REPEATED and RANDOM sentences to select the appropriatevariance covariance matrixs for modeling the relations, making use of pseudo likelihood function based on linear to estimate the model parameters. Results: GLMMs allow the index may be one of the exponential family (Contimuum distributions including Nomal ,beta distribution ,chi squareddistribution etc;Dispersedistributions includingBinomal ,Poisson and inverse Binomal etc), the vecor of expected means of the index is linked to the model parameters by a link function and model the linear equation, simple the calculator procedure. Conclusion: GLMMs can easily fit statistical models,the results are objective and reality, can strongly provide the abundant information for clinical curative effect evaluation. Key words generalized linear mixed models; clinical curative effects evaluation; categorical repeated measurement datas; GLIMMIX macro 临床疗效评价中常常需要对同一患者在不同时点进行多次观测并记录其疗效指标,当疗效指标为属性特征或类别时,称其为分类重复测量资料,如在治疗前、疗后4周、8周、12周等连续检测乙肝患者核心抗体,其结果有阴性、阳性两个水平;连续监测病人的治疗效果,反应变量为治愈、显效、好转、无效等。

simr包的中文名字:simr包:通过模拟计算线性混合模型的实力说明书

simr包的中文名字:simr包:通过模拟计算线性混合模型的实力说明书

Package‘simr’April13,2023Type PackageTitle Power Analysis for Generalised Linear Mixed Models by SimulationDescription Calculate power for generalised linear mixed models,usingsimulation.Designed to work with modelsfit using the'lme4'package.Described in Green and MacLeod,2016<doi:10.1111/2041-210X.12504>.Version1.0.7URL https:///pitakakariki/simrBugReports https:///pitakakariki/simr/issuesLicense GPL(>=2)Depends lme4(>=1.1-16)Importsbinom,iterators,pbkrtest,plotrix,plyr,RLRsim,stringr,stats,methods,utils,graphics,grDevices,car,lmerTest (>=3.0-0)Suggests testthat,knitr,rmarkdownLazyData trueRoxygenNote7.2.3VignetteBuilder knitrEncoding UTF-8NeedsCompilation noAuthor Peter Green[aut,cre](<https:///0000-0002-0238-9852>),Catriona MacLeod[aut],Phillip Alday[ctb]Maintainer Peter Green<********************>Repository CRANDate/Publication2023-04-1320:00:02UTC12doFit R topics documented:simr-package (2)doFit (2)doSim (3)doTest (3)extend (4)getData (5)lastResult (6)makeGlmer (6)modify (7)powerCurve (8)powerSim (9)print.powerSim (10)simdata (11)simrOptions (12)tests (13)Index16 simr-package simr:Simulation-based power calculations for mixed models.Descriptionsimr is a package that makes it easy to run simulation-based power analyses with lme4.doFit Fit model to a new response.DescriptionThis is normally an internal function,but it can be overloaded to extend simr to other packages.UsagedoFit(y,fit,subset,...)Argumentsy new values for the response variable(vector or matrix depending on the model).fit a previouslyfitted model object.subset boolean vector specifying how much of the data to use.If missing,the model isfit to all the data.This argument needs to be implemented for powerCurve towork....additional options.doSim3Valueafitted model object.doSim Generate simulated response variables.DescriptionThis is normally an internal function,but it can be overloaded to extend simr to other packages.UsagedoSim(object,...)Argumentsobject an object to simulate from,usually afitted model....additional options.Valuea vector containing simulated response values(or,for models with a multivariate response such asbinomial gl(m)m’s,a matrix of simulated response values).Suitable as input for doFit. doTest Apply a hypothesis test to afitted model.DescriptionThis is normally an internal function,but it can be overloaded to extend simr to other packages.UsagedoTest(object,test,...)Argumentsobject an object to apply a statistical test to,usually afitted model.test a test function,see tests....additional options.Valuea p-value with attributes describing the test.4extend extend Extend a longitudinal model.DescriptionThis method increases the sample size for a model.Usageextend(object,along,within,n,values)Argumentsobject afitted model object to extend.along the name of an explanatory variable.This variable will have its number of levels extended.within names of grouping variables,separated by"+"or",".Each combination of groups will be extended to n rows.n number of levels:the levels of the explanatory variable will be replaced by 1,2,3,..,n for a continuous variable or a,b,c,...,n for a factor.values alternatively,you can specify a new set of levels for the explanatory variable. Detailsextend takes"slices"through the data for each unique value of the extended variable.An extended dataset is built from n slices,with slices duplicated if necessary.ValueA copy of object suitable for doSim with an extended dataset attached as an attribute namednewData.Examplesfm<-lmer(y~x+(1|g),data=simdata)nrow(example)fmx1<-extend(fm,along="x",n=20)nrow(getData(fmx1))fmx2<-extend(fm,along="x",values=c(1,2,4,8,16))nrow(getData(fmx2))getData5 getData Get an object’s data.DescriptionGet the data associated with a model object.UsagegetData(object)getData(object)<-valueArgumentsobject afitted model object(e.g.an object of class merMod or lm).value a new data.frame to replace the old one.The new data will be stored in the newData attribute.DetailsLooks for data in the following order:1.The object’s newData attribute,if it has been set by simr.2.The data argument of getCall(object),in the environment of formula(object).ValueA data.frame with the required data.Exampleslm1<-lmer(y~x+(1|g),data=simdata)X<-getData(lm1)6makeGlmer lastResult Recover an unsaved simulationDescriptionSimulations can take a non-trivial time to run.If the user forgets to assign the result to a variable this method can recover it.UsagelastResult()See Also.Last.valueExamplesfm1<-lmer(y~x+(1|g),data=simdata)powerSim(fm1,nsim=10)ps1<-lastResult()makeGlmer Create an artificial mixed model objectDescriptionMake a merMod object with the specified structure and parameters.UsagemakeGlmer(formula,family,fixef,VarCorr,data)makeLmer(formula,fixef,VarCorr,sigma,data)Argumentsformula a formula describing the model(see glmer).family type of response variable(see family).fixef vector offixed effectsVarCorr variance and covariances for random effects.If there are multiple random ef-fects,supply their parameters as a list.data data.frame of explanatory variables.sigma residual standard deviation.modify7 modify Modifying model parameters.DescriptionThese functions can be used to change the size of a model’sfixed effects,its random effect vari-ance/covariance matrices,or its residual variance.This gives you more control over simulations from the model.Usagefixef(object)<-valuecoef(object)<-valueVarCorr(object)<-valuesigma(object)<-valuescale(object)<-valueArgumentsobject afitted model object.value new parameter values.DetailsNew values for VarCorr are interpreted as variances and covariances,not standard deviations and correlations.New values for sigma and scale are interpreted on the standard deviation scale.This means that both VarCorr(object)<-VarCorr(object)and sigma(object)<-sigma(object)leave object unchanged,as you would expect.sigma<-will only change the residual standard deviation,whereas scale<-will affect both sigma and VarCorr.These functions can be used to change the value of individual parameters,such as a singlefixed effect coefficient,using standard R subsetting commands.See AlsogetData if you want to modify the model’s data.Examplesfm<-lmer(y~x+(1|g),data=simdata)fixef(fm)fixef(fm)["x"]<--0.1fixef(fm)8powerCurve powerCurve Estimate power at a range of sample sizes.DescriptionThis function runs powerSim over a range of sample sizes.UsagepowerCurve(fit,test=fixed(getDefaultXname(fit)),sim=fit,along=getDefaultXname(fit),within,breaks,seed,fitOpts=list(),testOpts=list(),simOpts=list(),...)Argumentsfit afitted model object(see doFit).test specify the test to perform.By default,thefirstfixed effect in fit will be tested.(see:tests).sim an object to simulate from.By default this is the same as fit(see doSim).along the name of an explanatory variable.This variable will have its number of levels varied.within names of grouping variables,separated by"+"or",".Each combination of groups will be extended to n rows.breaks number of levels of the variable specified by along at each point on the power curve.seed specify a random number generator seed,for reproducible results.fitOpts extra arguments for doFit.testOpts extra arguments for doTest.simOpts extra arguments for doSim....any additional arguments are passed on to mon options in-clude:nsim:the number of simulations to run(default is1000).alpha:the significance level for the statistical test(default is0.05).progress:use progress bars during calculations(default is TRUE).powerSim9See Alsoprint.powerCurve,summary.powerCurve,confint.powerCurveExamples##Not run:fm<-lmer(y~x+(1|g),data=simdata)pc1<-powerCurve(fm)pc2<-powerCurve(fm,breaks=c(4,6,8,10))print(pc2)plot(pc2)##End(Not run)powerSim Estimate power by simulation.DescriptionPerform a power analysis for a mixed model.UsagepowerSim(fit,test=fixed(getDefaultXname(fit)),sim=fit,fitOpts=list(),testOpts=list(),simOpts=list(),seed,...)Argumentsfit afitted model object(see doFit).test specify the test to perform.By default,thefirstfixed effect in fit will be tested.(see:tests).sim an object to simulate from.By default this is the same as fit(see doSim).fitOpts extra arguments for doFit.testOpts extra arguments for doTest.simOpts extra arguments for doSim.seed specify a random number generator seed,for reproducible results.10print.powerSim ...any additional arguments are passed on to mon options in-clude:nsim:the number of simulations to run(default is1000).alpha:the significance level for the statistical test(default is0.05).progress:use progress bars during calculations(default is TRUE).See Alsoprint.powerSim,summary.powerSim,confint.powerSimExamplesfm1<-lmer(y~x+(1|g),data=simdata)powerSim(fm1,nsim=10)print.powerSim Report simulation resultsDescriptionDescribe and extract power simulation resultsUsage##S3method for class powerSimprint(x,alpha=x$alpha,level=0.95,...)##S3method for class powerCurveprint(x,...)##S3method for class powerSimsummary(object,alpha=object$alpha,level=0.95,method=getSimrOption("binom"),...)##S3method for class powerCurvesummary(object,alpha=object$alpha,level=0.95,method=getSimrOption("binom"),...simdata11 )##S3method for class powerSimconfint(object,parm,level=0.95,method=getSimrOption("binom"),alpha=object$alpha,...)##S3method for class powerCurveconfint(object,parm,level=0.95,method=getSimrOption("binom"),...)Argumentsx a powerSim or powerCurve objectalpha the significance level for the statistical test(default is that used in the call topowerSim).level confidence level for power estimate...additional arguments to pass to binom::binom.confint()alpha refers to the threshold for an effect being significant and thus directlydetermines the point estimate for the power calculation.level is the confidencelevel that is calculated for this point evidence and determines the width/coverageof the confidence interval for power.object a powerSim or powerCurve objectmethod method to use for computing binomial confidence intervals(see binom::binom.confint()) parm currently ignored,included for S3compatibility with stats::confintSee Alsobinom::binom.confint,powerSim,powerCurvesimdata Example dataset.DescriptionA simple artificial data set used in the tutorial.There are two response variables,a Poisson count zand a Gaussian response y.There is a continuous predictor x with ten values{1,2,...,10}and acategorical predictor g with three levels{a,b,c}.12simrOptionssimrOptions Options Settings for simrDescriptionControl the default behaviour of simr analyses.UsagesimrOptions(...)getSimrOption(opt)Arguments...a list of names to get options,or a named list of new values to set options.opt option name(character string).ValuegetSimrOption returns the current value for the option x.simrOptions returns1.a named list of all options,if no arguments are given.2.a named list of specified options,if a list of option names is given.3.(invisibly)a named list of changed options with their previous values,if options are set.Options in simrOptions that can be set with this method(and their default values).nsim default number of simulations(1000).alpha default confidence level(0.05).progress use progress bars during calculations(TRUE).binom method for calculating confidence intervals("exact").pbnsim number of simulations for parametric bootstrap tests using pbkrtest(100).pcmin minimum number of levels for the smallest point on a powerCurve(3).pcmax maximum number of points on the default powerCurve(10).observedPowerWarning warn if an unmodifiedfitted model is used(TRUE).carTestType type of test,i.e.type of sum of squares,for tests performed with car::Anova("II").lmerTestDdf approximation to use for denominator degrees of freedom for tests performed withlmerTest("Satterthwaite").Note that setting this option to"lme4"will reduce the lmerTestmodel to an lme4model and break functionality based on lmerTest.lmerTestType type of test,i.e.type of sum of squares,for F-tests performed with lmerTest::anova.lmerModLmerTest(2).Note that unlike the tests performed with car::Anova,the test type must be given as anumber and not a character.ExamplesgetSimrOption("nsim")oldopts<-simrOptions(nsim=5)getSimrOption("nsim")simrOptions(oldopts)getSimrOption("nsim")tests Specify a statistical test to applyDescriptionSpecify a statistical test to applyUsagefixed(xname,method=c("z","t","f","chisq","anova","lr","sa","kr","pb"))compare(model,method=c("lr","pb"))fcompare(model,method=c("lr","kr","pb"))rcompare(model,method=c("lr","pb"))random()Argumentsxname an explanatory variable to test(character).method the type of test to apply(see Details).model a null model for comparison(formula).Detailsfixed:Test a singlefixed effect,specified by xname.compare:Compare the current model to a smaller one specified by the formula model.fcompare,rcompare:Similar to compare,but only thefixed/random part of the formula needs to be supplied.random:Test the significance of a single random effect.ValueA function which takes afitted model as an argument and returns a single p-value.MethodsThe method argument can be used to specify one of the following tests.Note that"z"is an asymp-totic approximation for models notfitted with glmer and"kr"will only work with modelsfitted with lmer.z:Z-test for modelsfitted with glmer(or glm),using the p-value from summary.For modelsfit-ted with lmer,this test can be used to treat the t-values from summary as z-values,which is equivalent to assuming infinite degrees of freedom.This asymptotic approximation seems to perform well for even medium-sized data sets,as the denominator degrees of freedom are al-ready quite large(cf.Baayen et al.2008)even if calculating their exact value is analytically unsolved and computationally difficult(e.g.with Satterthwaite or Kenward-Roger approxi-mations).Setting alpha=0.045is roughly equal to the t=2threshold suggested by Baayen et al.(2008)and helps compensate for the slightly anti-conservative approximation.t:T-test for modelsfitted with lm.Also available for mixed models when lmerTest is installed,us-ing the p-value calculated using the Satterthwaite approximation for the denominator degrees of freedom by default.This can be changed by setting lmerTestDdf,see simrOptions.lr:Likelihood ratio test,using anova.f:Wald F-test,using car::eful for examining categorical terms.For modelsfitted with lmer,this should yield equivalent results to method= kr .Uses Type-II tests by default,this can be changed by setting carTestType,see simrOptions.chisq:Wald Chi-Square test,using car::Anova.Please note that while this is much faster than the F-test computed with Kenward-Roger,it is also known to be anti-conservative,especially for small es Type-II tests by default,this can be changed by setting carTestType, see simrOptions.anova:ANOV A-style F-test,using anova and lmerTest::anova.lmerModLmerTest.For‘lm‘, this yields a Type-I(sequential)test(see anova);to use other test types,use the F-tests pro-vided by car::Anova()(see above).For lmer,this generates Type-II tests with Satterthwaite denominator degrees of freedom by default,this can be changed by setting lmerTestDdf and lmerTestType,see simrOptions.kr:Kenward-Roger test,using KRmodcomp.This only applies to modelsfitted with lmer,and compares models with differentfixed effect specifications but equivalent random effects.pb:Parametric bootstrap test,using PBmodcomp.This test will be very accurate,but is also very computationally expensive.Tests using random for a single random effect call exactRLRT.ReferencesBaayen,R.H.,Davidson,D.J.,and Bates,D.M.(2008).Mixed-effects modeling with crossed random effects for subjects and items.Journal of Memory and Language,59,390–412.Exampleslm1<-lmer(y~x+(x|g),data=simdata)lm0<-lmer(y~x+(1|g),data=simdata)anova(lm1,lm0)compare(.~x+(1|g))(lm1)rcompare(~(1|g))(lm1)##Not run:powerSim(fm1,compare(.~x+(1|g)))Index.Last.value,6anova,14binom::binom.confint,11binom::binom.confint(),11car::Anova,12,14coef<-(modify),7compare(tests),13confint.powerCurve,9confint.powerCurve(print.powerSim),10 confint.powerSim,10confint.powerSim(print.powerSim),10 doFit,2,3,8,9doSim,3,4,8,9doTest,3,8,9exactRLRT,14extend,4family,6fcompare(tests),13fixed(tests),13fixef<-(modify),7getData,5,7getData<-(getData),5 getSimrOption(simrOptions),12glm,14glmer,6,14KRmodcomp,14lastResult,6lm,14lmer,14lmerTest,12,14lmerTest::anova.lmerModLmerTest,12,14 makeGlmer,6makeLmer(makeGlmer),6merMod,6modify,7PBmodcomp,14powerCurve,2,8,11,12powerSim,8,9,11print.powerCurve,9print.powerCurve(print.powerSim),10print.powerSim,10,10random(tests),13rcompare(tests),13scale<-(modify),7sigma<-(modify),7simdata,11simr-package,2simrOptions,8,10,12,14stats::confint,11summary,14summary.powerCurve,9summary.powerCurve(print.powerSim),10summary.powerSim,10summary.powerSim(print.powerSim),10tests,3,8,9,13VarCorr<-(modify),716。

Linear Mixed Effects Modeling in SPSS

Linear Mixed Effects Modeling in SPSS

Technical reportLinear Mixed-Effects Modelingin SPSS: An Introduction to the MIXED ProcedureTable of contentsIntroduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Data preparation for MIXED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1Fitting fixed-effects models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Fitting simple mixed-effects models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7Fitting mixed-effects models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .13Multilevel analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .16Custom hypothesis tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .18Covariance structure selection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .19Random coefficient models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .20Estimated marginal means. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .25References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .28About SPSS Inc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .28IntroductionThe linear mixed-effects models (MIXED) procedure in SPSS enables you to fit linear mixed-effects models to data sampled from normal distributions. Recent texts, such as those by McCulloch and Searle (2000) and Verbeke and Molenberghs (2000), comprehensively review mixed-effects models. The MIXED procedure fits models more general than those of the general linear model (GLM) procedure and it encompasses all models in the variance components (VARCOMP) procedure. This report illustrates the types of models that MIXED handles. We begin with an explanation of simple models that can be fitted using GLM and VARCOMP, to show how they are translated into MIXED. We then proceed to fit models that are unique to MIXED.The major capabilities that differentiate MIXED from GLM are that MIXED handles correlated data and unequal variances. Correlated data are very common in such situations as repeated measurements of survey respondents or experimental subjects. MIXED extends repeated measures models in GLM to allow an unequal number of repetitions. It also handles more complex situations in which experimental units are nested in a hierarchy. MIXED can, for example, process data obtained from a sample of students selected from a sample of schools in a district.In a linear mixed-effects model, responses from a subject are thought to be the sum (linear) of so-called fixed and random effects. If an effect, such as a medical treatment, affects the population mean, it is fixed. If an effect is associated with a sampling procedure (e.g., subject effect), it is random. In a mixed-effects model, random effects contribute only to the covariance structure of the data. The presence of random effects, however, often introduces correlations between cases as well. Though the fixed effect is the primary interest in most studies or experiments, it is necessary to adjust for the covariance structure of the data. The adjustment made in procedures like GLM-Univariate is often not appropriate because it assumes independence of the data.The MIXED procedure solves these problems by providing the tools necessary to estimate fixed and random effects in one model. MIXED is based, furthermore, on maximum likelihood (ML) and restricted maximum likelihood (REML) methods, versus the analysis of variance (ANOVA) methods in GLM. ANOVA methods produce an optimum estimator (minimum variance) for balanced designs, whereas ML and REML yield asymptotically efficient estimators for balanced and unbalanced designs. ML and REML thus present a clear advantage over ANOVA methods in modeling real data, since data are often unbalanced. The asymptotic normality of ML and REML estimators, furthermore, conveniently allows us to make inferences on the covariance parameters of the model, which is difficult to do in GLM.Data preparation for MIXEDMany datasets store repeated observations on a sample of subjects in “one subject per row” format. MIXED, however, expects that observations from a subject are encoded in separate rows. To illustrate, we select a subset of cases from the data that appear in Potthoff and Roy (1964). The data shown in Figure 1 encode, in one row, three repeated measurements of a dependent variable (“dist1” to “dist3”) from a subject observed at different ages (“age1” to“age3”).Figure 1. MIXED, however, requiresthat measurements at different ages be collapsed into one variable, so that each subject has three cases. The Data Restructure Wizard in SPSS simplifies the tedious data conversion process. We choose “Data->Restructure” from the pull-down menu, and select the option “Restructure selected variables into cases.” We then click the “Next” button to reach the dialog shownin Figure 2.Figure 2. We need to convert two groups ofvariables (“age” and “dist”) into cases. We therefore enter “2” and click “Next.” Thisbrings us to the “Select Variables” dialog box. Figure 3. In the “Select Variables” dialog box,we first specify “Subject ID [subid]” as the case group identification. We then enter the namesof new variables in the target variable drop-down list. For the target variable “age,” we drag “age1,” “age2,” and “age3” to the list box in the “Variables to be Transposed” group. We similarly associate variables “dist1,” “dist2,” and “dist3” with the target variable “distance.” We then drag variables that do not vary within a subject to the “Fixed Variable(s)” box. Clicking “Next” brings us to the “Create Index Variables” dialog box. We accept the default of one index variable, then click “Next” to arrive at the final dialog box.Figure 4. In the “Create One Index Variable”dialog box, we enter “visit” as the name of the indexing variable and click “Finish.” Figure 5. We now have three cases foreach subject.We can also perform the conversion using the following command syntax:VARSTOCASES/MAKE age FROM age1 age2 age3/MAKE distance FROM dist1 dist2 dist3/INDEX = visit(3)/KEEP = subid gender.The command syntax is easy to interpret—it collapses the three age variables into “age” and the three response variables into “distance.” At the same time, a new variable, “visit,” is created to index the three new cases within each subject. The last subcommand means that the two variables that are constant within a subject should be kept.Fitting fixed-effects modelsWith iid residual errorsA fitted model has the form , where is a vector of responses, is the fixed-effects design matrix, is a vector of fixed-effects parameters and is a vector of residual errors. In this model, we assume that is distributed as , where is an unknown covariance matrix. A common belief is that . We can use GLM or MIXED to fit a model with this assumption. Using a subset of the growth study dataset, we illustrate how to use MIXED to fit a fixed-effects model. The following command (Example 1) fits a fixed-effects model that investigates the effect of the variables “gender” and “age” on “distance,” which is a measure of the growth rate.Example 1: Fixed-effects model using MIXEDCommand syntax:MIXED DISTANCE BY GENDER WITH AGE/FIXED = GENDER AGE | SSTYPE(3)/PRINT = SOLUTION TESTCOV.Output:Figure 6Figure 7The command in Example 1 produces a “Type III Tests of Fixed Effects” table (Figure 6). Both “gender” and “age” are significant at the .05 level. This means that “gender” and “age” are potentially impor-tant predictors of the dependent variable. More detailed information on fixed-effects parameters may be obtained by using the subcommand /PRINT SOLUTION. The “Estimates of Fixed Effects” table (Figure 7) gives estimates of individual parameters,as well as their standard errors and confidence intervals.We can see that the mean distance for males is larger than that for females. Distance, moreover, increases with age. MIXED also produces an estimate of the residual error variance and its standard error. The /PRINT TESTCOV option gives us the Wald statistic and the confidence interval for the residual error variance estimate.Example 1 is simple—users familiar with the GLM procedure can fit the same model using GLM.Example 2: Fixed-effects model using GLM Command syntax:GLM DISTANCE BY GENDER WITH AGE /METHOD = SSTYPE(3)/PRINT = PARAMETER /DESIGN = GENDER AGE.Output:Figure 8Figure 9We see in Figure 9 that GLM and MIXED produced the same Type III tests and parameter estimates. Note, however, that in the MIXED “Type III Tests of Fixed Effects” table (Figure 6), there is no column for the sum of squares. This is because, for some complex models, the test statistics in MIXED may not be expressed as a ratio of two sums of squares. They are thus omitted from the ANOVA table.With non-iid residual errorsThe assumption may be violated in some situations. This often happens when repeated measurements are made on each subject. In the growth study dataset, for example, the response variable of each subject is measured at various ages. We may suspect that error terms within a subject are correlated. A reasonable choice of the residual error covariance will therefore be a block diagonal matrix, where each block is a first-order autoregressive (AR1) covariance matrix.Example 3: Fixed-effects model with correlated residual errors Command syntax:MIXED DISTANCE BY GENDER WITH AGE /FIXED GENDER AGE/REPEATED VISIT | SUBJECT(SUBID) COVTYPE(AR1)/PRINT SOLUTION TESTCOV R.Output:Figure 10Figure 11Figure 12Example 3 uses the /REPEATED subcommand to specify a more general covariance structure for the residual errors. Since there are three observations per subject, we assume that the set of three residual errors for each subject is a sample from a three-dimensional normal distribution with a first-order autoregressive (AR1) covariance matrix. Residual errors within each subject are therefore correlated, but are independent across subjects. The MIXED procedure, by default, usesthe REML method to estimate the covariance matrix. An alternative is to request ML estimates by using the /METHOD=ML subcommand.The command syntax in Example 3 also produces the “Residual Covariance (R) Matrix” (Figure 14), which shows the estimated covariance matrix of the residual error for one subject. We see from the “Estimates of Covariance Parameters” table (Figure 13) that the correlation parameter has a relatively large value (.729) and that the p-value of the Wald test is less than .05. The autoregressive structure may fit the data better than the model in Example 1.We also see that, for the tests of fixed effects, the denominator degrees of freedom are not integers. This is because these statistics do not have exact F distributions. The values for denominator degrees of freedom are obtained by a Satterthwaite approximation. We see in the new model that gender is not significant at the .05 level. This demonstrates that ignoring the possible correlations in your data may lead to incorrect conclusions. MIXED is therefore usually a better alternative to GLM and VARCOMP when data are correlated.Fitting simple mixed-effects models Balanced designMIXED, as its name implies, handles complicated models that involve fixed and random effects. Levels of an effect are, in some situations, only a sample of all possible levels. If we want to study the efficiency of workers in different environments, for example, we don’t need to include all workers in the study—a sample of workers is usually enough. The worker effect should be considered random, due to the sampling process. A mixed-effects model has, in general, the formwhere the extra term models the random effects. is the design matrix of random effects and is a vector of random-effects parameters. We can use GLM and MIXED to fit mixed-effects models. MIXED, however, fits a much wider class of models. To understand the functionality of MIXED, we first look at several simpler models that can be created in MIXED and GLM. We also look at the similarity between MIXED and VARCOMP in these models.Figure 13Figure 14In examples 4 through 6, we use a semiconductor dataset that appeared in Pinheiro and Bates (2000) to illustrate the similarity between GLM, MIXED, and VARCOMP. The dependent variable in this dataset is “current” and the predictor is “voltage.” The data are collected from a sample of ten silicon wafers. There are eight sites on each wafer and five measurements are taken at each site. We have, therefore, a total of 400 observations and a balanced design.Example 4: Simple mixed-effects model with balanced design using MIXEDCommand syntax:MIXED CURRENT BY WAFER WITH VOLTAGE/FIXED VOLTAGE | SSTYPE(3)/RANDOM WAFER/PRINT SOLUTION TESTCOV.Output:Figure 15Figure 16Figure 17Example 5: Simple mixed-effects model with balanced design using GLM Command syntax:GLM CURRENT BY WAFER WITH VOLTAGE /RANDOM = WAFER /METHOD = SSTYPE(3)/PRINT = PARAMETER /DESIGN = WAFER VOLTAGE.Output:Example 6: Variance components model with balanced design Command syntax:VARCOMP CURRENT BY WAFER WITH VOLTAGE /RANDOM = WAFER /METHOD = REML.Output:In Example 4, “voltage” is entered as a fixed effect and “wafer” is entered as a random effect. This example tries to model the relationship between “current” and “voltage” using a straight line, but the intercept of the regression line will vary from wafer to wafer according to a normal distribution. In the Type III tests for “voltage,” we see a significantrelationship between “current” and “voltage.” If we delve deeper into the parameter estimates table, the regression coefficient of “voltage” is 9.65. This indicates a positive relationship between “current” and “voltage.” In the “Estimates of Covariance Parameters” table (Figure 17), we have estimates for the residual error variance and the variance due to the sampling of wafers.Figure 18Figure 19Figure 20We repeat the same model in Example 5 using GLM. Note that MIXED produces Type III tests for fixed effects only, but GLM includes fixed and random effects. GLM treats all effects as fixed during computation and constructs F statistics by taking the ratio of the appropriate sums of squares. Mean squares of random effects in GLM are estimates of functions of the variance parameters of random and residual effects. These functions can be recovered from “Expected Mean Squares” (Figure 19). In MIXED, the outputs are much simpler because the variance parameters are estimated directly using ML or REML. As a result, there are no random-effect sums of squares.When we have a balanced design, as in examples 4 through 6, the tests of fixed effects are the same for GLM and MIXED. We can also recover the variance parameter estimates of MIXED by using the sum of squares in GLM. In MIXED, for example, the estimate of the residual variance is 0.175, which is the same as the MS(Error) in GLM. The variance estimate of random effect “wafer” is 0.093, which can be recovered in GLM using the “Expected Mean Squares” table (Figure 19) in Example 5:Var(WAFER) = [MS(WAFER)-MS(Error)]/40 = 0.093This is equal to MIXED’s estimate. One drawback of GLM, however, is that you cannot compute the standard error of the variance estimates.VARCOMP is, in fact, a subset of MIXED. These two procedures therefore always provide the same variance estimates, as seen in examples 4 and 6. VARCOMP only fits relatively simple models. It can only handle random effects that are iid. No statistics on fixed effects are produced. If your primary objective is to make inferences about fixed effects and your data are correlated, MIXED is a better choice.An important note: Due to the different estimation methods that are used, GLM and MIXED often do not produce the same results. The next section gives an example of situations in which they produce different results.Unbalanced designOne situation in which MIXED and GLM disagree is with an unbalanced design. To illustrate this, we removed some cases in the semiconductor dataset, so that the design is no longer balanced.Figure 21We then rerun examples 4 through 6 with this unbalanced dataset. The output is shown in examples 4a through 6a. We want to see whether the three methods—GLM, MIXED and VARCOMP—still agree with each other.Example 4a: Mixed-effects model with unbalanced design using MIXED Command syntax:MIXED CURRENT BY WAFER WITH VOLTAGE /FIXED VOLTAGE | SSTYPE(3)/RANDOM WAFER/PRINT SOLUTION TESTCOV.Output:Example 5a: Mixed-effects model with unbalanced design using GLM Command syntax:GLM CURRENT BY WAFER WITH VOLTAGE /RANDOM = WAFER /METHOD = SSTYPE(3)/PRINT = PARAMETER /DESIGN = WAFER VOLTAGE.Figure 22Figure 23Figure 24Output:Example 6a: Variance components model with unbalanced designCommand syntax:VARCOMP CURRENT BY WAFER WITH VOLTAGE /RANDOM = WAFER/METHOD = REML.Output:Since the data have changed, we expect examples 4a through 6a to differ from examples 4 through 6. We will focus instead on whether examples 4a, 5a, and 6a agree with each other.In Example 4a, the F statistic for the “voltage” effect is 67481.118, butExample 5a gives an F statistic value of 67482.629. Apart from the test of fixed effects, we also see a difference in covariance parameter estimates.Examples 4a and 6a, however, show that VARCOMP and MIXED can produce the same variance estimates, even in an unbalanced design. This is because MIXED and VARCOMP offer maximum likelihood or restricted maximum likelihood methods in estimation, while GLM estimates are based on the method-of-moments approach.MIXED is generally preferred because it is asymptotically efficient (minimum variance), whether or not the data are balanced. GLM, however, only achieves its optimum behavior when the data are balanced.Figure 25Figure 26Figure 27Fitting mixed-effects modelsWith subjectsIn the semiconductor dataset, “current” is a dependent variable measured on a batch of wafers. These wafers are therefore considered subjects in a study. An effect of interest (such as “site”) may often vary with subjects (“wafer”). One scenario is that the (population) means of “current” at separate sites are different. When we look at the current measured at these sites on individual wafers, however, they hover below or above the population mean according to some normal distribution. It is therefore common to enter an “effect by subject” interaction term in a GLM or MIXED model to account for the subject variations.In the dataset there are eight sites and ten wafers. The site*wafer effect, therefore, has 80 parameters, which can be denoted by , i=1...10 and j=1...8. A common assumption is that ’s are assumed to be iid normal with zero mean and an unknown variance. The mean is zero because ’s are used to model only the population variation. The mean of the population is modeled by entering “site” as a fixed effect in GLM and MIXED. The results of this model for MIXED and GLM are shown in examples 7 and 8.Example 7: Fitting random effect*subject interaction using MIXEDCommand syntax:MIXED CURRENT BY WAFER SITE WITH VOLTAGE/FIXED SITE VOLTAGE |SSTYPE(3)/RANDOM SITE*WAFER | COVTYPE(ID).Output:Figure 28Figure 29Example 8: Fitting random effect*subject interaction using GLMCommand syntax:GLM CURRENT BY WAFER SITE WITH VOLTAGE/RANDOM = WAFER/METHOD = SSTYPE(3)/DESIGN = SITE SITE*WAFER VOLTAGE.Output:Figure 30Figure 31Since the design is balanced, the results of GLM and MIXED in examples 7 and 8 match. This is similar to examples 4 and 5. We see from the results of Type III tests that “voltage” is still an important predictor of “current,” while “site” is not. The mean currents at different sites are thus not significantly different from each other, so we can use a simpler model without the fixed effect “site.” We should still, however, consider a random-effects model, because ignoring the subject variation may lead to incorrect standard error estimates of fixed effects or false significant tests.Up to this point, we examined primarily the similarities between GLM and MIXED. MIXED, in fact, has a much more flexible way of modeling random effects. Using the SUBJECT and COVTYPE options, Example 9 presents an equivalent form of Example 7.Example 9: Fitting random effect*subject interaction using SUBJECT specificationCommand syntax:MIXED CURRENT BY SITE WITH VOLTAGE/FIXED SITE VOLTAGE |SSTYPE(3)/RANDOM SITE | SUBJECT(WAFER) COVTYPE(ID).The SUBJECT option tells MIXED that each subject will have its own set of random parameters for the random effect “site.” The COVTYPE option will specify the form of the variance covariance matrix of the random parameters within one subject. The command syntax attempts to specify the distributional assumption in a multivariate form, which can be written as:Under normality, this assumption is equivalent to that in Example 7. One advantage of the multivariate form is that you can easily specify other covariance structures by using the COVTYPE option. The flexibility in specifying covariance structures helps us to fit a model that better describes the data. If, for example, we believe that the variances of different sites are different, we can specify a diagonal matrix as covariance type and the assumption becomes:The result of fitting the same model using this assumption is given in Example 10.Example 10: Using COVTYPE in a random-effects model Command syntax:MIXED CURRENT BY SITE WITH VOLTAGE /FIXED SITE VOLTAGE |SSTYPE(3)/RANDOM SITE | SUBJECT(WAFER) COVTYPE(DIAG)/PRINT G TESTCOV.Output:Figure 32Figure 33Figure 34In Example 10, we request one extra table, the estimated covariance matrix of the random effect “site.” It is an eight-by-eight diagonal matrix in this case. Note that changing the covariance structure of a random effect also changes the estimates and tests of fixed effects. We want, in practice, an objective method to select suitable covariance struc-tures for our random effects. In the section “Covariance Structure Selection,” we revisit examples 9 and 10 to show how to selectcovariance structures for random effects.Multilevel analysisThe use of the SUBJECT and COVTYPE options in /RANDOM and /REPEATED brings many options for modeling the covariance structures of random effects and residual errors. It is particularly useful when modeling data obtained from a hierarchy. Example 11 illustrates the simultaneous use of theseoptions in a multilevel model. We selected data from six schools from the Junior School Project of Mortimore, et al. (1988). We investi-gate below how the socioeconomic status (SES) of a student affects his or her math scores over a three-year period.Example 11: Multilevel mixed-effects model Command syntax:MIXED MATHTEST BY SCHOOL CLASS STUDENT GENDER SES SCHLYEAR /FIXED GENDER SES SCHLYEAR SCHOOL/RANDOM SES |SUBJECT(SCHOOL*CLASS) COVTYPE(ID)/RANDOM SES |SUBJECT(SCHOOL*CLASS*STUDENT) COVTYPE(ID)/REPEATED SCHLYEAR | SUBJECT(SCHOOL*CLASS*STUDENT) COVTYPE(AR1)/PRINT SOLUTION TESTCOV.Figure 35Figure 36In Example 11, the goal is to discover whether socioeconomic status (“ses”) is an important predictor for mathematics achievement (“mathtest”). To do so, we use the factor “ses” as a fixed effect. We also want to adjust for the possible sampling variation due to different classes and students. “Ses” is therefore also used twice as a random effect. The first random effect tries to adjust for the variation of the “ses” effect owing to class variation. In order to identify all classes in the dataset, school*class is specified in the SUBJECT option. The second random effect also tries to adjust for the variation of the “ses” effect owing to student variation. The subject specification is thus school*class*student. All of the students are followed for three years; the school year (“schlyear”) is therefore used as a fixed effect to adjust for possible trends in this period. The /REPEATED subcommand is also used to model the possible correlation of the residual errors within each student.We have a relatively small dataset. Since there are only six schools, we can only use school as a fixed effect while adjusting for possible differences between schools. In this example, there is only one random effect at each level. With SPSS 11.5 or later, you can specify more than one random effect in MIXED. If multiple random effects are specified on the same RANDOM subcommand, you can model their correlation by using a suitable COVTYPE specification. If the random effects are specified on separate RANDOM subcommands, they are assumed to be independent.Figure 38Figure 39In the Type III tests of fixed effects, in Example 11, we see that socioeconomic status does impact student performance. The parameter estimates of “ses” for students with “ses=1” (fathers have managerial or professional occupations) indicate that these students perform better than students at other socioeconomic levels. The effect “schlyear” is also significant in the model and the students’ performances increase with “schlyear.”From “Estimates of Covariance Parameters” (Figure 39), we notice that the estimate of the “AR1 rho” parameter is not significant, which means that a simple, scaled-identity structure may be used. For the variation of “ses” due to school* class, the estimate is very small compared to other sources of variance and the Wald test indicates that it is not significant. We can therefore consider removing the random effect from the model.We see from this example that the major advantages of MIXED are that it is able to look at different aspects of a dataset simultaneously and that all of the statistics are already adjusted for all effects in the model. Without MIXED, we must use different tools to study different aspects of the models. An example of this is using GLM to study the fixed effects andusing VARCOMP to study the covariance structure. This is not only time consuming, but the assumptions behind the statistics are usually violated.Custom hypothesis testsApart from predefined statistics, MIXED allows users to construct custom hypotheses on fixed- and random-effects parameters through the use of the /TEST subcommand. To illustrate, we use a dataset from Pinheiro and Bates (2000). The data consistof a CT scan on a sample of ten dogs. The dogs’ left and right lymph nodes were scanned and the intensity of each scan was recorded in the variable pixel. The following mixed-model command syntax tests whether there is a difference between the left and right lymph nodes.Example 12: Custom hypothesis testing in mixed-effects modelCommand syntax:MIXED PIXEL BY SIDE/FIXED SIDE/RANDOM SIDE | SUBJECT(DOG) COVTYPE(UN)/TEST(0) ‘Side (fixed)’ SIDE 1 -1/TEST(0) ‘Side (random)’ SIDE 1 -1 | SIDE 1 -1/PRINT LMATRIX.Output:Figure 40The output of the two /TEST subcommands is shown above. The first test looks at differences in the left and right sides in the general population (broad inference space). We should use the second test to test the differences between the left and right sides for the sample of dogs used in this particular study (narrow inference space). In the second test, the average differences of the random effects over the ten dogs are added to the statistics. MIXED automatically calculates the average over subjects. Note that the contrast coefficients for random effects are scaled by one/(number of subjects). Though the average difference for the random effect is zero, it affects the standard error of the statistic. We see that statistics of the two tests are the same, but the second has a smaller standard error. This means that if we make an inference on a larger population, there will be more uncertainty. This is reflected in the larger standard error of the test. The hypothesis in this example is not significant in the general population, but it is significant for the narrow inference. A larger sample size is therefore often needed to test a hypothesis about the general population.Covariance structure selectionIn examples 3 and 11, we see the use of Wald statistics in covariance structure selection. Another approach to testing hypotheses on covariance parameters uses likelihood ratio tests. The statistics are constructed by taking the differences of the -2 Log likelihoods of two nested models. Under the null hypothesis that the covariance parameters are 0 in the population, this difference follows a chi-squared distribution with degrees of freedom equal to the difference in the number of parameters of the models.To illustrate the use of the likelihood ratio test, we again look at the model in examples 9 and 10. In Example 9, we use a scaled identity as the covariance matrix of the random effect “site.” In Example 10, however, we use a diagonal matrix with unequal diagonal elements. Our goal is to discover which model better fits the data. We obtain the -2 Log likelihood values and other criteria about the two models from the information criteria tables shown on the next page.Figure 41Figure 42。

多水平统计分析模型(混合效应模型)

多水平统计分析模型(混合效应模型)

多⽔平统计分析模型(混合效应模型)⼀、概述普通的线性回归只包含两项影响因素,即固定效应(fixed-effect)和噪声(noise)。

噪声是我们模型中没有考虑的随机因素。

⽽固定效应是那些可预测因素,⽽且能完整的划分总体。

例如模型中的性别变量,我们清楚只有两种性别,⽽且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model?因为有些现实的复杂数据是普通线性回归是处理不了的。

例如我们对⼀些⼈群进⾏重复测量,此时存在两种随机因素会影响模型,⼀种是对某个⼈重复测试⽽形成的随机噪声,另⼀种是因为⼈和⼈不同⽽形成的随机效应(random effect)。

如果将⼀个⼈的测量数据看作⼀个组,随机因素就包括了组内随机因素(noise)和组间随机因素(random effect)。

这种嵌套的随机因素结构违反了普通线性回归的假设条件。

你可能会把⼈员(组间的随机效应)看作是⼀种分类变量放到普通线性回归模型中,但这样作是得不偿失的。

有可能这个factor的level很多,可能会⽤去很多⾃由度。

更重要的是,这样作没什么意义。

因为⼈员ID和性别不⼀样,我们不清楚它的意义,⽽且它也不能完整的划分总体。

也就是说样本数据中的路⼈甲,路⼈⼄不能完全代表总体的⼈员ID。

因为它是随机的,我们并不关⼼它的作⽤,只是因为它会影响到模型,所以不得不考虑它。

因此对于随机效应我们只估计其⽅差,不估计其回归系数。

混合模型中包括了固定效应和随机效应,⽽随机效应有两种⽅式来影响模型,⼀种是对截距影响,⼀种是对某个固定效应的斜率影响。

前者称为 Random intercept model,后者称为Random Intercept and Slope Model。

Random intercept model的函数结构如下Yij = a0 + a1*Xij + bi + eija0: 固定截距a1: 固定斜率b: 随机效应(只影响截距)X: 固定效应e: 噪声混合线性模型有时⼜称为多⽔平线性模型或层次结构线性模型由两个部分来决定,固定效应部分+随机效应部分,⼆、R语⾔中的线性混合模型可⽤包1、nlme包这是⼀个⽐较成熟的R包,是R语⾔安装时默认的包,它除了可以分析分层的线性混合模型,也可以处理⾮线性模型。

generalized additive mixed modeling

generalized additive mixed modeling

generalized additive mixed modeling1. 引言1.1 概述在统计建模中,回归模型是一种常见的分析工具,用于研究变量之间的关系。

然而,传统的回归模型通常对数据的线性关系做出了限制,无法很好地拟合复杂的非线性关系。

为了解决这个问题,广义可加混合模型(Generalized Additive Mixed Modeling, GAMM)应运而生。

GAMM是一种灵活而强大的统计建模方法,它结合了广义可加模型(Generalized Additive Model, GAM)和混合效应模型(Mixed Effects Model)。

通过引入非线性平滑函数和随机效应,GAMM能够更准确地描述变量之间的复杂关系,并考虑到数据中可能存在的随机变异。

本文将详细介绍GAMM的理论基础、模型框架和参数估计方法。

同时,我们还将探讨GAMM在各个领域中的应用,并与传统回归模型以及混合效应模型进行比较和评估。

最后,我们将总结目前对于GAMM方法的认识,并提出未来研究方向。

1.2 文章结构本文共分为五个部分。

首先,在引言部分概述了GAMM的背景和研究意义。

接下来,第二部分将介绍GAMM的理论基础、模型框架和参数估计方法。

第三部分将详细探讨GAMM在生态学、社会科学和医学研究中的应用案例。

第四部分将与其他回归模型和传统混合模型进行比较,并对GAMM方法的优缺点及局限性进行讨论。

最后,在第五部分中,我们将总结全文的主要内容,并提出对未来研究方向的建议。

1.3 目的本文旨在全面介绍广义可加混合模型(GAMM)这一统计建模方法,以及其在不同领域中的应用。

通过对GAMM的理论基础、模型框架和参数估计方法进行详细描述,读者可以了解到该方法如何解决传统回归模型无法处理非线性关系问题的局限性。

同时,通过实际案例研究,读者可以进一步了解GAMM在生态学、社会科学和医学研究等领域中的应用效果。

此外,通过与其他回归模型和传统混合模型进行比较,本文还旨在评估GAMM方法的优势和局限性。

Generalized Linear Mixed Models

Generalized Linear Mixed Models

glmm
Department of Biostatistics University of Copenhagen
Example 1: Smoking in School Children
library(foreign) library(lme4) dd <- read.xport("~/niss.xpt") names(dd) <- tolower(names(dd)) dd$ryger <- dd$rygvan<2 (fit <- lmer(ryger~koen + alder + factor(klastrin) + factor(skf + (1|skole)+(1|klassenr),family=binomial,data=dd)) anova(fit) waldtest <- function(fit,ix) {x <- fixef(fit)[ix]; V <- as.matrix(vcov(fit))[ix,ix]; x%*%solve(V,x)} waldtest(fit, 4:5)
glmm
Department of Biostatistics University of Copenhagen
ni
fij (yij |bi , β, φ)f (bi |D)dbi
j=1
They proceed to explain the main approximations as Approximating the integrand (Laplace) Approximating the data (PQL) Approximating the integral

Bayesian Linear Mixed-Effects Models 贝叶斯线性混合效应模型说

Bayesian Linear Mixed-Effects Models  贝叶斯线性混合效应模型说

Package‘blme’October12,2022Version1.0-5Date2020-12-28Title Bayesian Linear Mixed-Effects ModelsDepends R(>=3.0-0),lme4(>=1.0-6)Imports methods,stats,utilsSuggests expint(>=0.1-3),testthatDescription Maximum a posteriori estimation for linear and generalized linear mixed-effects mod-els in a Bayesian setting,implementing the meth-ods of Chung,et al.(2013)<doi:10.1007/s11336-013-9328-2>.Extends pack-age'lme4'(Bates,Maechler,Bolker,and Walker(2015)<doi:10.18637/jss.v067.i01>).License GPL(>=2)URL https:///vdorie/blmeBugReports https:///vdorie/blme/issuesNeedsCompilation noAuthor Vincent Dorie[aut,cre](<https:///0000-0002-9576-3064>), Douglas Bates[ctb](lme4non-modular functions,<https:///0000-0001-8316-9503>),Martin Maechler[ctb](lme4non-modular functions,<https:///0000-0002-8685-9910>),Ben Bolker[ctb](lme4non-modular functions,<https:///0000-0002-2127-0443>),Steven Walker[ctb](lme4non-modular functions,<https:///0000-0002-4394-9078>)Maintainer Vincent Dorie<****************>Repository CRANDate/Publication2021-01-0518:10:11UTCR topics documented:blme (2)bmerDist-class (5)bmerMod-class (7)12blme Index9blme Fit Bayesian Linear and Generalized Linear Mixed-Effects ModelsDescriptionMaximum a posteriori estimation for linear and generalized linear mixed-effects models in a Bayesian setting.Built off of lmer.Usageblmer(formula,data=NULL,REML=TRUE,control=lmerControl(),start=NULL,verbose=0L,subset,weights,na.action,offset,contrasts=NULL,devFunOnly=FALSE,cov.prior=wishart,fixef.prior=NULL,resid.prior=NULL,...)bglmer(formula,data=NULL,family=gaussian,control=glmerControl(),start=NULL,verbose=0L,nAGQ=1L,subset,weights,na.action,offset,contrasts=NULL,mustart,etastart,devFunOnly=FALSE,cov.prior=wishart,fixef.prior=NULL,...)Argumentscov.prior a BLME prior or list of priors with allowable distributions:wishart,invwishart, gamma,invgamma,or NULL.Imposes a prior over the covariance of the randomeffects/modeled coefficients.Default is wishart.The NULL argument imposesflat priors over all relevant parameters.fixef.prior a BLME prior of family normal,t,horseshoe,or NULL.Imposes a prior overthefixed effects/modeled coefficients.Default is NULL.resid.prior a BLME prior of family gamma,invamma,point or NULL.Imposes a prior overthe noise/residual variance,also known as common scale parameter or the con-ditional variance given the random effects.Default is NULL.start like the start arguments for lmer and glmer a numeric vector or named list.Unlike the aforementioned,list members of fixef and sigma are applicable tolinear mixed models provided that numeric optimization is required for theseparameters.formula,data,REML,family,control,verbose,nAGQ,mustart,etastart,devFunOnly,...model specification arguments as in lmer and glmer;see there for details.subset,weights,na.action,offset,contrastsfurther model specification arguments as in lm;see there for details.blme3 DetailsThe bulk of the usage for blmer and bglmer closely follows the functions lmer and glmer.Those help pages provide a good overview offitting linear and generalized linear mixed models.The primary distinction is that blmer and bglmer allow the user to do Bayesian inference or penalized maximum likelihood,with priors imposed on the different model components.For the specifics of any distribution listed below,see the distributions page.Covariance PriorThe cov.prior argument applies a prior over the covariance matrix of the random effects/modeled coefficients.As there is one covariance matrix for every named grouping factor-that is every element that appears to the right of a vertical bar("|")in the model formula-it is possible to apply as many different priors as there are said factors.The general formats of an argument to blmer or bglmer for such a prior are of the form:•cov.prior=~covariance.distribution(option1=value1,...)•cov.prior=list(fc.nm~dist1,fc.nm~dist2,...,default.distribution) If the“~”construct is ommitted,the prior is interpretted as a default and applied to all factors that lack specific priors of their own.Options are not required,but permitfine-tuning of the model.Supported distributions are gamma,invgamma,wishart,invwishart,NULL,and custom.The common.scale option,a logical,determines whether or not the prior applies to in the absolute-real world sense(value=FALSE),or if the prior is applied to the random effect covariance divided by the estimated residual variance(TRUE).As a practical matter,when false computation can be slower as the profiled common scale may no longer have a closed-form solution.As such,the default for all cases is TRUE.Other options are specified along with the specific distributions and defaults are explained in the blme distributions page.Fixed Effects PriorPriors on thefixed effects,or unmodeled coefficients,are specified in a fashion similar to that of covariance priors.The general format is•fixef.prior=multivariate.distribution(options1=value1,...) At present,the implemented multivariate distributions are normal,t,horseshoe,and NULL.t and horseshoe priors cannot be used when REML is TRUE,as that integral does not have a closed form solution.Residual Variance PriorThe general format for a residual variance prior is the same as for afixed effect prior.The supported distributions are point,gamma,invgamma.ValueAn object of class"bmerMod",for which many methods are available.See there for details.4blmeSee Alsolmer,glmer,merMod class,and lm.Examplesdata("sleepstudy",package="lme4")###Examples using a covariance prior###Here we are ignoring convergence warnings just to illustate how the package#is used:this is not a good idea in practice..control<-lmerControl(check.conv.grad="ignore")(fm1<-blmer(Reaction~Days+(0+Days|Subject),sleepstudy,control=control,cov.prior=gamma))(fm2<-blmer(Reaction~Days+(0+Days|Subject),sleepstudy,control=control,cov.prior=gamma(shape=2,rate=0.5,posterior.scale= sd ))) (fm3<-blmer(Reaction~Days+(1+Days|Subject),sleepstudy,control=control,cov.prior=wishart))(fm4<-blmer(Reaction~Days+(1+Days|Subject),sleepstudy,control=control,cov.prior=invwishart(df=5,scale=diag(0.5,2))))#Custom priorpenaltyFn<-function(sigma)dcauchy(sigma,0,10,log=TRUE)(fm5<-blmer(Reaction~Days+(0+Days|Subject),sleepstudy,cov.prior=custom(penaltyFn,chol=TRUE,scale="log")))###Examples using a fixed effect prior###(fm6<-blmer(Reaction~Days+(1+Days|Subject),sleepstudy,cov.prior=NULL,fixef.prior=normal))(fm7<-blmer(Reaction~Days+(1+Days|Subject),sleepstudy,cov.prior=NULL,fixef.prior=normal(cov=diag(0.5,2),common.scale=FALSE))) ###Example using a residual variance prior####This is the"eight schools"data set;the mode should be at the boundary#of the space.control<-lmerControl(check.conv.singular="ignore",check.nobs.vs.nRE="ignore",check.nobs.vs.nlev="ignore")y<-c(28,8,-3,7,-1,1,18,12)sigma<-c(15,10,16,11,9,11,10,18)g<-1:8(schools<-blmer(y~1+(1|g),control=control,REML=FALSE,resid.prior=point,cov.prior=NULL,weights=1/sigma^2))bmerDist-class Bayesian Linear Mixed-Effects Model Prior Representations andbmer*Dist MethodsDescriptionObjects created in the initialization step of a blme model that represent the type of prior being applied.Objects from the ClassObjects can be created by calls of the form new("bmerPrior",...)or,more commonly,as side effects of the blmer and bglmer functions.When using the main blme functions,the prior-related arguments can be passed what essentially are function calls with the distinction that they are delayed in evaluation until information about the model is available.At that time,the functions are defined in a special environment and then evaluated in an environment that directly inherits form the one in which blmer or bglmer was called.This is reflected in some of the prototypes of various prior-creating functions which depend on parameters not available in the top-level environment.Finally,if the trailing parentheses are omitted from a blmer/bglmer prior argument,they are simply added as a form of“syntactic sugar”.Prior DistributionsThis section lists the prototypes for the functions that are called to parse a prior during a modelfit.Fixed Effect Priors•normal(sd=c(10,2.5),cov,common.scale=TRUE)Applies a Gaussian prior to thefixed effects.Normal priors are constrained to have a mean of 0-non-zero priors are equivalent to shifting covariates.The covariance hyperparameter can be specified either as a vector of standard deviations, using the sd argument,a vector of variances using the cov argument,or the entire vari-ance/covariance matrix itself.When specifying standard deviations,a vector of length less than the number offixed effects will have its tail repeated,while thefirst element is assumed to apply only to the intercept term.So in the default of c(10,2.5),the intercept receives a standard deviation of10and the various slopes are all given a standard deviation of2.5.The common.scale argument specifies whether or not the prior is to be interpretted as being on the same scale as the residuals.To specify a prior in an absolute sense,set to FALSE.Argument is only applicable to linear mixed models.•t(df=3,mean=0,scale=c(10^2,2.5^2),common.scale=TRUE)The degrees of freedom-df argument-must be positive.If mean is of length1,it is repeated for everyfixed effect.Length2repeats just the second element for all slopes.Otherwise,the length must be equal to that of the number offixed effects.If scale is of length1,it is repeated along the diagonal for every component.Length2repeatsjust the second element for all slopes.Length equal to the number offixed effects sees thevector simply turned into a diagonal matrix.Finally,it can be a full scale matrix,so long as itis positive definite.t priors for linear mixed models require that thefixed effects be added to set of parameters thatare numerically optimized,and thus can substantially increase running time.In addition,whencommon.scale is TRUE,the residual variance must be numerically optimized as well.normalpriors on the common scale can be fully profiled and do not suffer from this drawback.At present,t priors cannot be used with the REML=TRUE argument as that implies an integralwithout a closed form solution.•horseshoe(mean=0,global.shrinkage=2.5,common.scale=TRUE)The horseshoe shrinkage prior is implemented similarly to the t prior,in that it requires addingthefixed effects to the parameter set for numeric optimization.global.shrinkage,also referred to asτ,must be positive and is on the scale of a standarddeviation.Local shrinkage parameters are treated as independent across allfixed effects andare integrated out.See Carvalho et al.(2009)in the references.Covariance Priors•gamma(shape=2.5,rate=0,common.scale=TRUE,posterior.scale="sd")Applicable only for univariate grouping factors.A rate of0or a shape of0imposes an im-proper prior.The posterior scale can be"sd"or"var"and determines the scale on which theprior is meant to be applied.•invgamma(shape=0.5,scale=10^2,common.scale=TRUE,posterior.scale="sd")Applicable only for univariate grouping factors.A scale of0or a shape of0imposes animproper prior.Options are as above.•wishart(df=level.dim+2.5,scale=Inf,common.scale=TRUE,posterior.scale="cov")A scale of Inf or a shape of0imposes an improper prior.The behavior for singular matriceswith only some infinite eigenvalues is undefined.Posterior scale can be"cov"or"sqrt",thelatter of which applies to the unique matrix root that is also a valid covariance matrix.•invwishart(df=level.dim-0.5,scale=diag(10^2/(df+level.dim+1),level.dim), common.scale=TRUE,posterior.scale="cov")A scale of0or a shape of0imposes an improper prior.The behavior for singular matriceswith only some zero eigenvalues is undefined.•custom(fn,chol=FALSE,common.scale=TRUE,scale="none")Applies to the given function(fn).If chol is TRUE,fn is passed a right factor of covariancematrix;FALSE results in the matrix being passed directly.scale can be"none","log",or"dev"corresponding to p(Σ),log p(Σ),and−2log p(Σ)respectively.Since the prior is may have an arbitrary form,setting common.scale to FALSE for a linearmixed model means that full profiling may no longer be possible.As such,that parameter isnumerically optimized.Residual Variance Priors•point(value=1.0,posterior.scale="sd")Fixes the parameter to a specific value given as either an"sd"or a"var".•gamma(shape=0,rate=0,posterior.scale="var")As above with different defaults.•invgamma(shape=0,scale=0,posterior.scale="var")As above with different defaults.Evaluating EnvironmentThe variables that the defining environment have populated are:•p aliased to n.fixef-the number offixed effects•n aliased to n.obs-the number of observations•q.k aliased to level.dim-for covariance priors,the dimension of the grouping factor/grouping level•j.k aliased to n.grps-also for covariance priors,the number of groups that comprise aspecific grouping factorMethodstoString Pretty-prints the distribution and its parameters.ReferencesCarvalho,Carlos M.,Nicholas G.Polson,and James G.Scott."Handling Sparsity via the Horse-shoe."AISTATS.V ol.5.2009.See Alsoblmer()and bglmer(),which produce these objects,and bmerMod-class objects which contain them.bmerMod-class Class"bmerMod"of Fitted Mixed-Effect ModelsDescriptionThe bmerMod class represents linear or generalized linear or nonlinear mixed-effects models with possible priors over model components.It inherits from the merMod class.Objects from the ClassObjects are created by calls to blmer or bglmer.SlotsA bmerMod object contains one additional slot beyond the base merMod class:priors:A named list comprised of covPriors,fixefPrior,and residPrior.In addition,the devcomp slot,element cmp includes the penalty item which is the computed de-viance for the priors.Add this to the regular deviance to obtain the value of the objective function that is used in optimization.See Alsoblmer and bglmer,which produce these objects.merMod,from which this class inherits.ExamplesshowClass("bmerMod")methods(class="bmerMod")Index∗GLMMblme,2∗NLMMblme,2∗classesbmerDist-class,5bmerMod-class,7∗methodsblme,2∗modelsblme,2bglmer,5,7,8bglmer(blme),2bglmerMod-class(bmerMod-class),7 blme,2blmer,5,7,8blmer(blme),2blmerMod-class(bmerMod-class),7 bmerDist(bmerDist-class),5bmerDist-class,5bmerMod,3bmerMod(bmerMod-class),7bmerMod-class,7distributions,3glmer,2–4lm,2,4lmer,2–4merMod,4,7,8print,bmerDist-method(bmerDist-class), 5print.bmerMod(bmerMod-class),7 print.summary.bmerMod(bmerMod-class),7 prior,2show,bmerDist-method(bmerDist-class),5show,bmerMod-method(bmerMod-class),7show.bmerMod(bmerMod-class),7summary.bmerMod(bmerMod-class),7summary.summary.bmerMod(bmerMod-class),7vcov.summary.bmerMod(bmerMod-class),7 9。

线性混合模型对作物育种无重复试验数据分析的 实证研究

线性混合模型对作物育种无重复试验数据分析的 实证研究

图 1 具有 99 个品系的试验 1 田间分布图 Fig .1 Field layout of trial 1 with 99 lines
3
CK2 S1 … S5 CK2 … … … CK2 S46 … S50 CK2
图 2 具有 50 个品系的试验 2 田间分布图 Fig .2 Field layout of trial 2 with 50 lines
[8-10] [6] [7]

近年来,随着线性混合模型理论的完善和统计分析软件的不断发展,为传统上不能进行统计分析的有关试 验提供了可能的分析途径。 根据线性混合模型分析原理,对以往不能进行统计分析的无重复试验数据提出统计分析方法,并利用 该方法和SAS软件对作物育种两种典型无重复试验设计的数据进行分析,展示对无重复试验应用线性混合 模型分析的程序和过程,为植物育种无重复试验分析的实现和广泛应用提供依据,进而达到改进无重复试 验在科学研究中应用成效的目的。解决植物育种无重复田间试验品系效应估计及其效应差异显著性统计测 验的问题。
Empirical Study on Analyzing Unreplicated Trials Data of Crop Breeding Based on Linear Mixed Model
Abstract;【Objective】The study provied a new method for analysing unreplicated datas ,and show how to use the PROC MIXED Program that provieded by SAS to analyse the unreplicated datas.【Method】Using linear mixed models to analyze unreplicated experimental data that come from crop breeding,the corresponding program and process for analyzing the data of unreplicated trials based on an international standard statistical software(SAS). The data characteristics of unreplicated trials and the drawbacks of classical analysis of variance being indicated,using the principle of covariance structures of linear mixed models to analyze unreplicated trial data , the information criteria of model fit was used for selecting the covariance structure model. .【Result】The empirical study showed that the linear mixed model analysis provided estimates and tests of line effects for the data of unreplicated plant breeding trials; The outcome in ranking and selection of the lines tested based on the estimates was different from that based on the observed value of the lines;The choice of covariance models had important impact on the results of analysis of unreplicated trials. 【Conclusion】The unreplicated experimental data could be analyzed using the PROC MIXED procedure in SAS based on the principle of the linear mixed model, so that the method can resolve the inferior comparable problems of unreplicated breeding lines and the problems of unreplicated trail data cannot be performed in significant test . Key words:Unreplicated trial; Model selection; Covariance structure; Information criterion

regularization paths for generalized linear models via coordinate descent

regularization paths for generalized linear models via coordinate descent
Jerome Friedman
Stanford University
Trevor Hastie
Stanford University
Rob Tibshirani
Stanford University
Abstract We develop fast algorithms for estimation of generalized linear models with convex penalties. The models include linear regression, two-class logistic regression, and multinomial regression problems while the penalties include 1 (the lasso), 2 (ridge regression) and mixtures of the two (the elastic net). The algorithms use cyclical coordinate descent, computed along a regularization path. The methods can handle large problems and can also deal efficiently with sparse features. In comparative timings we find that the new algorithms are considerably faster than competing methods.
JSS
Journal of Statistical Software
January 2010, Volume 33, Issue 1. /

An Introduction to Generalized Linear Model

An Introduction to Generalized Linear Model

AN INTRODUCTION TO GENERALIZED LINEAR MIXED MODELSStephen D.KachmanDepartment of Biometry,University of Nebraska–LincolnAbstractLinear mixed models provide a powerful means of predicting breeding values.However,for many traits of economic importance the assumptions of linear responses,constant variance,and normality are questionable.Generalized linear mixed modelsprovide a means of modeling these deviations from the usual linear mixed model.Thispaper will examine what constitutes a generalized linear mixed model,issues involvedin constructing a generalized linear mixed model,and the modifications necessary toconvert a linear mixed model program into a generalized linear mixed model program.1IntroductionGeneralized linear mixed models(GLMM)[1,2,3,6]have attracted considerable at-tention over the years.With the advent of SAS’s GLIMMIX macro[5],generalized linear mixed models have become available to a larger audience.However,in a typical breeding evaluation generic packages are too inefficient and implementations in FORTRAN or C are needed.In addition,GLMM’s pose additional challenges with some solutions heading for ±∞.The objective of this paper is to provide an introduction to generalized linear mixed models.In section2,I will discuss some of the deficiencies of a linear model.In sec-tion3,I will present the generalized linear mixed model.In section4,I will present the estimation equations for thefixed and random effects.In section5,I will present a set of estimating equations for the variance components.In section6,I will discuss some of the computational issues involved when these approaches are used in practice.2Mixed modelsIn this section I will discuss the linear mixed model and when the implied assumptions are not appropriate.A linear mixed model isy|u∼N(Xβ+Zu,R)where u∼N(0,G),X and Z are known design matrices,and the covariance matrices R and G may depend on a set of unknown variance components.The linear mixed model assumes that the relationship between the mean of the dependent variable y and thefixed and random effects can be modeled as a linear function,the variance is not afunction of the mean,and that the random effects follow a normal distribution.Any or all these assumptions may be violated for certain traits.A case where the assumption of linear relationships is questionable is pregnancy rate. Pregnancy is a zero/one trait,that is,at a given point an animal is either pregnant(1)or is not pregnant(0).For example,a change in management that is expected to increase pregnancy rate by.1in a herd with a current pregnancy rate of.5would be expected to have a smaller effect in a herd with a current pregnancy rate of.8;that is,a treatment,an environmental factor,or a sire would be expected to have a larger effect when the mean pregnancy rate is.5than when the pregnancy rate is.8.Another case where the assumption of a linear relationship is questionable is the anal-ysis of growth.Cattle,pigs,sheep,and mice have similar growth curves over time;that is,after a period of rapid growth they reach maturity and the growth rate is considerably slower.The relationship of time with weight is not linear,with time having a much larger effect when the animal is young and a very small effect when the animal is mature.The assumption of constant variance is also questionable for pregnancy rate.When the predicted pregnancy rate,µ,for a cow is.5the variance isµ(1−µ)=.25.If on the other hand the predicted pregnancy rate for a cow is.8the variance drops to.16.For some production traits the variance increases as mean level of production increases.The assumption of normality is also questionable for pregnancy rate.It is difficult to justify the assumption that the density function of a random variable which only takes on two values is similar to a continuous bell shaped curve with values ranging from−∞to +∞.A number of approaches have been taken to address the deficiencies of a linear mixed model.T ransformations have been used to stabilize the variance,to obtain a linear rela-tionship,and to normalize the distribution.However the transformation needed to stabilize the variance may not be the same transformation needed to obtain a linear relationship. For example a log transformation to stabilize the variance has the side effect that the model on the original scale is multiplicative.Linear and multiplicative adjustments are used to adjust to a common base and to account for heterogeneous variances.Multiple trait analysis can be used to account for heterogeneity of responses in different environ-ments.Separate variances can be estimated for different environmental groups where the environmental groups are based on the observed production.Afinal option is to ignore the deficiencies of the linear mixed model and proceed as if a linear model does hold.The above options have the appeal that they are relatively simple and cheap to im-plement.Given the robustness of the estimation procedures,they can be expected to produce reasonable results.However,these options sidestep the issue that the linear mixed model is incorrect.Specifically we have a set of estimation procedures which are based on a linear mixed model and manipulate the data to make itfit a linear mixed model. It seems more reasonable to start with an appropriate model for the data and use an es-timation procedure derived from that model.A generalized linear mixed model is a model which gives us extraflexibility in developing an appropriate model for the data[1].σFixedEffects h(η)InverseLink LinearPredictor GComponentsVarianceRandomEffectsβRMean Observations Figure 1:Symbolic representation of a generalized linear mixed model.3A Generalized Linear Mixed ModelIn this section I will present a formulation of a generalized linear mixed model.It differs from presentations such as [1]in that it focuses more on the inverse link function rather than the link function to model the relationship between the linear predictor and the conditional mean.Generalized linear mixed models also includes the nonlinear mixed models of [4].Figure 1provides a symbolic representation of a generalized linear mixed model.As in a linear mixed model,a generalized linear mixed model includes fixed effects,β,with (e.g.,management effect);random effects,u ∼N(0,G ),(e.g.,breeding values);design matrices X and Z ;and a vector of observations,y ,(e.g.,pregnancy status)for which the conditional distribution given the random effects has mean,µ(e.g.,mean pregnancy rate),and covariance matrix,R ,(e.g.,variance of pregnancy status is µ(1−µ)).In addition,a generalized linear mixed model includes a linear predictor,η,and a link and/or inverse link function.In addition,the conditional mean,µ,depends on the linear predictor through an inverse link function,h (·),and the covariance matrix,R ,depends on µthrough a variance function.For example,the mean pregnancy rate,µ,depends on the effect of management and the breeding value of the animal.The management effect and breeding value actadditively on a conceptual underlying scale.Their combined effect on the underlying scale is expressed as the linear predictor,η.The linear predictor is then transformed to the observed scale(i.e.,mean pregnancy rate)through an inverse link function,h(η).A typical transformation would beh(η)=eη1+eη.3.1Linear Predictor,ηAs with a linear mixed model,thefixed and random effects are combined to form a linear predictorη=Xβ+Zu.With a linear mixed model the model for the vector of observations y is obtained by adding a vector of residuals,e∼N(0,R),as followsy=η+e=Xβ+Zu+e.Equivalently,the residual variability can be modeled asy|u∼N(η,R).Unless y has a normal distribution,the formulation using e is clumsy.Therefore,a gen-eralized linear mixed model uses a second approach to model the residual variability.The relationship between the linear predictor and the vector of observations in a generalized linear mixed model is modeled asy|u∼(h(η),R)where the notation,y|u∼(h(η)R),specifies that the conditional distribution of y given u has mean,h(η),and variance,R.The conditional distribution of y given u will be referred to as the error distribution.Choice of whichfixed and random effects to include in the model will follow the same considerations as for a linear mixed model.It is important to note that the effect of the linear predictor is expressed through an inverse link function.Except for the identity link function,h(η)=η,the effect of a one unit change inηi will not correspond to a one unit change in the conditional mean;that is, predicted progeny difference will depend on the progeny’s environment through h(η).The relationship between the linear predictor and the mean response on the observed scale will be considered in more detail in the next section.3.2(Inverse)Link FunctionThe inverse link function is used to map the value of the linear predictor for observation i,ηi,to the conditional mean for observation i,µi.For many traits the inverse link function is one to one,that is bothµi andηi are scalars.For threshold models,µi is a t×1vector,T able1:Common link functions and variance functions for various distributionsDistribution Link Inverse Link v(µ)Normal Identityη1Binomial/n Logit eη/(1+eη)µ(1−µ)/nProbitΦ(η)Poisson Log eηµGamma Inverse1/ηµ2Log eηwhere t is the number of ordinal levels.For growth curve models,µi is an n i×1vector andηi is a p×1vector,where the animal is measured n i times and there are p growth curve parameters.For the linear mixed model,the inverse link function is the identity function h(ηi)=ηi. For zero/one traits a logit link functionηi=ln(µi/[1−µi])is often used,the corresponding.The logit link function,unlike the identity link function, inverse link function isµi=eηiηiwill always yield estimated means in the range of zero to one.However,the effect of a one unit change in the linear predictor is not constant.When the linear predictor is0 the corresponding mean is.5.Increasing the linear predictor by1to1increases the corresponding mean by.23.If the linear predictor is3,the corresponding mean is.95. Increasing the linear predictor by1to4increases the corresponding mean by only.03. For most univariate link functions,link and inverse link functions are increasing monotonic functions.In other words,an increase in the linear predictor results in an increase in the conditional mean,but not at a constant rate.Selection of inverse link functions is typically based on the error distribution.Table1 lists a number of common distributions along with their link functions.Other considera-tions include simplicity and the ability to interpret the results of the analysis.3.3Variance FunctionThe variance function is used to model non-systematic variability.Typically with a generalized linear model,residual variability arises from two sources.First,variability arises from the sampling distribution.For example,a Poisson random variable with mean µhas a variance ofµ.Second,additional variability,or over-dispersion,is often observed.Modeling variability due to the sampling distribution is straight forward.In Table1the variance functions for some common sampling distributions are given.Variability due to over-dispersion can be modeled in a number of ways.One approach is to scale the residual variability as var(y i|u)=φv(µi),whereφis an over-dispersion parameter.A second approach is to add a additional random effect,e i∼N(0,φ),to the linear predictor for each observation.A third approach is to select another distribution. For example,instead of using a one parameter(µ)Poisson distribution for count data,a two parameter(µ,φ)negative binomial distribution could be used.The three approaches all involve the estimation of an additional parameter,φ.Scaling the residual variability is the simplest approach,but can yield poor results.The addition of a random effect has theeffect of greatly increasing the computational costs.3.4The partsT o summarize,a generalized linear model is composed of three parts.First,a linear predictor,η=Xβ+Zu,is used to model the relationship between thefixed and random effects.The residual variability contained in the residual,e,of the linear mixed model equation is incorporated in the variance function of the generalized linear mixed model. Second,an inverse link function,µi=h(ηi),is used to model the relationship between the linear predictor and the conditional mean of the observed trait.In general,the link function is selected to be both simple and reasonable.Third,a variance function,v(µi,φ), is used to model the residual variability.Selection of the variance function is typically dictated by the error distribution that was chosen.In addition,observed residual variability is often greater than expected due to sampling and needs to be accounted for with an overdispersion parameter.4Estimation and PredictionThe estimating equations for a generalized linear mixed model can be derived in a number of ways.From a Bayesian perspective the solutions to the estimating equations are posterior mode predictors.The estimating equations can be obtained by using a Laplacian approximation of the likelihood.The estimating equations for thefixed and random effects areX H R−1HX X H R−1HZ Z H R−1HX Z H R−1HZ+G−1 βu= X H R−1y∗Z H R−1y∗(1)whereH=∂µ∂ηR=var(y|u)y∗=y−µ+Hη.Thefirst thing to note is the similarity to the usual mixed model equations.This is easier to see if we rewrite the equations in the following formX W X X W Z Z W X Z W Z+G−1 βu= X hryZ hry(2)whereW=H R−1Hhry=H R−1y∗.Unlike the mixed model equations,the estimating equations(1)for a generalized linear mixed must be solved iteratively.0.20.40.60.81-4-2024M e a n (µ)Linear Predictor (η)Figure 2:Inverse logit link function µ=e η/(1+e η).4.1Univariate LogitT o see how all the pieces fit together,we will examine a univariate binomial with a logit link function which would be an appropriate model for examining proportion data.Let y i be the proportion out of n for observation i ,a reasonable error distribution for ny i would be Binomial with parameters,n and µi .For a univariate model,the H ,R ,and W matricesare all diagonal with diagonal elements equal to ∂µi i ,ν(µi ),and W ii respectively.The variance function for a scaled Binomial random variable is ν(µi )=µi (1−µi )/n .The inverse link function is selected to model the nonlinear response of the means,µi ,tochanges in the linear predictor,ηi .The inverse link function selected is µi =e ηi 1+e ηi .The nonlinear relationship between the linear predictor and the mean can be seen in Figure 2.Changes in the linear predictor when the mean is close to .5have a larger impact than similar changes when the mean is close to 0or 1.For example,a change in the linear predictor from 0to .5changes the mean from 50%to 62%.However,a change from 2to2.5changes the mean from 88%to 92%.The weight,W ii,and scaled dependent variable,hry i,for observation i are∂µi ∂ηi =∂eηiηi∂ηi=µi(1−µi)W ii=∂µi∂ηi[ν(µi)]−1∂µi∂ηi=µi(1−µi) nµi(1−µi) µi(1−µi)=nµi(1−µi)hry i=∂µi∂ηi[ν(µi)]−1 y i−µi+∂µi∂ηiηi=n y i−µi+µi(1−µi)ηi .This can be translated into a FORTRAN subroutine as followsSUBROUTINE LINK(Y,WT,ETA,MU,W,HRY)REAL*8Y,WT,ETA,MU,W,R,VAR,Hmu=exp(eta)/(1.+exp(eta))!h(ηi)=eηi/(1+eηi)h=mu*(1.-mu)!∂µi/∂ηi=µi(1−µi)var=mu*(1.-mu)/wt!ν(µi)=µi(1−µi)/nW=(H/VAR)*H!W=Diag(H R−1H)HRY=(H/VAR)*(Y-MU)+W*ETA!hry=H R−1[y−µ+Hη]RETURNENDIn the subroutine Y=y i,WT=n,and ETA= ηi are passed to the subroutine.The sub-routine returns W=W ii,MU=µi,and HRY=hry i.In the subroutine the linesmu=exp(eta)/(1.+exp(eta))h=mu*(1.-mu)would need to be changed if a different link function was selected.The linevar=h/wtwould need to be changed if a different variance function was selected.The changes to the LINK subroutine to take into account boundary conditions will be discussed in sec-tion6.For each round of solving the estimating equations(2),a new set of linear predictors needs to be calculated.This can be accomplished during the construction of the LHS and RHS assuming solutions are not destroyed.The FORTRAN code for thisDO I=1,N!Loop to read in the N records.Read in recordETA=0!Calculate the risk factorηi=ETA based on the DO J=1,NEFF!solution for the NEFFfixed and random effects ETA=ETA+X*SOL!SOL and the design matrix X.END DOCALL LINK(Y,WT,ETA,MU,W,HRY)!Calculate µi=MU,W ii=W,and hry i=HRY!based on y i=Y,n=Wt,and ηi=ETA.Build LHS and RHSEND DO4.2Numerical exampleThe data in T able2have been constructed to illustrate the calculations involved for binomial random variables with a logit link function.A linear predictor of pregnancy rate for a cow in herd h on diet d isµhd=µ+Herd h+Diet dwhere Herd h is the effect of herd h and Diet d is the effect of diet d.For thefirst round we will use an initial estimate forηi of0.The contributions of each observation for round one is given in T able3.The estimating equations and solutions for round1are37.512.52515202.512.512.5057.50250251012.52.5155101500207.512.502002.502.5002.5µ Herd1Herd2Diet ADiet BDiet C=4053512253µ Herd1Herd2Diet ADiet BDiet C=1.2−1.04416−0.05194810.441558and η=X β=0.1038960.5974031.148051.641561.2.The new estimates of the linear predictor are used to obtain the contributions of each observation for round two given in T able4.Table2:Example data.Number PregnancyHerd Diet of Cows Rate1A20.501B30.662A40.802B50.902C10.80T able3:Contributions to the estimating equations for round one.Herd Diet y iηiµi n i W ii hry i1A0.500.520501B0.¯600.5307.552A0.800.54010122B0.900.55012.5202C0.800.510 2.53T able4:Contributions to the estimating equations for round two. Herd Diet y iηiµi n i W ii hry i1A0.50.1038960.52595120 4.98653-0.000932566 1B0.¯60.5974030.64506230 6.86871 4.751532A0.8 1.148050.759155407.3135510.03012B0.9 1.641560.83774750 6.7963514.26932C0.8 1.20.76852510 1.77894 2.449495Variance Component EstimationThe code given above assumes that the variance components are known.In this section I will discuss how estimates of the variance components can be obtained.Con-ceptually the variance component problem can be broken into two parts.Thefirst part is the estimation of the variance components associated with the random effects.The second part is the estimation of the variance components associated with the error dis-tribution.Before examining the modifications for a GLMM we will briefly review the major features of a variance component estimation program for a mixed model.5.1Mixed modelDerivative based programs for estimation of variance components under a mixed model involve the computation of quadratic forms of the predicted random effects, u Q i u, along with functions of the elements of a generalized inverse of the left hand sides,f ij(C), where C is a generalized inverse of the left hand sides in(2).For example,the univariate Fisher scoring quadratic form of the REML estimator of variance component i isu i I q i u iσiwhere u i is the q i×1vector of predicted random effects for i th set of random effects.The function of the left hand sides aref ii(C)=1σ4iq i−2tr(C ii)σ2i+tr(C ii C ii)σ4if ij(C)=1σiσjtr(C ij C ji)σiσjwhereC=C00C01 0C10C11 (1)...C r0C r1...C rris the partitioned generalized inverse of the left hand sides of(2).For the residual variance component,the quadratic form is(y−X β−Z u) I N(y−X β−Z u)σwhere N is the number of observations.The functions of the left hand sides aref i0(C)=1σ2iσ2tr(C ii)σ2i−rj=1tr(C ij C ji)σ2iσ2jf00(C)=1σ4N−p∗−q+ri=1rj=1tr(C ij C ji)σ2iσ2j.(3)5.2Modifications for GLMMThe variance components can be estimated using the approximate REML quasi-likelihood [1]which after some algebra isql (β,σ)=−12ln |V |−12ln |X H V −1HX |−12(y ∗−HXβ) V −1(y ∗−HXβ)(4)where σis the vector of variance component and V =R +HZGZ H .For the variance components associated with the random effects in G the estimating equations remain thesame except uand C are obtained using (2)instead of the usual mixed model equations.Estimation of the variance components associated with the error distribution is more problematic.The quadratic form becomes(y − µ) R −1∂R ∂φR −1(y − µ).(5)However,the corresponding functions for the left hand side in (3)for the linear mixedmodel assumes that R =I σ2o .The functions of the left hand sides for φaref 00(C )=[tr(Φ)−2tr(ΩΦ)+tr(ΩΦΩΦ)]f i 0(C )= tr(C i X H ΦHX X H ΦHZ Z H ΦHX Z H ΦHZC i ) where C i =(C i 0C i 1...C ir ),Φ=R −1∂R ∂φR −1andΩ= HX HZ CX H Z H .6Some Computational Issues While the mathematics are “straight forward,”implementation in practice is often chal-lenging.Many of the difficulties associated with generalized linear mixed models are related to either estimates going to infinity or divide by zero problems.Consider the uni-variate logit model.The calculations involve 1i i .Provided 0<µ<1this quantity is well defined.If µapproaches either zero or one,then the quantity 1µi (1−µi )approaches infinity.Estimates of µof zero or one occur when a contemporary group consists entirely of zeros or ones.Several options exist for handling this situation.One approach is to remove from the data any troublesome contemporary groups.While mathematically sound,an additional edit is needed to remove legitimate data values.A second approach is to treat contem-porary group effects as random effects.However,the decision to treat a set of effects asfixed or random should be decided from a modeling standpoint and not as an artifact of the estimation procedure.A third approach is to adjust y i away from0and1.For exam-ple,one could use(y i+∆)/(1+2∆)instead of y i.A fourth approach is based on the examination of what happens to the quantities W ii and hry i whenηi approaches infinity.Asηi→±∞,W ii→0andhry i→n(y i−µi)In the limit W ii and hry i are both well defined.One could then recode the link function as SUBROUTINE LINK(Y,WT,ETA,MU,W,HRY)REAL*8Y,WT,ETA,MU,W,R,VAR,Hif(abs(eta)>BIG)then!Check for ηi→±∞if(eta>0)then! ηi→∞⇒ µi→1mu=1.else! ηi→−∞⇒ µi→0mu=0end ifw=0! ηi→±∞⇒W ii→0hry=wt*(y-mu)!and hry i→n(y i− µi)returnend ifmu=exp(eta)/(1.+exp(eta))h=mu*(1.-mu)var=h/wtW=(H/VAR)*HHRY=(H/VAR)*(Y-MU)+W*ETARETURNENDwhere BIG is a sufficiently large,but not too large a number.For example BIG=10might be a reasonable choice.However,using W ii=0will usually result in a zero diagonal in the estimating equations.When these equations are solved,thefixed effect which was approaching infinity will be set to zero.This substitution can result in interesting convergence problems.A solution would be to“fix”the offensivefixed effect by making the following changesSUBROUTINE LINK(Y,WT,ETA,MU,W,HRY)REAL*8Y,WT,ETA,MU,W,R,VAR,Hif(abs(eta)>BIG)then!Check for ηi→±∞if(eta>0)then!−BIG≤ ηi≤BIGeta=BIGelseeta=-BIGend ifmu=exp(eta)/(1.+exp(eta))h=mu*(1.-mu)var=h/wtW=(H/VAR)*HHRY=(H/VAR)*(Y-MU)+W*ETARETURNEND6.1Additional ParametersThreshold models with more than two classes and growth curve models require more than one linear predictor per observation.In the case of a threshold model the value of an observation is determined by the usual linear predictor and a threshold linear predictor.With growth curve models animal i has n i observations.The linear predictor for animal i is a p ×1vector,where p is the the number of growth curve parameters.Unlike the univariate case,the matrices H ,W ,and R may not be diagonal.Typical structures for these matrices are:R =Diag(R i )with R i a n i ×n i residual covariance matrixH =Diag(H i )with H i a n i ×p matrix of partial derivatives H i = ∂µi 1∂ηi 1...∂µi 1∂ηip ...∂µin i i 1...∂µin i ip W =Diag(W i )with W i =H i R −1i H i .7ConclusionsGeneralized linear mixed models provide a flexible way to model production traits which do not satisfy the assumptions of a linear mixed model.This flexibility allows the researcher to focus more on selecting an appropriate model as opposed to finding manip-ulations to make the data fit a restricted class of models.The flexibility of a generalized linear mixed model provides an extra challenge when selecting an appropriate model.As a general rule the aim should be to select as simple a model as possible which does a reasonable job of modeling the data.While generalized linear mixed model programs are not as readily available as linear mixed model programs,modifications needed for a linear mixed model program should be minor.The two changes that are needed are to add a Link subroutine and to solve iteratively the estimating equations.In constructing a link subroutine it is important to handle boundary conditions robustly.When selecting a link function it is important to remember that when you get the results from the analysis,you will need to both interpret the results and to present the results in a meaningful manner.Generalized linear mixed models provide us with a very powerfulReferences[1]N.E.Breslow and D.G.Clayton.Approximate inference in generalized linear mixedmodels.J.Amer.Statist.Assoc.,88:9–25,1993.[2]B.Engel and A.Keen.A simple approach for the analysis of generalized linear mixedmodels.Statist.Neerlandica,48(1):1–22,1994.[3]D.Gianola and J.L.Foulley.Sire evaluation for ordered categorical data with a thresh-old model.G´en´et.S´el.Evol.,15(2):210–224,1983.[4]M.J.Lindstrom and D.M.Bates.Nonlinear mixed effects models for repeated mea-sures data.Biometrics,46:673–687,1990.[5]Ramon C.Littell,George liken,Walter W.Stroup,and Russell D.Wolfinger.SASSystem for Mixed Models.SAS Institute Inc.,Cary,NC,1996.[6]I.Misztal,D.Gianola,and puting aspects of a nonlinear method ofsire evaluation for categorical data.J.Dairy Sci.,72:1557–1568,1989.。

SAS_PROC_MIXED

SAS_PROC_MIXED

SAS PROC MIXED 1SAS PROC MIXEDhttp://www.id.unizh.ch/software/unix/statmath/sas/sasdoc/stat/chap41/index.htmOverviewThe MIXED procedure fits a variety of mixed linear models to data and enables you to use these fitted models to make statistical inferences about the data. A mixed linear model is a generalization of the standard linear model used in the GLM procedure, the generalization being that the data are permitted to exhibit correlation and nonconstant variability. The mixed linear model, therefore, provides you with the flexibility of modeling not only the means of your data (as in the standard linear model) but their variances and covariances as well. The primary assumptions underlying the analyses performed by PROC MIXED are as follows: • • • The data are normally distributed (Gaussian). The means (expected values) of the data are linear in terms of a certain set of parameters. The variances and covariances of the data are in terms of a different set of parameters, and they exhibit a structure matching one of those available in PROC MIXED.Since Gaussian data can be modeled entirely in terms of their means and variances/covariances, the two sets of parameters in a mixed linear model actually specify the complete probability distribution of the data. The parameters of the mean model are referred to as fixed-effects parameters, and the parameters of the variance-covariance model are referred to as covariance parameters. The fixed-effects parameters are associated with known explanatory variables, as in the standard linear model. These variables can be either qualitative (as in the traditional analysis of variance) or quantitative (as in standard linear regression). However, the covariance parameters are what distinguishes the mixed linear model from the standard linear model. The need for covariance parameters arises quite frequently in applications, the following being the two most typical scenarios: • • The experimental units on which the data are measured can be grouped into clusters, and the data from a common cluster are correlated. Repeated measurements are taken on the same experimental unit, and these repeated measurements are correlated or exhibit variability that changes.The first scenario can be generalized to include one set of clusters nested within another. For example, if students are the experimental unit, they can be clustered into classes, which in turn can be clustered into schools. Each level of this hierarchy can introduce an additional source of variability and correlation. The second scenario occurs in longitudinal studies, where repeated measurements are taken over time. Alternatively, the repeated measures could be spatial or multivariate in nature. PROC MIXED provides a variety of covariance structures to handle the previous two scenarios. The most common of these structures arises from the use of random-effects parameters, which are additional unknown random variables assumed to impact the variability of the data. The variances of the random-effects parameters, commonly known as variance components, become the covariance parameters for this particular structure. Traditional mixed linear models contain both fixed- and random-effects parameters, and, in fact, it is the combination of these two types of effects that led to the name mixed model. PROC MIXED fits not only these traditional variance component models but numerous other covariance structures as well. PROC MIXED fits the structure you select to the data using the method of restricted maximum likelihood (REML), also known as residual maximum likelihood. It is here that the Gaussian assumption for the data is exploited. OtherSAS PROC MIXED 2estimation methods are also available, including maximum likelihood and MIVQUE0. The details behind these estimation methods are discussed in subsequent sections. Once a model has been fit to your data, you can use it to draw statistical inferences via both the fixed-effects and covariance parameters. PROC MIXED computes several different statistics suitable for generating hypothesis tests and confidence intervals. The validity of these statistics depends upon the mean and variance-covariance model you select, so it is important to choose the model carefully. Some of the output from PROC MIXED helps you assess your model and compare it with others.Basic FeaturesPROC MIXED provides easy accessibility to numerous mixed linear models that are useful in many common statistical analyses. In the style of the GLM procedure, PROC MIXED fits the specified mixed linear model and produces appropriate statistics. Some basic features of PROC MIXED are • • • • • • • covariance structures, including variance components, compound symmetry, unstructured, AR(1), Toeplitz, spatial, general linear, and factor analytic GLM-type grammar, using MODEL, RANDOM, and REPEATED statements for model specification and CONTRAST, ESTIMATE, and LSMEANS statements for inferences appropriate standard errors for all specified estimable linear combinations of fixed and random effects, and corresponding t- and F-tests subject and group effects that enable blocking and heterogeneity, respectively REML and ML estimation methods implemented with a Newton-Raphson algorithm capacity to handle unbalanced data ability to create a SAS data set corresponding to any tablePROC MIXED uses the Output Delivery System (ODS), a SAS subsystem that provides capabilities for displaying and controlling the output from SAS procedures. ODS enables you to convert any of the output from PROC MIXED into a SAS data set. See the "Changes in Output" section.Notation for the Mixed ModelThis section introduces the mathematical notation used throughout this chapter to describe the mixed linear model. You should be familiar with basic matrix algebra (refer to Searle 1982). A more detailed description of the mixed model is contained in the "Mixed Models Theory" section. A statistical model is a mathematical description of how data are generated. The standard linear model, as used by the GLM procedure, is one of the most common statistical models:In this expression, y represents a vector of observed data, known design matrix X, andis an unknown vector of fixed-effects parameters with . Theis an unknown random error vector modeling the statistical noise aroundSAS PROC MIXED 3focus of the standard linear model is to model the mean of y by using the fixed-effects parameters . The residual errors are assumed to be independent and identically distributed Gaussian random variables with mean 0 and variance .The mixed model generalizes the standard linear model as follows:Here, is an unknown vector of random-effects parameters with known design matrix Z, and is an unknown random error vector whose elements are no longer required to be independent and homogeneous.To further develop this notion of variance modeling, assume that and are Gaussian random variables that are uncorrelated and have expectations 0 and variances G and R, respectively. The variance of y is thus V = ZGZ' + R Note that, when and Z = 0, the mixed model reduces to the standard linear model.You can model the variance of the data, y, by specifying the structure (or form) of Z, G, and R. The model matrix Z is set up in the same fashion as X, the model matrix for the fixed-effects parameters. For G and R, you must select some covariance structure. Possible covariance structures include • • • • • • • variance components compound symmetry (common covariance plus diagonal) unstructured (general covariance) autoregressive spatial general linear factor analyticBy appropriately defining the model matrices X and Z, as well as the covariance structure matrices G and R, you can perform numerous mixed model analyses.PROC MIXED Contrasted with Other SAS ProceduresPROC MIXED is a generalization of the GLM procedure in the sense that PROC GLM fits standard linear models, and PROC MIXED fits the wider class of mixed linear models. Both procedures have similar CLASS, MODEL, CONTRAST, ESTIMATE, and LSMEANS statements, but their RANDOM and REPEATED statements differ (see the following paragraphs). Both procedures use the nonfull-rank model parameterization, although the sorting of classification levels can differ between the two. PROC MIXED computes only Type I -Type III tests of fixed effects, while PROC GLM offers Types I - IV. The RANDOM statement in PROC MIXED incorporates random effects constituting the vector in the mixed model. However, in PROC GLM, effects specified in the RANDOM statement are still treated as fixed as far as the model fit is concerned, and they serve only to produce correspondingSAS PROC MIXED 4expected mean squares. These expected mean squares lead to the traditional ANOVA estimates of variance components. PROC MIXED computes REML and ML estimates of variance parameters, which are generally preferred to the ANOVA estimates (Searle 1988; Harville 1988; Searle, Casella, and McCulloch 1992). Optionally, PROC MIXED also computes MIVQUE0 estimates, which are similar to ANOVA estimates. The REPEATED statement in PROC MIXED is used to specify covariance structures for repeated measurements on subjects, while the REPEATED statement in PROC GLM is used to specify various transformations with which to conduct the traditional univariate or multivariate tests. In repeated measures situations, the mixed model approach used in PROC MIXED is more flexible and more widely applicable than either the univariate or multivariate approaches. In particular, the mixed model approach provides a larger class of covariance structures and a better mechanism for handling missing values (Wolfinger and Chang 1995). PROC MIXED subsumes the VARCOMP procedure. PROC MIXED provides a wide variety of covariance structures, while PROC VARCOMP estimates only simple random effects. PROC MIXED carries out several analyses that are absent in PROC VARCOMP, including the estimation and testing of linear combinations of fixed and random effects. The ARIMA and AUTOREG procedures provide more time series structures than PROC MIXED, although they do not fit variance component models. The CALIS procedure fits general covariance matrices, but it does not allow fixed effects as does PROC MIXED. The LATTICE and NESTED procedures fit special types of mixed linear models that can also be handled in PROC MIXED, although PROC MIXED may run slower because of its more general algorithm. The TSCSREG procedure analyzes time-series cross-sectional data, and it fits some structures not available in PROC MIXED.SyntaxThe following statements are available in PROC MIXED. PROC MIXED < options > ; BY variables ; CLASS variables ; ID variables ; MODEL dependent = < fixed-effects > < / options > ; RANDOM random-effects < / options > ; REPEATED < repeated-effect > < / options > ; PARMS (value-list) ... < / options > ; PRIOR < distribution > < / options > ; CONTRAST 'label' < fixed-effect values ... > < | random-effect values ... > , ... < / options > ; ESTIMATE 'label' < fixed-effect values ... > < | random-effect values ... >< / options > ; LSMEANS fixed-effects < / options > ; MAKE 'table' OUT=SAS-data-set ; WEIGHT variable ; Items within angle brackets ( < > ) are optional. The CONTRAST, ESTIMATE, LSMEANS, MAKE, and RANDOM statements can appear multiple times; all other statements can appear only once. The PROC MIXED and MODEL statements are required, and the MODEL statement must appear after the CLASS statement if a CLASS statement is included. The CONTRAST, ESTIMATE, LSMEANS, RANDOM, and REPEATED statements must follow the MODEL statement. The CONTRAST and ESTIMATE statements must also follow any RANDOM statements.SAS PROC MIXED 5Table 41.1 summarizes the basic functions and important options of each PROC MIXED statement. The syntax of each statement in Table 41.1 is described in the following sections in alphabetical order after the description of the PROC MIXED statement. Table 41.1: Summary of PROC MIXED Statements Statement PROC MIXED BY Description invokes the procedure performs multiple PROC MIXED analyses in one invocation declares qualitative variables that create indicator variables in design matrices lists additional variables to be included in predicted values tables Important Options DATA= specifies input data set, METHOD= specifies estimation method noneCLASSnoneIDnoneMODELspecifies dependent variable S requests solution for fixed-effects parameters, DDFM= specifies and fixed effects, setting up X denominator degrees of freedom method, OUTP= outputs predicted values to a data set specifies random effects, setting up Z and G sets up R SUBJECT= creates block-diagonality, TYPE= specifies covariance structure, S requests solution for random-effects parameters, G displays estimated G SUBJECT= creates block-diagonality, TYPE= specifies covariance structure, R displays estimated blocks of R, GROUP= enables between-subject heterogeneity, LOCAL adds a diagonal matrix to R HOLD= and NOITER hold the covariance parameters or their ratios constant, PDATA= reads the initial values from a SAS data set NSAMPLE= specifies the sample size, SEED= specifies the starting seed E displays the L matrix coefficients CL produces confidence limits DIFF computes differences of the least squares means, ADJUST= performs multiple comparisons adjustments, AT changes covariates, OM changes weighting, CL produces confidence limits, SLICE= tests simple effects none. Has been superceded by the Output Delivery System (ODS) noneRANDOMREPEATEDPARMSspecifies a grid of initial values for the covariance parameters performs a sampling-based Bayesian analysis for variance component models constructs custom hypothesis tests constructs custom scalar estimates computes least squares means for classification fixed effectsPRIORCONTRAST ESTIMATE LSMEANSMAKE WEIGHTconverts any displayed table into a SAS data set specifies a variable by which to weight RSAS PROC MIXED 6PROC MIXED StatementPROC MIXED < options >; The PROC MIXED statement invokes the procedure. You can specify the following options. ABSOLUTE makes the convergence criterion absolute. By default, it is relative (divided by the current objective function value). See the CONVF, CONVG, and CONVH options in this section for a description of various convergence criteria. ALPHA=number requests that confidence limits be constructed for the covariance parameter estimates with confidence level 1-number. The value of number must be between 0 and 1; the default is 0.05. ASYCORR produces the asymptotic correlation matrix of the covariance parameter estimates. It is computed from the corresponding asymptotic covariance matrix (see the description of the ASYCOV option, which follows). For ODS purposes, the label of the "Asymptotic Correlation" table is "AsyCorr." ASYCOV requests that the asymptotic covariance matrix of the covariance parameters be displayed. By default, this matrix is the observed inverse Fisher information matrix, which equals 2H-1, where H is the Hessian (second derivative) matrix of the objective function. See the "Covariance Parameter Estimates" section for more information about this matrix. When you use the SCORING= option and PROC MIXED converges without stopping the scoring algorithm, PROC MIXED uses the expected Hessian matrix to compute the covariance matrix instead of the observed Hessian. For ODS purposes, the label of the "Asymptotic Covariance" table is "AsyCov." CL<=WALD> requests confidence limits for the covariance parameter estimates. A Satterthwaite approximation is used to construct limits for all parameters that have a default lower boundary constraint of zero. These limits take the formwhere , Z is the Wald statistic , and the denominators are quantiles of the distribution with degrees of freedom. Refer to Milliken and Johnson (1992) and Burdick and Graybill (1992) for similar techniques. For all other parameters, Wald Z-scores and normal quantiles are used to construct the limits. The optional =WALD specification requests Wald limits for all parameters. The confidence limits are displayed as extra columns in the "Covariance Parameter Estimates" table. The confidence level is by default; this can be changed with the ALPHA= option.CONVF<=number> requests the relative function convergence criterion with tolerance number. The relative function convergence criterion isSAS PROC MIXED 7where fk is the value of the objective function at iteration k. To prevent the division by |fk|, use the ABSOLUTE option. The default convergence criterion is CONVH, and the default tolerance is 1E-8. CONVG <=number> requests the relative gradient convergence criterion with tolerance number. The relative gradient convergence criterion iswhere fk is the value of the objective function, and gjk is the jth element of the gradient (first derivative) of the objective function, both at iteration k. To prevent division by |fk|, use the ABSOLUTE option. The default convergence criterion is CONVH, and the default tolerance is 1E-8. CONVH<=number> requests the relative Hessian convergence criterion with tolerance number. The relative Hessian convergence criterion iswhere fk is the value of the objective function, gk is the gradient (first derivative) of the objective function, and Hk is the Hessian (second derivative) of the objective function, all at iteration k. If Hk is singular, then PROC MIXED uses the following relative criterion:To prevent the division by |fk|, use the ABSOLUTE option. The default convergence criterion is CONVH, and the default tolerance is 1E-8. COVTEST produces asymptotic standard errors and Wald Z-tests for the covariance parameter estimates.SAS PROC MIXED 8DATA=SAS-data-set names the SAS data set to be used by PROC MIXED. The default is the most recently created data set. DFBW has the same effect as the DDFM=BW option in the MODEL statement. EMPIRICAL computes the estimated variance-covariance matrix of the fixed-effects parameters by using the asymptotically consistent estimator described in Huber (1967), White (1980), Liang and Zeger (1986), and Diggle, Liang, and Zeger (1994). This estimator is commonly referred to as the "sandwich" estimator, and it is computed as follows:Here, , S is the number of subjects, and matrices with an i subscript are those for the ith subject. You must include the SUBJECT= option in either a RANDOM or REPEATED statement for this option to take effect. When you specify the EMPIRICAL option, PROC MIXED adjusts all standard errors and test statistics involving the fixed-effects parameters. This changes output in the following tables (listed in Table 41.7): Contrast, CorrB, CovB, Diffs, Estimates, InvCovB, LSMeans, MMEq, MMEqSol, Slices, SolutionF, Tests1 -Tests3. The OUTP= and OUTPM= data sets are also affected. Finally, the Satterthwaite and Kenward-Roger degrees of freedom methods are not available if you specify EMPIRICAL. IC displays a table of various information criteria. Four different criteria are computed in four different ways, producing 16 values in all. Table 41.2 displays the four criteria in both larger-is-better and smaller-is-better forms. Table 41.2: Information Criteria Criteria Larger-is-better Smaller-is-better Reference AIC l- d HQIC l- d loglogn BIC l- d/2 logn CAIC l- d(logn + 1)/2 -2l+ 2d -2l+ 2d loglogn -2l+ d logn -2l+ d(logn + 1) Akaike (1974) Hannan and Quinn (1979) Schwarz (1978) Bozdogan (1987)Here l denotes the maximum value of the (possibly restricted) log likelihood, d the dimension of the model, and n the number of effective observations. In Version 6 of SAS/STAT software, n equals the number of valid observations for maximum likelihood estimation and n-p for restricted maximum likelihood estimation, where p equals the rank of X. In later versions, n equals the number of effective subjects as displayed in the "Dimensions" table, unless this value equals 1, in which case n reverts to the Version 6 values. PROC MIXED evaluates the criteria for both forms using d equal to both q and q+p, where q is theSAS PROC MIXED 9effective number of estimated covariance parameters. In Version 6, when a parameter estimate lies on a boundary constraint, then it is still included in the calculation of d, but in later versions it is not. The most common example of this behavior is when a variance component is estimated to equal zero. For ODS purposes, the name of the "Information Criteria" table is "InfoCrit." INFO is a default option. The creation of the "Model Information" and "Dimensions" tables can be suppressed using the NOINFO option. Note that, in Version 6, this option displays the "Model Information" and "Dimensions" tables. ITDETAILS displays the parameter values at each iteration and enables the writing of notes to the SAS log pertaining to "infinite likelihood" and "singularities" during Newton-Raphson iterations. LOGNOTE writes periodic notes to the log describing the current status of computations. It is designed for use with analyses requiring extensive CPU resources. MAXFUNC=number specifies the maximum number of likelihood evaluations in the optimization process. The default is 150. MAXITER=number specifies the maximum number of iterations. The default is 50. METHOD=REML METHOD=ML METHOD=MIVQUE0 METHOD=TYPE1 METHOD=TYPE2 METHOD=TYPE3 specifies the estimation method for the covariance parameters. The REML specification performs residual (restricted) maximum likelihood, and it is the default method. The ML specification performs maximum likelihood, and the MIVQUE0 specification performs minimum variance quadratic unbiased estimation of the covariance parameters. The METHOD=TYPEn specifications apply only to variance component models with no SUBJECT= effects and no REPEATED statement. An analysis of variance table is included in the output, and the expected mean squares are used to estimate the variance components (refer to Chapter 30, "The GLM Procedure," for further explanation). The resulting method-of-moment variance component estimates are used in subsequent calculations, including standard errors computed from ESTIMATE and LSMEANS statements. For ODS purposes, the new table names are "Type1," "Type2," and "Type3," respectively. MMEQ requests that coefficients of the mixed model equations be displayed. These areassuming thatis nonsingular. Ifis singular, PROC MIXED produces the following coefficientsSAS PROC MIXED 10See the "Estimating and in the Mixed Model" section for further information on these equations. MMEQSOL requests that a solution to the mixed model equations be produced, as well as the inverted coefficients matrix. Formulas for these equations are provided in the preceding description of the MMEQ option. When is singular, and a generalized inverse of the left-hand-side coefficient matrix are transformed is a generalized inverse of the left-hand-sideto produce and , respectively, where using coefficient matrix of the original equations.NAMELEN<=number> specifies the length to which long effect names are shortened. The default and minimum value is 20. NOBOUND has the same effect as the NOBOUND option in the PARMS statement. NOCLPRINT<=number> suppresses the display of the "Class Level Information" table if you do not specify number. If you do specify number, only levels with totals that are less than number are listed in the table. NOINFO suppresses the display of the "Model Information" and "Dimensions" tables. NOITPRINT suppresses the display of the "Iteration History" table. NOPROFILE includes the residual variance as part of the Newton-Raphson iterations. This option applies only to models that have a residual variance parameter. By default, this parameter is profiled out of the likelihood calculations, except when you have specified the HOLD= or NOITER option in the PARMS statement. ORD displays ordinates of the relevant distribution in addition to p-values. The ordinate can be viewed as an approximate odds ratio of hypothesis probabilities. ORDER=DATA ORDER=FORMATTED ORDER=FREQ ORDER=INTERNAL specifies the sorting order for the levels of all CLASS variables. This ordering determines which parameters in the model correspond to each level in the data, so the ORDER= option may be useful when you use CONTRAST or ESTIMATE statements. The default is ORDER=FORMATTED, and its behavior has been modified for Version 8. Now, for numeric variables for which you have supplied no explicit format (that is, for which there is noSAS PROC MIXED 11corresponding FORMAT statement in the current PROC MIXED run or in the DATA step that created the data set), the levels are ordered by their internal (numeric) value. In releases previous to Version 8, numeric class levels with no explicit format were ordered by their BEST12. formatted values. In order to revert to the previous method you can specify this format explicitly for the CLASS variables. The change was implemented because the former default behavior for ORDER=FORMATTED often resulted in levels not being ordered numerically and required you to use an explicit format or ORDER=INTERNAL to get the more natural ordering. The following table shows how PROC MIXED interprets values of the ORDER= option. Value of ORDER= Levels Sorted By DATA FORMATTED order of appearance in the input data set external formatted value, except for numeric variables with no explicit format, which are sorted by their unformatted (internal) value FREQ descending frequency count; levels with the most observations come first in the order INTERNAL unformatted valueFor FORMATTED and INTERNAL, the sort order is machine dependent. For more information on sorting order, see the chapter on the SORT procedure in the SAS Procedures Guide and the discussion of BYgroup processing in SAS Language Reference: Concepts. RATIO produces the ratio of the covariance parameter estimates to the estimate of the residual variance when the latter exists in the model. RIDGE=number specifies the starting value for the minimum ridge value used in the Newton-Raphson algorithm. The default is 0.3125. SCORING<=number> requests that Fisher scoring be used in association with the estimation method up to iteration number, which is 0 by default. When you use the SCORING= option and PROC MIXED converges without stopping the scoring algorithm, PROC MIXED uses the expected Hessian matrix to compute approximate standard errors for the covariance parameters instead of the observed Hessian. The output from the ASYCOV and ASYCORR options is similarly adjusted. SIGITER is an alias for the NOPROFILE option. UPDATE is an alias for the LOGNOTE option.BY StatementBY variables ;SAS PROC MIXED 12You can specify a BY statement with PROC MIXED to obtain separate analyses on observations in groups defined by the BY variables. When a BY statement appears, the procedure expects the input data set to be sorted in order of the BY variables. The variables are one or more variables in the input data set. If your input data set is not sorted in ascending order, use one of the following alternatives: • • Sort the data using the SORT procedure with a similar BY statement. Specify the BY statement options NOTSORTED or DESCENDING in the BY statement for the MIXED procedure. The NOTSORTED option does not mean that the data are unsorted but rather means that the data are arranged in groups (according to values of the BY variables) and that these groups are not necessarily in alphabetical or increasing numeric order. Create an index on the BY variables using the DATASETS procedure (in base SAS software).•Since sorting the data changes the order in which PROC MIXED reads observations, the sorting order for the levels of the CLASS variable may be affected if you have specified ORDER=DATA in the PROC MIXED statement. This, in turn, affects specifications in the CONTRAST statements. For more information on the BY statement, refer to the discussion in SAS Language Reference: Concepts. For more information on the DATASETS procedure, refer to the discussion in the SAS Procedures Guide.CLASS StatementCLASS variables ;The CLASS statement names the classification variables to be used in the analysis. If the CLASS statement is used, it must appear before the MODEL statement. Classification variables can be either character or numeric. The procedure uses only the first 16 characters of a character variable. Class levels are determined from the formatted values of the CLASS variables. Thus, you can use formats to group values into levels. Refer to the discussion of the FORMAT procedure in the SAS Procedures Guide and to the discussions of the FORMAT statement and SAS formats in SAS Language Reference: Dictionary. You can adjust the display order of CLASS variable levels with the ORDER= option in the PROC MIXED statement.CONTRAST StatementCONTRAST 'label' < fixed-effect values ...> < | random-effect values ...> , ...< / options > ; The CONTRAST statement provides a mechanism for obtaining custom hypothesis tests. It is patterned after the CONTRAST statement in PROC GLM, although it has been extended to include random effects. This enables you to select an appropriate inference space (McLean, Sanders, and Stroup 1991).You can test the hypothesis , where L' = (K' M') and , in several inference spaces. The inference space corresponds to the choice of M. When M = 0, your inferences apply to the entire population from which the random effects are sampled; this is known as the broad inference space. When all elements of M are nonzero, your inferences apply only to the observed levels of the random effects. This is known as the narrow inference space, and you can also choose it by specifying all of the random effects as fixed. The GLM procedure uses the narrow inference space. Finally, by zeroing portions of M corresponding to selected main effects and interactions, you can choose intermediate inference spaces. The broad inference space is usually the most appropriate, and it is used when you do not specify any random effects in the CONTRAST statement. In the CONTRAST statement, label。

Finite element model updating in structural dynamics by using the response surface method

Finite element model updating in structural dynamics by using the response surface method

Engineering Structures32(2010)2455–2465Contents lists available at ScienceDirectEngineering Structures journal homepage:/locate/engstructFinite element model updating in structural dynamics by using the response surface methodWei-Xin Ren∗,Hua-Bing ChenDepartment of Civil Engineering,Central South University,Changsha,410075,ChinaNational Engineering Laboratory for High Speed Railway Construction,Changsha,410075,Chinaa r t i c l e i n f o Article history:Received3February2009 Received in revised form6April2010Accepted7April2010 Available online4May2010 Keywords:Response surfaceFinite elementModel updatingOptimizationDesign of experiment Regression analysis a b s t r a c tFast-running response surface models that approximate multivariate input/output relationships of time-consuming physical-based computer models enable effective finite element(FE)model updating analyses. In this paper,a response surface-based FE model updating procedure for civil engineering structures in structural dynamics is presented.The key issues to implement such a model updating are discussed such as sampling with design of experiments,selecting the significant updating parameters and constructing a quadratic polynomial response surface.The objective function is formed by the residuals between analytical and measured natural frequencies.Single-objective optimization with equal weights of natural frequency residual of each mode is used for optimization computation.The proposed procedure is illustrated by a simulated simply supported beam and a full-size precast continuous box girder bridge tested under operational vibration conditions.The results have been compared with those obtained from the traditional sensitivity-based FE model updating method.The real application to a full-size bridge has demonstrated that the FE model updating process is efficient and converges fast with the response surface to replace the original FE model.With the response surface at hand,an optimization problem is formulated explicitly.Hence,no FE calculation is required in each optimization iteration.The response surface-based FE model updating can be easily implemented in practice with available commercial FE analysis packages.©2010Elsevier Ltd.All rights reserved.1.IntroductionNowadays the finite element(FE)method has become an important and practical numerical analysis tool.It is commonly used in almost all areas of engineering analysis.However,the FE model of a structure is normally constructed on the basis of highly idealized engineering blueprints and designs that may not truly represent all the aspects of an actual structure.As a result, the analytical predictions from a FE model often differ from the results of a real structure.These discrepancies originate from the uncertainties in simplifying assumptions of structural geometry, materials as well as inaccurate boundary conditions.It is often required to update or calibrate the uncertain parameters of a FE model that leads to the better predictions of the responses of an actual structural.Finite element model updating is such a procedure that mod-ifies or updates the uncertainty parameters in the initial finite element model based on the experimental results so that a more realistic or refined model can be achieved[1].In other words,FE ∗Corresponding author at:Department of Civil Engineering,Central South University,Changsha,410075,China.Tel.:+8673182654349;fax:+86731 85571736.E-mail addresses:renwx@,renwx@(W.-X.Ren).model updating is the process of using experimental results to re-fine a mathematical model of a physical structure.Basically,FE model updating is an inverse problem to identify or correct the uncertain parameters of FE models.It is usually posed as an op-timization problem.Setting-up of an objective function,selecting updating parameters and using robust optimization algorithm are the three crucial steps in FE model updating.In a model updating process,not only the satisfactory correlation is required between analytical and experimental results,but also the updated param-eters should preserve the physical significance.The updated FE models are used in many applications for civil engineering struc-tures such as damage detection,health monitoring,structural con-trol,structural evaluation and assessment.Finite element model updating is a topic of significant interest in the field of structural dynamics.A number of FE model updating methods in structural dynamics have been proposed.The direct updating methods compute a closed-form solution for the global stiffness and/or mass matrices using the structural equations of motion and the orthogonality equations.These non-iterative methods that directly update the elements of stiffness and mass matrices are one-step procedures[2,3].The resulting updated matrices reproduce the measured structural modal properties well but do not generally maintain structural connectivity and the corrections suggested are not always physically meaningful.0141-0296/$–see front matter©2010Elsevier Ltd.All rights reserved. doi:10.1016/j.engstruct.2010.04.0192456W.-X.Ren,H.-B.Chen/Engineering Structures32(2010)2455–2465The iterative parameter updating methods involve using the sensitivity of the parameters to find their changes.The sensitivity-based FE model updating methods are often posed as optimization problems.These methods set the errors of the structural response features between analytical and experimental results as an objective function and try to minimize the objective function by making changes to the pre-selected set of physical parameters of the FE model.Link[4]gave a clear overview of the sensitivity-based updating methods in structural dynamics.It is noted that the mathematical model used in the model updating is usually ill posed and the special attention is required for an accurate solution[5].Jaishi and Ren[6–8]used either single-objective or multi-objective optimization technique to update the FE models of civil engineering structures in structural dynamics.However,the sensitivity-based method involving in the determination of local gradients at points may cause not only computational intensive, but also convergence difficulty.If the structure of interest is represented by,e.g.a large finite element model,the large number of computations involved can rule out many approaches due to the expense of carrying out many runs.For such a large FE model where so many elements are involved,there are huge of both geometric and physical possible parameters to be updated.In addition,there are now many commercial finite element analysis packages available such as ANSYS,ABAQUS and SAP2000et al..The structural FE models are often constructed by using these packages.In all the iterative parameter updating methods,each iteration needs to go back to run the finite element analysis package with any parameter updated,which limits the popular applications of structural FE model updating in practice.One way of circumnavigating the time-consuming and FE analysis package-related problems during the sensitivity-based model updating is to replace the FE model by an approximate surrogate/replacement meta-model that is fast-running and less parameters involved[9].Such a meta-model is easier to compute with,because it is controlled only by a few explanatory variables. The FE model updating is carried out on the meta-model instead of the analytical FE model.Response surface is one of the commonly used meta-models.Response surface methodology is originally an experimental design approach for selecting design parameters for experiments with the objective of optimizing some function of a response[10–12].It provides a mechanism for guiding experimentation in search of optimal settings for design parameters or optimal values of unknown response.Many additional applications(largely a consequence of the increased use of computational analyses)have broadened the range of application of response surface methods in the statistical and engineering literature.Recent literature has addressed more flexible functional forms for modeling the response,new methods of response surface construction[13],alternate approaches to updating the surface estimate[14],new sampling methods[15], etc.In many fields of engineering,the term‘‘response surface’’is used synonymously with‘‘meta-model’’or‘‘surrogate model’’, which refer to any relatively simple mathematical relationship between parameters and a response,often based on limited data[16].In the case of structural finite element model updating, once the response surface of a structure has been constructed, updating the model is reduced to the task of finding the smallest value on the response surface.The parameter values that correspond to this smallest value are those that are used to update the model.Recently,the response surface that is represented by a simple least-squares multinomial model has been adopted in structural FE model updating,verification and validation[17–20].The response surface method for damage detection and reliability analysis is not quite new[21,22].However,the response surface method for structural finite element model updating is somewhat new,especially with the civil engineering communities. This paper is intended to present a response surface-based finite element model updating procedure in structural dynamics and to take advantages of response surfaces for the FE model updating of civil engineering structures in practice.Its purpose is to estimate the values of structural parameters(moment of inertia,cross-sectional area and modulus of elasticity)based on the measured response quantities(natural frequencies).The proposed procedure is based on constructing quadratic response surfaces.Those surfaces represent an estimate for the relation between the unknown parameters of the finite element model and response quantities of interest.With the response surface at hand,an optimization problem,whose solution is the estimate for the values of the structural parameters,is formulated explicitly. Hence,no further finite element simulations are required.The objective function is formed by the residuals between analytical and measured natural frequencies.Single-objective optimization with equal weights of each natural frequency is implemented for optimization computation.The presented procedure is illustrated by a simulated simply supported beam and a full-size precast continuous box girder bridge tested under operational vibration conditions.The results have been compared with those obtained from the traditional sensitivity-based FE model updating method. The real application to a full-size bridge has demonstrated that the response surface-based FE model updating procedure is simple and fast so that it can be easily implemented in practice.2.Response surface-based finite element model updatingResponse surface-based finite element model updating is an approach to achieve the global approximations of the structural response feature objectives and constrains based on functional evaluations at various points in the design space.It often involves experimental strategies,mathematical methods,probability and statistical inference that enable an experimenter to make efficient empirical exploration of the structure of interest.The flowchart of response surface-based finite element model updating is shown in Fig.1.The main steps include the following.•The selection of the structural parameters and the definition ofa number of‘‘level’’for each selected parameters by using thedesign of experiments(DOE)techniques.•In design space,the response features are calculated by carrying out finite element analysis(FEA).•Performing the final regression followed by a regression error analysis to create the response surface model of the structure.•The response features of the structure are measured and corresponding objective functions(feature residuals)to be minimized are constructed.•The finite element model updating(iteration)is carried out within the established response surface model.Updated parameters are obtained and transferred to the original finite element model.2.1.Sampling and parameter selectionTo create a response surface that will serve as a surrogate for the FE simulation model,the basic process is one of calculating predicted values of the response features at various sample points in the parameter space by performing a experiment at each of those points.A number of feature values from the experiment ran across the parameter domain are fit with a response surface. The key is to select the parameters carefully,to minimize the number of dimensions in the parameter space,and then to selectW.-X.Ren,H.-B.Chen /Engineering Structures 32(2010)2455–24652457Fig.1.Flowchart of response surface-based FE model updating.the combinations of parameter values where the experiment is performed.The term experiment herein refers to either physical experi-ments or computer experiments.The planning of experimentation is further referred to design of experiments (DOE).The selection of sample points is related to the accuracy and cost of a response sur-face to be constructed.Less sample points may reduce the surface accuracy,while more sample points may improve the surface accu-racy but increase the work load.In the real application,the sample points mainly depend on the problem to be solved,the response feature values of interest and the selected method of DOE.DOE plays an important role in constructing a response surface.Two efficient DOE methods are commonly used.They are central composite design (CCD)method and D-optimal design method.There are several other DOE methods that can be used in constructing a response surface such as box-behnken design,Latin Square design,fractional factor design [12].For the purpose of structural finite element model updating in structural dynamics,Guo and Zhang [18]found that both CCD and D-optimal design methods can achieve almost the same accuracy in the creation of polynomial surfaces.In the current study,the CCD method in DOE is used in the paper as it is simple in constructing the response surfaces of a polynomial type.Central composite design uses the orthogonal table to perform the experimentation to determine the sample points of selected parameters.It contains a fractional factorial design 2k (levels are ±1and k is quantity of factors)with central points that are augmented with a group of 2k star points that allow for the estimation of curvature.2k star points are (±α,0,...,0),(0,±α,...,0),...,(0,0,...,±α).For the purpose of constructing high-order surfaces,another 2k star points (±α1,0,...,0),(0,±α1,...,0),...,(0,0,...,±α1)are appended.The precise values of αand α1rely on certain properties such as orthogonality and rotatability desired for the design and on the number of factors involved [12].The selection of updating parameters to calculate the response features of a structure is an important issue in FE model updating.The structural parameters selected for updating should be able to clarify the ambiguity of the model,and in that case it is necessary for the model response to be sensitive to these parameters.The problem always arises:How many parameters should be selected and which parameters from many possible parameters are used in the FE model updating?If too many parameters are included in the FE model updating in structural dynamics,the optimized problem may appear ill-conditioned because only limited vibration modes can be effectively identified from the field vibration measurements.The parameter selection requires a considerable physical insight into the target structure,and trial-and-error approaches are often used with different set of selected parameters for complicated structures.In the frame of response surface-based FE model updating,the selected updating parameters can be evaluated from the sampled data by performing a parameter effect analysis (hypothesis testing)based on statistical variance (square of the standard deviation)analysis.Compared with tests on means for a hypothesis testing,tests on variances are rather sensitive to the normality assumption.Variance analysis formulates all individual square of deviations of sample data.The basic idea of variance analysis is to divide the total square deviation of sampling features into two parts:S A —a square of deviation caused by design parameter A (system deviation)and S e —a square of deviation caused by the experiment.The F-test method [12]can be used to carry out the hypothesis testing to check the significance of parameter A:F A =S A /f A S e /f e∼F (f A ,f e )(1)where f A and f e are the degrees of freedom of S A and S e respectively.For a given significance level p ,if F A ≥F 1−p (f a ,f e ),the effect of parameter A is significant.2.2.Response surface regressionThe family of response surface forms selected to represent a response can have a substantial impact on the results of an analysis.The selected response surface form should be capable of attaining surfaces that meet specific smoothness requirements of an application.There is often a balance between assumptions and data requirements.Different surface families may be preferred for different applications.In the case of structural finite element model updating in structural dynamics,polynomials are popular forms representing a response surface because the calculations are simple and the resulting function is closed-form algebraic expression.For example,a quadratic polynomial response surface has the form:y =β0+k i =1βi x i +k i =1k j =1βij x i x j(2)where β0,βi ,βij are the regression coefficients to be estimatedfrom the experimental data.In such a way,a response surface y is represented as a function of the parameters or variables x of the design space (generally some subset of Euclidian k -space,R k ).In the case of k =2,the vector of surface parameters is six dimensional.The number of sampling points n must be greater than or equal to the number of terms in the polynomial.When there are more sampling points,the equation is over-determined and the regression techniques are required to fit the response surface to the sampled data.In this case,the surface does not,in general,exactly match the response values at every sample points.The method2458W.-X.Ren,H.-B.Chen/Engineering Structures32(2010)2455–2465 Table1of least-square fitting is usually used in the coefficient estimationprocess to create a response surface.Before the regressed response surface is put forward to be usedin structural FE model updating,it should be verified to checkwhether the regressed surface has enough accuracy.If not,theparameters are adjusted(based on the response data)to achievea balance between variability and bias of the regressed surface.R2(ranged from0.0to1.0)criterion and empirical integrated squaredError(EISE)criterion can be used in response surface verification:R2=1−Ni=1[y RS(j)−y(j)]2Nj=1[y(j)−¯y]2(3)EISE=1N∗Nj=1[y RS(j)−y(j)]2(4)where y RS is the response value of the confirmation samples;y is the true value of the confirmation samples;¯y is the mean of all true values.The larger the value of R2,the more accurate the regressed response surface.On the contrary,the smaller the value of EISE,the closer the fit is to the data.Another criterion is the root mean squared error(RMSE).It is the square root of the mean square error(MSE)where the distance vertically from the point is taken to the corresponding point on the curve fit(the error).The squaring is done so negative values do not cancel positive values.The RMSE is thus the distance,on average,of a data point from the fitted line,measured along a vertical line.The smaller the value of root mean squared error,with 0corresponding to the ideal,the closer the fit is to the data.In the real application,the R2criterion and EISE or the RMSE criterion can be used in a complementary way to check the accuracy of the regressed response surface.2.3.Model updatingIn the finite element model updating in structural dynamics,the structural response features of interest are often eigen solutions related to such as natural frequencies and mode shapes.In this study,the structural natural frequency is employed as a response feature.Therefore,the optimized objective function is formulated in terms of the residuals between analytical and measured natural frequenciesΠ(x)=mi=1w i(λai−λei)0≤w i≤1(5)whereλai andλei are the analytical(finite element calculated) and measured natural frequencies of the i th mode respectively. w i is the weight factor to impose to the different order of naturalfrequencies.m is the number of modes involved and x is the design set.The FE model updating can then be posed as a constrained minimization problem to find the satisfied design set such that Min Π(x) 22(6)subjected tox lk≤x k≤x uk,k=1,2,3,...where x lk and x uk are the upper and lower bounds on the k th design variable to be set.N is the total number of design variables (updated parameters).Single-objective optimization with equal weight of each natural frequency is used in this paper for optimization computation.3.Numerical verificationA simulated simply supported beam without damage and with an assumed damage element is considered to demonstrate the procedure of response surface-based FE model updating.The beam of6m length is equally divided into15two-dimensional beam elements as shown in Fig.2.The density and elastic modulus of the beam are2500kg/m3and3.2×1010Pa,respectively,while the area and moment of inertia of cross-section are0.05m2and 1.66×10−4m4,respectively.FE modal analysis is first carried out on the beam without damage to get the analytical natural frequencies.To simulate the actual(target)beam,one damage location is assumed at beam element10where the element bending stiffness is reduced by50%. The FE modal analysis is again carried out on this damaged beam to get the simulated measured modal parameters.The initial values of the first10natural frequencies selected and corresponding differences are shown in Table1.The maximum and minimum errors that appeared in the natural frequency are15.63%and1.89% respectively.Updating of the FE model of the undamaged beam is to achieve a goal to correlate the natural frequencies with the damaged beam. In this study,moment of inertia I10and elastic modulus E10of individual element10are chosen as the updating parameters that affect the independent variable,natural frequencies.The range selected for each parameter should reflect the change that one expects to observe for the domain of the prediction of interest with 1.920≤E10≤3.465(×1010Pa)and0.863≤I10≤1.797 (×10−4m4).The central composite design(CCD)method in design of experiments is then used to get the sample points of the parameter values selected.A total of9runs of experiments are carried out.The sampled parameter values and corresponding natural frequencies calculated from FE model are listed in Table2.By using least-square fitting with these sample values,a quadratic polynomial response surface(Eq.(2))can be regressed where parameter x1refers to E10while parameter x2refers to I10.Fig.3illustrates the regressed response surfaces of the first and second natural frequencies.The horizontal axes are parameters selected while the vertical axis gives the response(natural frequency)at any point or location.To evaluate the above-selected parameters of the FE model,a statistics-based parameter effects analysis is performed.Parts of the general regression significance of selected parameters(x1and x2)calculated from Eq.(1)on each mode of natural frequencies comparing with the significance level of p=0.05are shown in Fig.4.The results show that both elastic modulus E10and moment of inertia I10are high significance on the structural natural frequencies.It is observed that the cross-quadratic term x1x2inW.-X.Ren,H.-B.Chen /Engineering Structures 32(2010)2455–24652459Fig.2.Simulated simply supported beam.Table 2Central composite design and sample values.9.598.58F r e q u e n c y7.576.53.532.5x 1×1010x 2×10–4221.51.510.5F r e q u e n c yx 1×1011.522.533.5x 2×10–40.511.523836343230(a)Surface of first natural frequency.(b)Surface of second natural frequency.Fig.3.Typical regressed response surfaces.Table 3the surface expression is less significant on the structural natural frequencies.For the simplification of a quadratic polynomial response surface,the cross-quadratic term can be omitted in practice to increase the efficiency of structural FE model updating.To check the accuracy of the regressed surface,the RMSE is calculated for each mode as shown in Table 3.It is demonstrated that all RMSE values are close to zero,which indicates that the created response surface has a high regression accuracy.Now it is time to carry out the FE model updating where the FE model is replaced by the regressed quadratic polynomial.The residuals between analytical (undamaged beam)and measured (damaged beam)natural frequencies are used as the optimized objective function.Single-objective optimization algorithm with equal weight of each natural frequency is implemented to achieve the best minimization of natural frequency residuals.The optimization algorithm used to minimize the objective function is a standard Trust Region Newton method in MATLAB.The tuning minimization process is over when the tolerances are achieved or pre-defined number of iterations is reached.The updated natural frequencies and differences of the simulated simply supported beam are also summarized in Table 1.It can be observed that a good agreement of natural frequencies has been achieved after carrying out the response surface-based FE model updating.The2460W.-X.Ren,H.-B.Chen/Engineering Structures32(2010)2455–2465(a)Parameter significance on the first natural frequency.(b)Parameter significance on the second natural frequency.Fig.4.The regression significance of selected parameters on the naturalfrequencies.Fig.5.Convergence of objective function of response surface-based updating.final updated results on the parameters are E10=1.57×1010Pa(true value E10=1.60×1010Pa)and I10=0.85×10−4m4(truevalue I10=0.83×10−4m4).Fig.5illustrates the convergence of the objective functionduring each optimization iteration.The horizontal axis is thenumber of iteration while the vertical axis is the square sum ofnatural frequency residuals of each mode.To demonstrate theefficiency of current response surface-based FE model updating,the objective function convergence of the traditional sensitivity-based FE model updating is shown in Fig.6.For such a numericalsimple beam,it is shown that the convergence of the objectivefunction is fast and iteration number is dramatically reduced toreach the same residual level of the objective function by using theresponse surface method.4.A precast concrete bridge tested in the fieldA real case study is carried out on the Hongtang Bridge,locatedin Fuzhou city,the capital of Fujian province,China.The bridge is amulti-span continues-deck precast concrete highway bridge.Theconstruction was completed in December1990.The total lengthof the bridge is1843m,with a span layout of(16+27+4∗30+60+120+60+31∗40+8∗25)m.It includesthreeFig.6.Convergence of objective function of sensitivity-basedupdating.Fig.7.Potion of Hongtang bridge.types of spans:simply supported spans,precast truss supportedspans and precast continuous girder spans.The photograph ofthe part of bridge at present condition is shown in Fig.7.Forthe purpose of dynamics-based condition assessment,one portionof six continuous girder spans with a span length of40m weretested on the site under operational vibration conditions.On-sitedynamic testing of a structure provides an accurate and reliabledescription of its current dynamic characteristics.。

fitting linear mixed-effects models using lme4

fitting linear mixed-effects models using lme4
The lme4 package for R provides functions to fit and analyze linear mixed models, generalized linear mixed models and nonlinear mixed models. In each of these names, the term “mixed” or, more fully, “mixed effects”, denotes a model that incorporates both fixed- and randomeffects terms in a linear predictor expression from which the conditional mean of the response can be evaluated. In this paper we describe the formulation and representation of linear mixed models. The techniques used for generalized linear and nonlinear mixed models will be described separately, in a future paper. The development of general software for fitting mixed models remains an active area of research with many open problems. Consequently, the lme4 package has evolved since it was first released, and continues to improve as we learn more about mixed models. However, we recognize the need to maintain stability and backward compatibility of lme4 so that it continues to be broadly useful. In order to maintain stability while continuing to advance mixed-model computation, we have developed several additional frameworks that draw on the basic ideas of lme4 but modify its structure or implementation in various ways. These descendants include the MixedModels package in Julia, the lme4pureR package in R, and the flexLambda development branch of lme4. The current article is largely restricted to describ-

lsmeans包用户说明说明书

lsmeans包用户说明说明书

Package‘lsmeans’October13,2022Type PackageTitle Least-Squares MeansVersion2.30-0Date2018-11-02Depends emmeans(>=1.3),methods,R(>=3.2)SuggestsByteCompile yesDescription Obtain least-squares means for linear,generalized linear,and mixed pute contrasts or linear functions ofleast-squares means,and comparisons of slopes.Plots and compact letter displays.Least-squares means were proposed inHarvey,W(1960)``Least-squares analysis of data with unequal subclass numbers'',Tech Report ARS-20-8,USDA National Agricultural Library,and discussedfurther in Searle,Speed,and Milliken(1980)``Population marginal meansin the linear model:An alternative to least squares means'',The American Statistician34(4),216-221<doi:10.1080/00031305.1980.10483031>.NOTE:lsmeans now relies primarily on code in the'emmeans'package.'lsmeans'will be archived in the near future.License GPL-2|GPL-3NeedsCompilation noAuthor Russell Lenth[aut,cre,cph]Maintainer Russell Lenth<***********************>Repository CRANDate/Publication2018-11-0223:50:11UTCR topics documented:lsmeans-package (2)auto.noise (2)ref.grid (3)ref.grid-class (4)transition (4)12auto.noise Index5 lsmeans-package Least-squares meansDescriptionThis package provides methods for obtaining so-called least-squares means for factor combinations in a variety offitted linear models.It can also compute contrasts or linear combinations of these least-squares means,(several standard contrast families are provided),and in addition can estimate and contrast slopes of trend lines.Popular adjustments for multiple-comparisons are provided,as well as graphical ways of displaying the results.Almost the entire codebase for lsmeans now resides in the emmeans package(named for the more general term,“estimated marginal means”).lsmeans exists only as a transitional entity for the few remaining packages that depend on it.Author(s)Russell V.Lenth(author),Maxime Hervé(contributor)Maintainer:Russ Lenth<***********************>ReferencesRussell V.Lenth(2016)Least-Squares Means:The R Package lsmeans.Journal of Statistical Software,69(1),1-33.doi:10.18637/jss.v069.i01Searle S.R.Speed liken G.A.(1980)Population marginal means in the linear model:An alternative to least squares means.The American Statistician34(4),216-221.auto.noise Data setsDescriptionThe datasets‘auto.noise’,‘feedlot’,‘fiber’,‘MOats’,‘nutrition’,and‘oranges’are provided in casea user customarily loads the data from lsmeans.But the same datasets are provided in the emmeanspackage,and they are documented there.Usageauto.noiseAuthor(s)Russell V.Lenthref.grid3 ref.grid Create a reference grid from afitted modelDescriptionThese functions are provided in lsmeans because they have been renamed in emmeansUsageref.grid(object,...)recover.data(object,...)lsm.basis(object,...)Argumentsobject A model object in a supported class....Additional arguments passed to companion functions in the emmeans package.Valuelsmeans now passes all its computations to emmeans,and the return values are thus what is returned by the corresponding functions ref_grid,recover_data,and emm_basis,respectively.Author(s)Russell V.LenthExamplesfiber.lm<-lm(strength~machine+diameter,data=fiber)rg<-ref.grid(fiber.lm,at=list(diameter=c(20,24,28)))rg#Note this is an emmGrid object defined in emmeans.The old"ref.grid"#class is now an extension of this:r.g.<-new("ref.grid",rg)lsmeans(r.g.,"machine")4transition ref.grid-class The ref.grid and lsmobj classesDescriptionThe codebase for lsmeans is now mostly in emmeans.These two classes are simple extensions of the emmGrid class defined in emmeans,and are provided as support for objects created in older versions of lsmeans.For details,see emmGrid-class.Author(s)Russell V.Lenthtransition Transition to emmeansDescriptionThe lsmeans package is being deprecated and further development will take place in its successor, ers may use emmeans in almost exactly the same way as lsmeans,but a few function names and internal details are changed.DetailsIn transitioning to emmeans,users willfind that the vignettes are constructed quite differently and that,in those and in the documentation,emphasis is placed on“estimated marginal means”rather than“least-squares means”.The term“estimated marginal means”is broader and more appropriate for use with some models,e.g.ordinal regression,that don’t really involve least-squares methods.That is the reason for the change.Accordingly,emmeans users are encouraged to use the functions emmeans(),emtrends(),emmip(), etc.in lieu of lsmeans(),etc.The latter functions are still available in emmeans;they run the corresponding emmxxxx function and relabel the results.The emmeans package provides some functions that help convert scripts and R Markdownfiles con-taining lsmeans code so they will work in emmeans.There is also a function to convert ref.grid and lsmobj objects to the emmGrid objects used in emmeans.More extensive information is given in vignette("transition-from-lsmeans",package="emmeans").Author(s)Russell V.LenthIndex∗datasetsauto.noise,2∗htestlsmeans-package,2∗modelslsmeans-package,2ref.grid,3∗packagelsmeans-package,2∗regressionlsmeans-package,2ref.grid,3auto.noise,2contrast(ref.grid),3emm_basis,3emmeans-transition(transition),4feedlot(auto.noise),2fiber(auto.noise),2lsm.basis(ref.grid),3lsmeans(ref.grid),3lsmeans-package,2lsmobj-class(ref.grid-class),4MOats(auto.noise),2nutrition(auto.noise),2oranges(auto.noise),2recover.data(ref.grid),3recover_data,3ref.grid,3ref.grid-class,4ref_grid,3summary.ref.grid(ref.grid),3transition,45。

asremlPlus与asreml混合模型分析:实例说明书

asremlPlus与asreml混合模型分析:实例说明书

Using asremlPlus,in conjunction with asreml,to do a linear mixed model analysis of a wheat experiment using hypothesis testsChris Brien05November,2023This vignette shows how to use asremlPlus(Brien,2023),in conjunction with asreml(Butler et al.,2020), to employ hypothesis tests to select the terms to be included in a mixed model for an experiment that involves spatial variation.It also illustrates diagnostic checking and prediction production and presentation for this experiment.Here,asremlPlus and asreml are packages for the R Statistical Computing environment(R Core Team,2023).It is divided into the following main sections:1.Set up the maximal model for this experiment2.Perform a series of hypothesis tests to select a linear mixed model for the data3.Diagnostic checking using residual plots and variofaces4.Prediction production and presentation1.Set up the maximal model for this experimentlibrary(knitr)opts_chunk$set("tidy"=FALSE,comment=NA)suppressMessages(library(asreml,quietly=TRUE))##Offline License checked out Sun Nov513:26:552023packageVersion("asreml")##[1]'4.2.0.276'suppressMessages(library(asremlPlus))packageVersion("asremlPlus")##[1]'4.4.20'suppressMessages(library(qqplotr,quietly=TRUE))options(width=100)Get data available in asremlPlusThe data are from a1976spring wheat experiment and are taken from Gilmour et al.(1995).An analysis is presented in the asreml manual by Butler et al.(2020,Section7.6),although they suggest that it is a barley experiment.data(Wheat.dat)Fit the maximal modelIn the following a model isfitted that has the terms that would be included for a balanced lattice.In addition,a term WithinColPairs has been included to allow for extraneous variation arising between pairs of adjacentlanes.Also,separable ar1residual autocorrelation has been included.This model represents the maximalanticipated model,current.asr<-asreml(yield~WithinColPairs+Variety,random=~Rep/(Row+Column)+units,residual=~ar1(Row):ar1(Column),maxit=30,data=Wheat.dat)Warning in asreml(yield~WithinColPairs+Variety,random=~Rep/(Row+:Some components changedby more than1%on the last iterationThe warning from asreml is probably due to a bound term.Initialize a testing sequence by loading the currentfit into an asrtests objectA label and the information criteria based on the full likelihood(Verbyla,2019)are included in thetest.summary stored in the asrtests object.current.asrt<-as.asrtests(current.asr,NULL,NULL,label="Maximal model",IClikelihood="full")Warning in infoCriteria.asreml(asreml.obj,IClikelihood=ic.lik,bound.exclusions=bound.exclusions): RepWarning in asreml(fixed=yield~WithinColPairs+Variety,random=~Rep/(Row+:Log-likelihoodnot convergedCheck for and remove any boundary termscurrent.asrt<-rmboundary(current.asrt,IClikelihood="full")Warning in infoCriteria.asreml(asreml.obj,IClikelihood=ic.lik):The following bound terms were discou RepWarning in asreml(fixed=yield~WithinColPairs+Variety,random=~Rep/(Row+:Log-likelihoodnot convergedsummary(current.asrt$asreml.obj)$varcompcomponent std.error z.ratio bound%chRep:Row 4.293282e+033.199458e+03 1.3418779P0.0Rep:Column 1.575689e+021.480357e+030.1064398P0.7units 5.742689e+031.652457e+03 3.4752438P0.0Row:Column!R 4.706787e+042.515832e+04 1.8708669P0.0Row:Column!Row!cor7.920301e-011.014691e-017.8056280U0.0Row:Column!Column!cor8.799559e-017.370402e-0211.9390486U0.0print(current.asrt,which="testsummary")####Sequence of model investigations(If a row has NA for p but not denDF,DF and denDF relate to fixed and variance parameter numbers)terms DF denDF p AIC BIC action1Maximal model266NA1646.1291742.469Starting model2Rep1NA NA1646.1291742.469BoundaryRep has been removed because it has been constrained to zero.Following the recommendation ofLittel et al.(2006,p.150),the bound on all variance components is set to unconstrained(U)using setvariances.asreml so as to avoid bias in the estimate of the residual variance.Alternatively,one couldmove Rep to thefixed model.Unbind Rep,Row and Column components and reload into an asrtests objectcurrent.asr<-setvarianceterms(current.asr$call,terms=c("Rep","Rep:Row","Rep:Column"),bounds="U")Warning in asreml(fixed=yield~WithinColPairs+Variety,random=~Rep/(Row+:Some components changed by more than1%on the last iterationcurrent.asrt<-as.asrtests(current.asr,wald.tab=NULL,test.summary=current.asrt$test.summary,IClikelihood="full",label="Max model&Unbound components") current.asrt<-rmboundary(current.asrt)summary(current.asrt$asreml.obj)$varcompcomponent std.error z.ratio bound%chRep-2458.34858411.197491e+03-2.0529167U0.0Rep:Row5008.71514863.401335e+03 1.4725732U0.0Rep:Column916.46411981.699576e+030.5392309U0.2units5959.02208171.609649e+03 3.7020634P0.0Row:Column!R46637.63034292.724392e+04 1.7118545P0.0Row:Column!Row!cor0.81505901.000281e-018.1483012U0.0Row:Column!Column!cor0.88568247.492514e-0211.8208968U0.0print(current.asrt,which="testsummary")####Sequence of model investigations(If a row has NA for p but not denDF,DF and denDF relate to fixed and variance parameter numbers)terms DF denDF p AIC BIC action1Maximal model266NA1646.1291742.469Starting model2Rep1NA NA1646.1291742.469Boundary3Max model&Unbound components267NA1647.1931746.544Starting modelprint(current.asrt,which="pseudoanova")####Pseudo-anova table for fixed termsWald tests for fixed effects.Response:yieldDf denDF F.inc Pr(Intercept)1 1.7153.5000.0115WithinColPairs115.6 2.5450.1307Variety2476.110.1100.0000Now the Rep component estimate is negative.The test.summary output has been extended,by supplying the previous test.summary to as.asrtests,toshow that there is a new starting model.The pseudo-anova table shows that Varieties are highly significant(p<0.001)2.Perform a series of hypothesis tests to select a linear mixed model for the dataThe hypothesis tests in this section are Wald tests forfixed terms,with denominator degrees of freedom calculated using the Kenward-Rogers adjustment(Kenward and Rogers(1997),and Restricted Maximum Likelihood Ratio Tests(REMLRT)for random terms.Check the term for within Column pairs(a post hoc factor)The information criteria based on the full likelihood(Verbyla,2019)is also included in the test.summarystored in the asrtests object.current.asrt<-testranfix(current.asrt,term="WithinColPairs",drop.fix.ns=TRUE,IClikelihood="full")Warning in asreml(fixed=yield~Variety,random=~Rep/(Row+Column)+:Some components changed by more than1%on the last iterationWarning in asreml(fixed=yield~Variety,random=~Rep/(Row+Column)+:Some components changed by more than1%on the last iterationprint(current.asrt)####Summary of the fitted variance parameterscomponent std.error z.ratio bound%chRep-2391.94899391.194581e+03-2.0023338U0.4Rep:Row5035.53110543.406006e+03 1.4784269U0.3Rep:Column761.95356221.612103e+030.4726458U1.2units5933.21337941.610805e+03 3.6833848P0.1Row:Column!R45970.83830272.635124e+04 1.7445415P0.0Row:Column!Row!cor0.81016159.995498e-028.1052641U0.1Row:Column!Column!cor0.88469707.503039e-0211.7911827U0.0####Pseudo-anova table for fixed termsWald tests for fixed effects.Response:yieldDf denDF F.inc Pr(Intercept)1 1.7158.900.0112Variety2476.810.270.0000####Sequence of model investigations(If a row has NA for p but not denDF,DF and denDF relate to fixed and variance parameter numbers)terms DF denDF p AIC BIC action1Maximal model26 6.0NA1646.1291742.469Starting model2Rep1NA NA1646.1291742.469Boundary3Max model&Unbound components267.0NA1647.1931746.544Starting model4WithinColPairs115.60.13071645.3251741.666DroppedIt is clear in the call to testranfix that the model is being changed by dropping the withinColPairsterm,which could also be achieved using update.asreml.However,an asremlPlus model-changing functionoperates on an asrtests object,that includes an asreml object,and,except for changeTerms.asrtests,results in an asrtests object that may contain the changed model or the supplied model depending onthe results of hypothesis tests or comparisons of information criteria.In addition,the result of the test or comparison will be added to a test.summary data.frame stored in the new asrtests object and,if the modelwas changed,the wald.tab in the new asrtests object will have been updated for the new model.In this case,as can be seen from the summary of current.asrt after the call,the p-value for the withinColPairs was greater than0.05and so now the model stored in current.asrt does not include withinColPAirs.The wald.tab has been updated for the new model.Test the nugget termThe nugget term represents non-spatial variance,such as measurement error.It isfitted using the asremlreserved word units.current.asrt<-testranfix(current.asrt,"units",positive=TRUE,IClikelihood="full")Warning in asreml(fixed=yield~Variety,random=~Rep+Rep:Row+Rep:Column,:Some components changed by more than1%on the last iterationWarning in asreml(fixed=yield~Variety,random=~Rep+Rep:Row+Rep:Column,:Some components changed by more than1%on the last iterationTest Row autocorrelationWe begin testing the autocorrelation by dropping the Row autocorrelation.Because of messages about the instability of thefit,iterate.asrtests is used to execute extra iterations of thefitting process.current.asrt<-testresidual(current.asrt,"~Row:ar1(Column)",label="Row autocorrelation",simpler=TRUE,IClikelihood="full")Warning in asreml(fixed=yield~Variety,random=~Rep/(Row+Column)+:Some components changed by more than1%on the last iterationWarning in asreml(fixed=yield~Variety,random=~Rep/(Row+Column)+:Some components changed by more than1%on the last iterationcurrent.asrt<-iterate(current.asrt)Test Column autocorrelation(depends on whether Row autocorrelation retained)The function getTestPvalue is used to get the p-value for the Row autocorrelation test.If it is significantthen the Column autocorrelation is tested by by dropping the Column autocorrelation,while retainingthe Row autocorrelation.Otherwise the model with just Row autocorrelation,whosefit is returned viacurrent.asrt after the test,is compared to one with no autocorrelation.(p<-getTestPvalue(current.asrt,label="Row autocorrelation"))[1]4.676754e-06{if(p<=0.05)current.asrt<-testresidual(current.asrt,"~ar1(Row):Column",label="Col autocorrelation",simpler=TRUE,IClikelihood="full")elsecurrent.asrt<-testresidual(current.asrt,"~Row:Column",label="Col autocorrelation",simpler=TRUE,IClikelihood="full")}Warning in DFdiff(bound.h1,bound.h0,DF=DF,bound.exclusions=bound.exclusions):There were a total The following bound terms occur in only one of the models compared and so were discounted:Row:Column!Row!corOutput the resultsprint(current.asrt,which="test")####Sequence of model investigations(If a row has NA for p but not denDF,DF and denDF relate to fixed and variance parameter numbers)terms DF denDF p AIC BIC action1Maximal model26 6.0NA1646.1291742.469Starting model2Rep1NA NA1646.1291742.469Boundary3Max model&Unbound components267.0NA1647.1931746.544Starting model4WithinColPairs115.60.13071645.3251741.666Dropped5units1NA0.00061645.3251741.666Retained6Row autocorrelation1NA0.00001645.3251741.666Unswapped7Col autocorrelation2NA0.00001645.3181741.658UnswappedprintFormulae(current.asrt$asreml.obj)####Formulae from asreml objectfixed:yield~Varietyrandom:~Rep/(Row+Column)+unitsresidual:~ar1(Row):ar1(Column)summary(current.asrt$asreml.obj)$varcompcomponent std.error z.ratio bound%chRep -2385.86975511.211207e+03-1.9698276U 0.0Rep:Row 5027.71232533.415391e+03 1.4720753U 0.0Rep:Column 753.59135361.609865e+030.4681086U 0.6units 5920.35470381.611274e+03 3.6743304P 0.0Row:Column!R 45870.09715952.623601e+04 1.7483638P 0.0Row:Column!Row!cor 0.80987861.001805e-018.0841906U 0.0Row:Column!Column!cor 0.88457687.510598e-0211.7777144U 0.0The test.summary shows is that the model with Row and without Column autocorrelation failed to converge.The asreml.obj in current.asrt contains the model selected by the selection process,which has been printed using printFormulae.asrtests .It is clear that no changes were made to the variance terms.3.Diagnosing checking using residual plots and variofacesGet current fitted asreml object and update to include standardized residuals current.asr <-current.asrt $asreml.objcurrent.asr <-update (current.asr,aom=TRUE )Wheat.dat $res <-residuals (current.asr,type ="stdCond")Wheat.dat $fit <-fitted (current.asr)Do diagnostic checkingDo residuals-versus-fitted values plotwith (Wheat.dat,plot (fit,res))1200140016001800−3−2−1012fit r e sPlot variofacesvariofaces (current.asr,V=NULL ,units="addtores",maxiter=50,update =FALSE ,ncores =parallel ::detectCores ())The variofaces are the lag1plots of the sample semivariogram with simulated confidence envelopes(Stefanova et al.,2009).Plot normal quantile plotThe plot is obtained using the ggplot function with extensions available from the qqplotr package(Almeida, A.,Loy,A.and Hofmann,H.,2023).4.Prediction production and presentationGet Variety predictions and all pairwise prediction differences and p-values Var.diffs<-predictPlus(classify="Variety",asreml.obj=current.asr,error.intervals="halfLeast",wald.tab=current.asrt$wald.tab,sortFactor="Variety",tables="predictions")####Predictions for yield from VarietyNotes:-The predictions are obtained by averaging across the hypertable calculated from model terms constructed solely from factors inthe averaging and classify sets.-Use'average'to move ignored factors into the averaging set.-The ignored set:Rep,Row,Column,units-Variety is included in this prediction-(Intercept)is included in this prediction-units is ignored in this predictionVariety predicted.value standard.error upper.halfLeastSignificant.limit 1101168.989120.47681228.315 211242.750119.81041302.076 391257.137119.97081316.463 4161285.718119.94001345.045 5141293.526119.92271352.853 *******.653120.29291372.979 7111322.159120.19641381.485 871374.447120.24071433.773 931394.070120.40321453.396 1041410.980120.10551470.306 11121444.557120.60341503.883 1281453.396120.59401512.723 13151458.383120.43461517.709 1451473.782120.44551533.108 15171487.828120.28961547.154 1661498.294120.11891557.620 17211517.121120.22621576.447 1821520.466119.63221579.792 19241533.769120.29951593.095 20181541.148120.36641600.474 21251575.795120.51421635.121 22221610.482120.32811669.808 23131610.762120.45751670.088 24201627.971120.23281687.297 25191652.992120.34351712.318 lower.halfLeastSignificant.limit est.status11109.663Estimable21183.424Estimable31197.811Estimable41226.392Estimable51234.200Estimable61254.327Estimable71262.832Estimable81315.120Estimable91334.743Estimable101351.653Estimable111385.231Estimable121394.070Estimable131399.057Estimable141414.456Estimable151428.501Estimable161438.968Estimable171457.795Estimable181461.140Estimable191474.443Estimable201481.821Estimable211516.468Estimable221551.156Estimable231551.436Estimable241568.645Estimable251593.666EstimableLSD valuesminimum LSD=114.0128mean LSD=118.6523maximum LSD=123.3578(sed range/mean sed=0.0788)We have set error.intervals to halfLeast so that the limits for so that the limits for each prediction±(0.5LSD)are calculated.When these are plotted overlapping error bars indicate predictions that are not significant,while those that do not overlap are significantly different(Snee,1981).Also set was sortFactor,so that the results would be ordered for the values of the predictions for Variety. The function predictPlus returns an alldiffs object,a list consisting of the following components:•predictions:the predictions,their standard errors and error intervals;•vcov:the variance matrix of the predictions;•differences:all pairwise differences between the predictions,•p.differences:p-values for all pairwise differences between the predictions;•sed:the standard errors of all pairwise differences between the predictions;•LSD:the mean,minimum and maximum LSDs.Plot the Variety predictions,with halfLSD intervals,and the p-values plotPredictions(Var.diffs$predictions,classify="Variety",y="predicted.value",error.intervals="half")ReferencesAlmeida,A.,Loy,A.and Hofmann,H.(2023)qqplotr:Quantile-Quantile plot extensions for‘ggplot2’, Version0.0.6.https:///package=qqplotr/or https:///aloy/qqplotr/.Brien,C.J.(2023)asremlPlus:Augments ASReml-R infitting mixed models and packages generally in exploring prediction differences.Version4.4.20.https:///package=asremlPlus/or /rpackages/.Butler,D.G.,Cullis,B.R.,Gilmour,A.R.,Gogel,B.J.and Thompson,R.(2023).ASReml-R Reference Manual Version4.2.VSN International Ltd,https:///.Gilmour,A.R.,Thompson,R.,&Cullis,B.R.(1995).Average Information REML:An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models.Biometrics,51,1440–1450.Kenward,M.G.,&Roger,J.H.(1997).Small sample inference forfixed effects from restricted maximum likelihood.Biometrics,53,983-997.Littell,R.C.,Milliken,G.A.,Stroup,W.W.,Wolfinger,R.D.,&Schabenberger,O.(2006).SAS for Mixed Models(2nd ed.).Cary,N.C.:SAS Press.R Core Team(2023)R:A language and environment for statistical computing.Vienna,Austria:R Foundation for Statistical Computing.https:///.Snee,R.D.(1981).Graphical Display and Assessment of Means.Biometrics,37,835–836.Stefanova,K.T.,Smith,A.B.&Cullis,B.R.(2009)Enhanced diagnostics for the spatial analysis offield trials.Journal of Agricultural,Biological,and Environmental Statistics,14,392–410.Verbyla,A.P.(2019).A note on model selection using information criteria for general linear models estimated using REML.Australian&New Zealand Journal of Statistics,61,39-50.https:///10.1111/anzs.12254/.。

一般混合线性模型SAS的MIXED过程实现_混合线性模型及其SAS软件实现_一_

一般混合线性模型SAS的MIXED过程实现_混合线性模型及其SAS软件实现_一_

一般混合线性模型SAS的M IXED过程实现———混合线性模型及其SAS软件实现(一)山西医科大学卫生统计教研室(030001) 张岩波 何大卫 刘桂芬 王琳娜 郭明英 【提 要】 目的 系统结构数据在医学领域广泛存在,其统计分析方法各异,可统称之为混合模型。

本文研讨其实现方法。

方法 以多水平模型例证一般混合线性模型的SAS M IX ED实现过程。

结果 以JSP数据为实例显示SAS的拟合结果与M Ln相一致。

结论 SAS M IXED可灵活地拟合包括多水平模型的各类混合模型。

【关键词】 系统结构数据 混合线性模型 多水平模型 M IX ED过程 近些年,国内医学统计学界对系统结构数据有了较多的认识,并进行了大量实效的研究和应用。

徐勇勇教授对系统结构数据做了全面的表述〔1〕。

由于常规的统计方法分析这类数据时忽略了误差结构,因此分析方法多采用以下模型:混合线性模型(Mixed lin-ear,M LM)、分层线性模型(Hierarchical linear, H LM)、广义线性混合模型(Generalized linear mixed, GLM M)、分层广义线性模型(Hierarchical generalized linear,HGLM)、多水平模型(Multilevel,M LM)、方差成分模型(Variance components,VCM)、随机系数模型(Random coefficients,RCM)等,以下且统称之为混合模型。

分析模型相应的软件有自行开发的软件(如陈长生博士针对重复测量数据自行开发的REP软件)及国外开发的专业软件,如M Ln(或M lw iN)软件,其他还有BUGS、H LM、VARCL等软件。

由于至今各种方法仍处于发展完善阶段,加之工具软件的限制,大大制约了此类方法的实际应用。

目前国内SAS软件已相当普及,其新增的M IXED模块及宏程序GLIM-M IX、NLINM IX可以有效、灵活地拟合各类混合模型,无疑为上述数据提供了有力的分析工具〔2,3〕。

linear mixed model random effect analysis

linear mixed model random effect analysis

linear mixed model random effect analysis
线性混合模型(Linear Mixed Model,LMM)是一种统计模型,用于分析具有嵌套结构或重复测量数据的数据集。

它结合了固定效应和随机效应,以便更好地解释和预测数据中的变异。

在LMM中,固定效应是指可以量化的、稳定的效应,通常用于解释数据中的变异。

随机效应则是指由于随机因素引起的效应,通常用于解释数据中的随机误差或变异。

在随机效应分析中,我们通常关注的是随机效应的估计和解释,而不是固定效应。

随机效应分析可以帮助我们理解数据中的重复测量或嵌套结构,并确定这些结构对总体平均值和变异的影响。

在LMM中,随机效应通常包括组内随机效应和组间随机效应。

组内随机效应是指在同一个组内的不同观测值之间的随机变异,而组间随机效应则是指不同组之间的随机变异。

这些随机效应可以影响模型的估计和预测,因此在进行LMM分析时,我们需要考虑这些效应的影响。

在进行LMM分析时,我们通常使用统计软件(如R、Stata或SAS)来拟合模型并进行推断。

这些软件提供了各种工具和函数,以便我们进行模型选择、参数估计、假设检验和模型诊断等操作。

通过拟合LMM,我们可以更好地理解数据中的变异和关系,并进行更准确的预测和推断。

广义线性混合效应模型试验设计若干

广义线性混合效应模型试验设计若干

广义线性混合效应模型试验设计若干IntroductionThe Generalized Linear Mixed Effects Model (GLMM) is a powerful statistical tool used for analyzing data in various fields including ecology, social sciences, economics, and psychology. GLMM combines both fixed effects and random effects to account for within-group and between-group variability. It is a flexible model that can handle non-normal data types and non-linear relationships. In this article, we will discuss the experimental design considerations for GLMMs.Experimental Design Considerations for GLMMs1. Choosing the Appropriate Response VariableThe first step in designing a GLMM is to choose an appropriate response variable. The response variable should be continuous and non-negative. When using count data, the response variable should follow a Poisson distribution. If the data is binary or categorical, then logistic regression or multinomial regression should be used.2. Determining Fixed and Random EffectsFixed effects are variables that are assumed to have a systematic effect on the response variable. Examples of fixed effects include treatment group, time, and sex. Random effects are variables that contribute to within-group variation. Examples of random effects include subject ID, site, and block. It is important to carefully consider which variables to include as fixed or random effects.3. Sample Size ConsiderationsSample size considerations are important in any experimental design. In a GLMM, it is important to have enough observations and groups to capture the random effects. A rule of thumb is to have at least 5-10 groups per random effect included in the model. Additionally, a larger sample size will improve the power of the model.4. Choosing the Appropriate Link Function and DistributionThe link function connects the response variable to the linear predictor in the GLMM. Common link functions include the identity, logit, and log-link. The choice of the link function depends on the distribution of the response variable. The distribution of the response variable should also be carefully considered. Common distributions used in GLMMs include the Gaussian, Poisson, and binomial distributions.5. Controlling for CovariatesCovariates are variables that may influence the response variable but are not of interest in the study. It is important to control for covariates in the GLMM to prevent confounding. Examples of covariates include age, weight, and baseline measurements.6. Model Selection and ValidationSelecting an appropriate model is key to obtaining meaningful results. It is important to consider both statistical and biological significance of the model. Model selection can be done through stepwise selection or model comparison techniques such as AIC or BIC. The final model should be validated through checking model assumptions, assessing the fit of the model, and examining any influential observations.ConclusionGLMM is a powerful statistical tool used in various fields to analyze data. In this article, we have discussed the experimental design considerations for GLMMs. Careful consideration of the response variable, fixed and random effects, sample size, link function and distribution, controlling for covariates, and model selection and validation are key to obtaining valid results.。

五、混合能源模型(MixedEnergyModel)

五、混合能源模型(MixedEnergyModel)

五、混合能源模型(Mixed Energy Model)描述:对能源系统(从能源的开采、转化、运输、市场到最终能源需求)的模拟,通过系统仿真来预测各部门能源的供应能力、能源价格、需求量以及宏观经济参数,从而为国家制定能源战略和决策提供信息支持,因此既包括自顶向下的宏观经济模型,又包括自底向上的能源供应、需求模型。

这类模型是对整个能源-经济-环境系统的模拟和仿真,是一个复杂巨系统。

研究范围多是全球的、区域的或国家的,结构上也多是包括经济、供应、转化、需求、环境等模块的综合集成模型。

目前中国的应用及研究相对较少,由于模型涉及的技术和领域非常广泛,所以必须有足够的、专业的研究人员和时间做保证,才能完成这项复杂的系统工程。

代表性的混合能源模型是美国能源部(DOE)开发的NEMS模型和奥地利国际应用系统分析研究所(the International Institute for Applied Systems Analysis,IIASA)与世界能源委员会(the World Energy Council,WEC)合作开发的IIASA-WEC E3模型。

表4 典型综合模型的特点1.NEMS最具代表性的混合能源模型是NEMS (National Energy Modeling Systems)模型,由美国EIA/DOE于1993年开发的能源经济区域模型。

NEMS综合考虑了宏观经济、财政因素、世界能源市场、资源可获得性和成本、行为和技术选择标准、能源技术的成本和运行特性以及人口统计资料,反应了能源的生产、进口、转化、消费以及价格的情况。

但从文献上看,用“NEMS模型”当做关键词在期刊网上检索,没有一篇文献。

(1)NEMS的功能EIA把NEMS用来模拟在美国能源政策和能源市场上不同假设下的能源、经济、环境以及安全之间的影响。

NEMS通过制定能源产品的生产、转换、消费的经济决策,清晰地描述了美国国内能源市场,同时NEMS 还描述了能源技术。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

SELECTING THE BEST LINEAR MIXED MODEL USING PREDICTIVE APPROACHESByJun WangA project submitted to the faculty ofBrigham Young Universityin partial fulfillment of the requirements for the degree ofMaster of ScienceDepartment of StatisticsBrigham Young UniversityApril 2007BRIGHAM YOUNG UNIVERSITYGRADUATE COMMITTEE APPROVALof a project submitted byJun WangThis project has been read by each member of the following graduate committee and by majority vote has been found to be satisfactory._______________________________ ____________________________________ Date G. Bruce Schaalje, Chair_______________________________ ____________________________________ Date Gilbert W. Fellingham_______________________________ ____________________________________ Date H. Dennis TolleyBRIGHAM YOUNG UNIVERSITYAs chair of the candidate’s graduate committee, I have read the project of Jun Wang in itsfinal form and have found that (1) its format, citations, and bibliographical style areconsistent and acceptable and fulfill university and department style requirements; (2) itsillustrative materials including figures, tables, and charts are in place; and (3) the finalmanuscript is satisfactory to the graduate committee and is ready for submission to theuniversity library._______________________________ ____________________________________Date G. Bruce SchaaljeChair, Graduate CommitteeAccepted for the Department____________________________________ Scott D. GrimshawGraduate CoordinatorAccepted for the College____________________________________ Thomas W. SederbergAssociate Dean, College of Physical andMathematical SciencesABSTRACTSELECTING THE BEST LINEAR MIXED MODELUSING PREDICTIVE APPROACHESJun WangDepartment of StatisticsMaster of ScienceThe linear mixed model is widely implemented in the analysis of longitudinal data. Inference techniques and information criteria are available and well-studied for goodness-of-fit within the linear mixed model setting. Predictive approaches such as 2adj R , PRESS ,and adj CCC are available for the linear mixed model but require more research (Edward,2005). This project used simulation to investigate the performance of 2adj R , PRESS,adj CCC , Pseudo F-test and information criterion for goodness-of-fit within the linear mixed model framework. Marginal and conditional approaches for these predictive statistics were studied under different variance-covariance structures.For compound symmetry structure, the success rates for all 17 statistics (marginal and conditional PRESS, CCC, and 2R , F test, AIC and BIC) were high. The study suggested using marginal rather than conditional residuals for PRESS, CCC and 2R . Itsuggested using REML likelihood function which has the term ||log 211'∑=mi i i X X for AICand BIC. For CCC, 2R , and the information criterion, there was no difference for the various parameter number adjustments.For autoregressive order 1 plus random effect, the study suggested using conditional residuals for PRESS, marginal residuals for CCC and 2R , and using REMLfunction with the term ||log 211'∑=mi i i X X for AIC and BIC. Also there was no differencefor the different parameter number adjustments.The F-test performed well for all covariance structures. The study also indicated that characteristics of the data, such as the covariance structure, parameter values, and sample size, can greatly impact performance of various statistics. No one criterion is consistently better than the others in terms of selection performance in the simulation study.AcknowledgementsI would like to thank my advisor, Dr. Schaalje for his time and effort throughout this project. I am also grateful for the members of my graduate committee, and the Department of Statistics for all the help and support. Additionally, I would like to thank my family for all the love, support and encouragement they have given me these years.ContentsChapter1I n t r o d u c t i o n (1)2 Literature Review (3)2.1. Linear Mixed Models for Repeated measures (3)2.2. Model selection in linear mixed models (4)2.2.1 AIC and BIC (5)2.2.2 Pseudo F-tests and Likelihood Ratio Tests (7)2.2.3 Predictive Approaches (8)3 Methodology (12)3.1. Data Generation (12)3.2. Modeling the Data (14)4 Results (18)4.1. Statistics (18)4.2. Sample Size (20)4.3. Covariance Structure (20)5 Conclusions (34)6 Bibliography (36)Appendix A: Sample SAS Code (38)xiiList of Tables3.1: Parameter Values for Covariance Structures Used In the Simulations (13)3.2: Values of β Used In the Simulation (15)3. 3: Example of SAS PROC MIXED code (15)4.1: Success Rate of Choosing the Correct Model Using Different Statistics forC o v a r i a n c e S t r u c t u r e=C S(0.5) (21)4.2: Success Rate of Choosing the Correct Model Using Different Statistics forC o v a r i a n c e S t r u c t u r e=C S(0.9) (23)4.3: Success Rate of Choosing the Correct Model Using Different Statistics forC o v a r i a n c e S t r u c t u r e=A R R E(0.5) (25)4.4: Success Rate of Choosing the Correct Model Using Different Statistics forC o v a r i a n c e S t r u c t u r e=A R R E(0.9) (27)4.5: Average Success Rate of Choosing the Correct Model Using Different Statistics (33)xiiiList of FiguresFigure 1: Correct model selection rates for 17 criteria for the compound symmetry structure involving within subject correlation of 0.5 (29)Figure 2: Correct model selection rates for 17 criteria for the compound symmetry structure involving within subject correlation of 0.9 (30)Figure 3: Correct model selection rates for 17 criteria for the autoregressive order 1 plus random effect structure with autocorrelation of 0.5 (31)Figure 4: Correct model selection rates for 17 criteria for the autoregressive order 1 plus random effect structure with autocorrelation of 0.9 (32)xivxv1Chapter 1IntroductionThe linear mixed model is a powerful tool for the analysis of longitudinal data.Longitudinal data involves repeated measurements on the same individual or experimental unit. For the purposes of estimation and inference, the linear mixed model extends the ordinary linear model with independently and identically distributed Gaussian errors to a wide variety of correlated Gaussian data.Information criteria such as AIC (Akaike, 1973) and BIC (Schwarz, 1978) are available as goodness-of-fit statistics for model selection in the linear mixed model setting (Ferron et al., 2002; Gomez et al., 2005). Pseudo F-tests for linear mixed model selection are also available using various small sample approximations. The performance of these tests for complex mixed linear models was studied by Schaalje et al. (2002). Predictive criteria such as the 2R , PRESS (Allen, 1974), and CCC (Lin, 1989)statistics have been well established in model selection for standard linear models buthave received little study for linear mixed models (Edwards et al., 2005). However, many of these statistics can be adapted to the mixed model setting (Vonesh et al., 1996). Forexample, 2R =∑∑==−−−ni i ii ni iy yyy 1221)()ˆ(1 is defined in linear regression as the proportion of total variation in y that is explained by the model. 2R can be modified in the longitudinal setting by using either marginal or conditional predicted values (Schabenberger, 2004). CCC, an alternative measure of agreement, was modified by Vonesh et al. (1992) for usein mixed-effects regression settings. The PRESS statistic, a cross-validative measure of goodness-of-fit, was modified for the mixed linear model setting by Schabenberger (2004) and Christensen et al. (1992).This study uses simulations to investigate and compare the performance of AIC, BIC, pseudo F, 2R, PRESS and CCC for model selection within the linear mixed model framework. Marginal and conditional versions for these statistics (Vonesh et al., 1996) are studied under different variance-covariance structures.23Chapter 2Literature Review2.1. Linear Mixed Models for Repeated MeasuresIn the context of repeated measures data, let i y be an )1(×i n vector of observed data on the th i subject where },...,1{m i ∈, and m is the number of independent units (subjects). The mixed model can be written as:i i i i i u Z X y εβ++= , (1)where i X is an )(p n i × known constant design matrix for the th i subject, β is a )1(×p vector of unknown coefficients of the fixed effects, i Z is a )(q n i × known constant design matrix for the th i subject with rank q , i u is a )1(×q vector of unknownindividual-specific random effects, and i ε is a )1(×i n vector of random errors (Gurka, 2006). The vector of coefficients of the random effects, i u , is assumed to follow the normal distribution with mean vector 0 and covariance matrix G . i ε is assumed to follow a normal distribution with mean vector 0 and covariance matrix i R . The vectors i u and i ε are assumed independent; consequently,0),cov(=i i u εand i i i i i V R GZ Z y =+=')var(. Let ∑==mi i n N 1, and denote the unique parameters of i R and G by the )1(×k vector θ.The log-likelihood function for the linear mixed model is then∑∑==−−−−−−=m i m i i i i i i ML X y V X y V Nl 111')()()(21))(log(21)2log(2),(βθβθπθβ . (2)4The maximum likelihood estimator (MLE) is produced by numerically maximizing),(θβML l with respect to the unknown parameters. Since the MLE of θ is biased (Gurka,2006), the Restricted Maximum Likelihood estimator (REML) of θ is typically used.Assuming that rank ('1X L '2X 'm X ) is p , the REML estimator of θ is calculated by maximizing the likelihood function of a transformation of the original N observations to a new set of N – p observations (Gurka, 2006). The transformation is constructed such that the restricted log-likelihood function can be written as)ˆ()()'ˆ(21|)(|log 21|)(|log 21|'|log 21)2log(2)(1111'11βθβθθπθi i m i i i i i i mi i mi i mi i i REML X y V X y X V X V X X p N l −−−−−+−=∑∑∑∑=−−=== , (3)where βˆ is ∑∑=−−=−=mi i i i mi i i i y V X X V X 11111)ˆ(ˆ'))ˆ(ˆ'(ˆθθβ , (4)and i V ˆis the estimated covariance matrix iV . In practice, iterative procedures are used to maximize (2) or (3).It is interesting to note that SAS PROC MIXED (2003) computes the REMLestimates with a version of equation (3) that leaves out the term |'|log 211∑=mi i i X X (Gurka,2006). The same estimates are obtained, but the information criteria (to be discussed in section 2.2) are different. The log-likelihood function used by SAS will be denoted2REML l .52.2 Model selection in linear mixed models2.2.1 AIC and BICInformation criteria such as AIC (Akaike, 1973), BIC (Schwarz, 1978), AICC (Hurvich Tsai, 1989), and CAIC (Bozdogan, 1987) are available for selecting the covariance structure (i Z , G , and i R ) and the mean model (i X and β) in linear mixed models. In general, these information criteria are functions of both the maximized log-likelihood for a given model (l ) and a penalty term based on the number of parameters (s ) in the model (Gurka, 2006). The usual formulas for AIC and BIC are:s l AIC −= , (5)2/)(log *N s l BIC −= , (6) where *N is a function of N , p or m (usually *N =N or m under ML, and *N =p N − or m under REML), and k p s +=. When using ML estimation, generally k p s +=, the total number of parameters in the model. However, under REML estimation, k s =, the number of covariance parameters only (Gurka, 2006). A larger value of the information criterion for a given model indicates a preference for that model.Keselman et al. (1998) compared the effectiveness of AIC and BIC for detecting various population covariance structures. In particular, the criteria were compared in unbalanced nonspherical repeated measures designs having equal/unequal group sizes, various covariance matrices, and normally/nonnormally distributed responses. Their results show that neither approach uniformly selected the correct covariance structure. Indeed, for most of the structures investigated in the study, AIC and particularly BIC more frequently picked the wrong covariance structure than the correct covariance structure. Keselman’s results show that AIC chose the correct covariance structure 47%of the time while BIC was right 35% of the time. The negative results obtained for BICcould be due to the severe penalty imposed for the number of parameters (Keselman etal., 1998).Ferron et al. (2002) also studied the performance of AIC and BIC for choosing acovariance structure. Data were generated following an AR(1) structure with differentsample sizes, numbers of repeated measures, and levels of autocorrelation. AIC and BICwere computed for the true covariance structure, AR(1), and for the independent σ. The results show that AIC was more successful than BIC for every structure, I2combination of sample size, number of repeated measures, and level of autocorrelation.AIC chose the right covariance structure 79% of the time versus 66% for BIC. AIC andBIC both performed better when the sample size, the number of the repeated measures,and the level of autocorrelation were higher. The number of the repeated measures wasthe most influential factor and the effect of sample size was greater when the number ofrepeated measures was smaller.Gomez et al. (2005) studied the performance of AIC and BIC criteria in selectingthe true covariance structure. They generated data using 15 covariance structures. Theresults show that the success rate of AIC and BIC in choosing the correct covariancestructure depends greatly on the sample size and complexity of the covariance structures.Success rates were generally low, but they were higher for larger sample sizes andsimpler covariance structures. Pairing (a positive pairing refers to the case in whichthe largest sample size was associated with the covariance matrix containing the largestvalues) had only a small effect on the performance of AIC and BIC.6Gurka (2006) compared the performance of AIC and BIC for ML l , REML l , and2REML l (see section 2.1). He assessed linear mixed model selection under three scenarios:(1) selection of the fixed effects when the covariance structure is known, (2) selection ofthe random effects when the fixed effects are known, and (3) simultaneous selection ofthe fixed and random effects. The results show that all versions of the two criteriaselected the proper random effects structure over 90% of the time, no matter the truevariance or correlation. BIC performed slightly better than AIC. The results also showedthat BIC performed better overall when using REML l , while for AIC 2REML l selection wassuperior to REML l as well as ML l .Keselman et al. (1998) investigated the performance of AIC and BIC in thecontext of covariance structures common in social science. Their results provide aninteresting comparison to those reported by Gurka (Keselman et al. 2006). Thedifferences between Gurka’s results and Keselman’s results show that the performance ofAIC and BIC depends on the conditions investigated.2.2.2 Pseudo F-tests and Likelihood Ratio TestsIf C is a contrast matrix of rank q , a commonly used test statistic for Ho:0=βC isthe Wald F statistic q W F ddf q /,=, where)ˆ()')ˆ'(()'ˆ(111ββC C X V X C C W −−−= (7) and ddf is the denominator degrees of freedom. Under the null hypothesis, the Wald Fapproximately follows the F distribution with q and ddf degrees of freedom (Schaalje etal., 2002). There are two commonly used methods to compute ddf for the Pseudo F -test.The Fai and Cornelius (FC) method computes the ddf using spectral decomposition of theapproximate covariance matrix for the contrasts, together with repeated application of a method for single degree of freedom tests (Fai and Cornelius, 1996; Giesbrecht and Burns, 1985; Schaalje et al., 2002). Kenward and Roger (1997) suggested an improved method (KR) to calculate the Pseudo F-test. The KR method adjusts the estimated covariance matrix of the parameter estimates, computes a scale adjustment to the test statistic and computes the degrees of freedom of the Pseudo F test.The performance of these tests for selecting the model for the means (“mean structure”) in complex mixed linear models was studied by Schaalje et al. (2002). The result shows that Type I error rates of both methods are affected by covariance structure complexity, sample size, and imbalance. The KR method performs well in situations with complicated covariance structures when sample sizes are moderate to small and the design is balanced. The KR method is preferable to the FC method, although it also inflates Type I error rates for complex covariance structures combined with small sample sizes.Gomez et al. (2005) investigated Type I error rates after using AIC and BIC to select the covariance structures. The results show that Type I error rates of KR method hypothesis tests for best AIC and best BIC models were always higher than the target values. Type of selection criterion, number of subjects per treatment, number of repeated measures, effect tested, covariance structure, and pairing all affected the Type I error rates.2.2.3 Predictive ApproachesThere are two main distinctions between mixed models and linear models. First, in amixed model, the data can be considered in a conditional and a marginal sense. Theconditional distribution of u y | uses a particular realization of the random effects. Themarginal distribution of y uses quantities averaged over all possible values of the randomeffects. Correspondingly, there are conditional and marginal predictive approaches formixed models. Second, the estimates of the fixed effects β depend on the estimates ofthe covariance parameters (Schabenberger, 2004). Predictive approaches for modelselection in mixed models reflect these distinctions. R², the squared correlation coefficient, is defined in the simple linear model as, (8)where i i i yy r ˆ−=, i y is the ith observation, i y ˆ is the ith predicted value and y is the grand mean. The squared correlation coefficient can also be calculated in the mixedmodel setting, but differs depending on whether marginal residuals m r or conditionalresiduals c r are used (Vonesh, 1996). A marginal residual, βˆ'i i mi x y r −=, is the difference between the observed data and the estimated (marginal) mean. A conditionalresidual, u z x y r i i i ci ˆˆ'−−=β, is the difference between the observed data and thepredicted value of the observation. Mixed model versions of 2R can be adjusted for thenumber of parameters to obtain:, (9) ∑∑==−−=n i i n i iy y r R 12122)(1)1(12*2R k R adj −−=where p N N k −=* or )(*k p N N k +−=. The PRESS residual (Allen 1974) is the difference between the observed valueand the predicted value, where the predicted value is obtained without the observation inquestion. The PRESS statistic for the simple linear model can be computed without re-fitting the model by using the update formula (Schabenberger, 2004). In the longitudinalsetting such update formulas are not available because the observations are correlated(Schabenberger, 2004). The PRESS statistic can also be calculated for mixed models, butdiffers depending on whether the marginal or conditional residual is used. The concordance correlation coefficient (CCC) was described by Lin (1989) andmodified by Vonesh (1992) for use in mixed-effects regression settings. Given bivariatemeasurements 1Y and 2Y with means 1µ, 2µ, variances 21σ, 22σ, and correlation ρ, Lin(1989) defined the concordance correlation between 1Y and 2Y to be 2212221212212221221)(2)(])[(1µµσσσρσµµσσρ−++=−++−−=Y Y E c . (10) CCC is a measure of the degree of accuracy and precision to which pairs ofobservations correspond to the line of identity. The point estimate of c ρ based on arandom sample of n observations is ∑∑∑===−+−+−−−=n i n i i i ni i i c Y Y n Y Y Y YY Y112212222111221)()()()(1ˆρ. (11)Vonesh (1992) suggested adapting CCC to measure model fit. He defined the “model CCC” as ∑∑∑===−+−−+−−−−−=n i n i i i i i i i ni i i c y y N j y y j y y j y y j y yy y yy r 112''1')ˆ()ˆˆ()ˆˆ()()()ˆ()ˆ(1. (12) There are several advantages of using c r as a measure of model fit. First, c r is directly interpretable as a concordance correlation coefficient between observed and predicted values. Second, it does not require specification of a null model since the line of identity serves as a point of reference. In the mixed model setting, we can compute a marginal or conditional model CCC by using marginal or conditional estimates (Vonesh, 1996). Also, CCC can be adjusted for the number of parameters as )1(1*CCC k CCC adj −−=, (13) where p N N k −=* or )(*k p N N k +−=.Chapter 3MethodologyThis chapter describes the methods used for the simulation study. First, data were generated with a specified model for the means and with repeated measurements on the same individual, correlated according to some known covariance structure. Then, two mixed models (a full model and a reduced model) were fit to the data, and several predictive approaches were used to select one of the models. The proportion of runs in which the correct model was selected was recorded.An explanation of the data generation process is outlined in the first section of this chapter. The second section discusses the use of simulations to investigate the performance of 2adjR , PRESS,adj CCC , information criteria, and the F test for model selection. Marginal and conditional approaches for these statistics were studied under different variance-covariance structures. Data were simulated using PROC IML of SAS and analyzed with PROC MIXED.3.1. Data GenerationThe simulated data followed a repeated measures design. Each data set involved m =10 or 20 subjects, and a constant number of observations per subject 3(==n n i or 5 for 1=i to m ). The repeated measures followed one of four covariance structures, and the fixed effects followed one of six possible models. Thus, 2×2×4×6=96 simulation studies were carried out.The covariance structures were the same for each subject (V V i = for 1=i to m ). The four structures are denoted CS (0.5), CS (0.9), ARRE (0.5), and ARRE (0.9). Details on these covariance structures when the number of repeated measures was five are in Table 1. The upper left 3×3 matrix was used when the number of repeated measures was three.Table 3.1: Parameter Values for Covariance Structures Used In the SimulationsCompound symmetry (CS (0.5))151551555155551..........⎛⎞⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠Compound symmetry (CS (0.9))⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛19.9.9.9.19.9.9.19.9.19.1 Autoregressive order 1 plus randomeffect (ARRE(ρ=0.5))⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛16.4.3.25.16.4.3.16.4.16.1Autoregressive order 1 plus random effect (ARRE(ρ=0.9)) ⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛192.848.7832.72488.192.848.7832.192.848.192.1Data for the four described covariance structures were generated following the two step method of Ripley (1987). In the first step, a random multivariate normal vector was generated with E(y ) = 0 and Var(y ) =I , where I is the identity matrix. In the secondstep, this vector was multiplied by the Cholesky decomposition of the covariance matrix corresponding to the covariance structure in question. The resulting random vector had a mean vector of zero and covariance matrix of V (Table 1).3.2. Modeling the DataThe fixed effects part of the models were defined by a combination of designmatrices ([]210X X X X i =) and parameter vectors (⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=210ββββ). In the studies with n =5, i X was the same for all subjects, and was given by ⎟⎟⎟⎟⎟⎟⎠⎞⎜⎜⎜⎜⎜⎜⎝⎛−−=231121*********i X . In the studies with n=3, ⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−=031101031i X . These design matrices were combined with valuesof the parameter vector given in Table 2. For 6.02=β, some of the expected values change by as much as two standard deviations.Two mixed models were fitted to each data set using REML as implemented in PROC MIXED. For both models, the correct covariance structure was specified using RANDOM and REPEATED statements. Also, two models for the fixed effects were specified using the MODEL statement, the full model, and the reduced model. For the full model, the MODEL statement was specified as 1X y = 2X , whereas for the reduced model the MODEL statement was 1X y =.Table 3.2: Values of β Used In the Simulation⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=011β⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=05.11β ⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=1.11β ⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=2.11β⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=4.11β ⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=6.11βTable 3.3: Example of SAS PROC MIXED code *SAS code for compound symmetry structure;proc mixed data =datanew;by iter;class subj;model y=x1 x2/ S influence(iter=5 est) outpred=pred_full(keep=y iter x1 x2 subj Pred Resid)outpredm=predm_full(keep=y iter x1 x2 subj Pred Resid)DDFM =kenwardroger;random subj/s;ods output fitstatistics=cs;ods output tests3=pval_cs;ods output SolutionF=coef_fixed (keep=estimate iter) ;ods output SolutionR=coef_random;ods output influence=press (keep=iter pressres); run ; *SAS code for autoregressive order 1 plus random effect structure; proc mixed data =datanew; by iter;class subj;model y=x1 x2/ S influence(iter=5 est) outpred=pred_full(keep=y iter x1 x2 subj Pred Resid)outpredm=predm_full(keep=y iter x1 x2 subj Pred Resid)DDFM =kenwardroger;random subj/s;repeated /type =ar(1) subject =subj;ods output fitstatistics=arre;ods output tests3=pval_arre;ods output SolutionF=coef_fixed (keep=estimate iter) ;ods output SolutionR=coef_random;ods output influence=press (keep=iter pressres); run ;Several versions of PRESS,adj CCC , 2adjR , the pseudo F-test and information criteria (AIC and BIC) were calculated and recorded for each model. Some of the information criteria and the F-test were given by default by PROC MIXED. The Kenward-Roger method for the F-test was specified by using the ddfm=kenwardrogeroption in the MODEL statement. For PRESS, adj CCC and 2adj R , both marginal andconditional residuals were considered. The marginal residual βˆ'i i m x y r −= is the difference between the observed data and the estimated (marginal) mean. The conditionalresidual u z x y r i i i c ˆˆ'−−=β is the difference between the observed data and the predictedvalue of the observation. Marginal and conditional residuals were given by using outpredm and outpred options in the model statement of PROC MIXED. The ODS (output delivery system) was used to capture the fitted statistics and the residuals.PRESS, adj CCC and 2adj R were then calculated using PROC IML.The PRESS statistic is the sum of the squared differences between the observed values and the predicted values, where each predicted value is obtained without theobservations in question. Two versions of PRESS were calculated, the marginal andconditional PRESS. A PROC MIXED option was used to compute the marginal PRESS residuals. A SAS MACRO (Appendix A) was used to calculate the conditional PRESS statistics.adj CCC was calculated as∑∑∑===−+−−+−−−−−=n i n i i i i i i i n i i iy y N j y y j y y j y y j y y y y yy CCC 112''1')ˆ()ˆˆ()ˆˆ()()()ˆ()ˆ(1 and。

相关文档
最新文档