Generalized Linear Mixed Models

合集下载

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

从结果看到，考虑了不同中心的差异，Ａ、Ｂ两种药物的副作用的发生有差异（ β１＝－０．９２９８，Ｐ＝０．０３５４，Ｂ药更易发生副作用），不同实验中心的间的方差为１．５８０９，Ｐ＝０．１８４２，差异无统计学意义。
设随机效应ｕｉ的密度函数为：ｆｕ（ｕｉ，Ｇ），有边际似然函数：
& Ｌｉ ! β， γ" ＝Ｌｉ ! β，ｕｉ " ｆｕ ! ｕｉ， γ" ｄｕｉｕｉ
&#% $ ｎｉ
＝
ｕｉ
ｆｙ ! ｙｉｊｕｉ，
ｊ＝１
Ｘｉｊ，
β"ｆｕ ! ｕｉ，
γ" ｄｕｉ
（１．５）
γ是Ｕｉ的方差协方差矩阵，是Ｇ的参数估计值。
得到似然函数：
% Ｌ! β， γ" ＝Ｌｉ ! β， γ" （１．６）ｉ
从上式可以看到，计算似然函数比线性混合效应模型复杂
· ２１０４ ·
现代预防医学２００７年第３４卷第１１期ＭｏｄｅｒｎＰｒｅｖｅｎｔｉｖｅＭｅｄｉｃｉｎｅ，２００７，Ｖｏｌ．３４，ＮＯ．１１
得多，需要解决随机效应ｕｉ的高维积分的问题，许多最大化似然函数的近似的推断方法被提出，目前积分近似方法主要有Ｌａｐｌａｃｅ近似（ＬｉｕａｎｄＰｉｅｒｃｅ，１９９３），ＡｄａｐｔｉｖｅＧａｕｓｓｉａｎ积分，一阶泰勒（ｆｉｒｓｔ－ｏｒｄｅｒＴａｙｌｏｒ）序列展开近似。２实例分析
为了研究Ａ、Ｂ两种药物的的副作用情况，研究者随机选取了１５个中心做临床实验，在每个中心中，随机抽取一定数量的病人，其中ｎＡ个病人接受Ａ药物，ｎＢ个病人接受Ｂ药物。数据格式见表１。

重复测量诊断试验的ROC曲线广义线性混合效应模型

－1 珔 “有病” 时患者诊断试验结果 T1 大于或等于 F 0， x （ p）的概率。也就是说ＲOC 曲线是ＲOC x （ p ）与 p 构建的
积标准正态分布函数的逆函数。ＲOC 曲线下面积（ AUC ）是评价诊断试验最常用（阳性）和“不的一个指标，它表示诊断系统中“患病 ” （阴性）诊断结果分布与“金标准 ” 的重叠程度，患病” 体现了诊断试验的价值，面积越大诊断价值越高。理论上 AUC 的取值范围为 0. 5 1 ，两端点分别表示完全无价值的诊断及完善的诊断。一般ＲOC 曲线面积在 0. 50 0. 70 ，表示诊断准确度低，在 0. 70 0. 90 ，表 0. 90 以上，示诊断准确度中等，认为诊断准确度较高。 2. 参数估计方法本文采用马尔科夫蒙特卡洛（ M arkov chain M onte M CM C ）的贝叶斯方法来估计ＲOC 曲线广义线 Carlo ，性混合效应模型。相比于极大似然估计方法，该方法且用于估计的随机效应变量个数可以更加灵活准确，是任意的。本文用 WinBUGS 软件来进行计算。具体［3 － 5 ］：的参数估计过程如下第一步，给定假阳性率集， Γ = （ p）基于 Alonzo 和 Pepe （ 2001 ）的模拟研究，可选择 50 个等间距的 FPＲs， …， 50 /51 ），即 Γ = （ 1 /51 ，不仅可以保证模型参数估计的有效性和稳健性，还能节约模型估计的运算时间，实现用较小的假阳性率集获得较大假阳性率集的统计功效。 D 的先验分布第二步，确定参数 γ， β，（ γ | A） U（－ A， A）（ β | B） U （－ B， B）－1 （ D | V， v ） Wishart（ V － 1 ， v） M CM C 过程第三步， D 的后验分布， β，此处采用 Gibbs 抽样为获得 γ， k = 1， …， n，方法，迭代更新方法如下，对于第 k 步，直至收敛： p， y）（ 1 ）（ γ | λ （ k），（ k）给定 λ ，给定假阳性率集 Γ = （ p ），对于大样本，（ k）（ k）（ γ|λ ， p， y ）近似服从均值为极大似然估计值 ^ γ ，

Bayesian covariance selection in generalized linear mixed models

determine which factors vary in their eﬀects across women. For example, do the eﬀects of aging, recent oral contraceptive use, and smoking vary? Addressing such hypotheses is equivalent to assessing whether random eﬀects can be excluded from the model by eﬀectively setting certain elements of Σ to 0, a diﬃcult problem given the constraints on Σ. Potentially, one can select a preferred GLMM by repeatedly ﬁtting the model for all possible choices of xit and zit and then applying a standard criterion, such as the AIC or BIC. Such an approach is not feasible unless the number of candidate predictors is modest, and there is no general consensus on what the penalty for model complexity should be in a model with random eﬀects. In addition, it is often not enough to say which model is preferred, one wants to report a weight of evidence (e.g., that the eﬀect of smoking on fecundability varies among women). To assess whether one or more random eﬀects should be included in the model, several authors have proposed frequentist score tests (Commenges and Jacqmin-Gadda, 1997; Lin, 1997; Hall and Praestgaard, 2001). In the Bayesian literature, Albert and Chib (1997) proposed an approach for testing whether a random intercept should be included, Sinharay and Stern (2001) developed a more general approach for calculating Bayes factors for comparing GLMMs, and Chen et al. (2003) proposed a class of informative priors for model selection in GLMMs. These methods focus on comparing two models at a time, and do not provide a general approach for searching for promising subsets of candidate predictors. In the setting of linear mixed models for normal data, Chen and Dunson (2003) proposed a Bayesian approach for random eﬀects selection based on using variable selection priors for the components in a special decomposition of the random eﬀects covariance. Related approaches have been used in graphical (or covariance structure) modeling for multivariate normal data (refer to Wong, Carter and Kohn, 2003; Liechty, Liechty and Muller, 2004 for recent references). Bayesian variable selection in conventional GLMs has also received a lot of interest in the literature. Raftery (1996) proposed an approximate Bayes factor approach, Meyer and Laud (2002) considered predictive variable selection, Nott and Leonte (2004) developed an innovative sampling algorithm and Ntzoufras, Dellaportas, and Foster (2003) developed methods for joint variable and link selection.

r2glmm 0.1.2 软件包文档说明书

Package‘r2glmm’October14,2022Type PackageTitle Computes R Squared for Mixed(Multilevel)ModelsDate2017-08-04Version0.1.2Description The model R squared and semi-partial R squared for the linear and generalized linear mixed model(LMM and GLMM)are computed with conﬁdence limits.The R squared measure from Edwards et.al(2008)<DOI:10.1002/sim.3429> is extended to the GLMM using penalized quasi-likelihood(PQL)estimation(see Jaeger et al.2016<DOI:10.1080/02664763.2016.1193725>).Three methods of computation are provided and described as follows.First,TheKenward-Roger approach.Due to some inconsistency between the'pbkrtest'package and the'glmmPQL'function,the Kenward-Roger approach in the'r2glmm'package is limited to the LMM.Second,The method introducedby Nakagawa and Schielzeth(2013)<DOI:10.1111/j.2041-210x.2012.00261.x>and later extended by Johnson(2014)<DOI:10.1111/2041-210X.12225>.The'r2glmm'package only computes marginal R squared for the LMM and doesnot generalize the statistic to the GLMM;however,conﬁdence limits andsemi-partial R squared forﬁxed effects are useful stly,anapproach using standardized generalized variance(SGV)can be used forcovariance model selection.Package installation instructions can be foundin the readmeﬁle.Imports mgcv,lmerTest,Matrix,pbkrtest,ggplot2,afex,stats,MASS,gridExtra,grid,data.table,dplyrSuggests lme4,nlme,testthatLicense GPL-2LazyData TRUERoxygenNote6.0.1URL https:///bcjaeger/r2glmmBugReports https:///bcjaeger/r2glmm/issuesNeedsCompilation noAuthor Byron Jaeger[aut,cre]12calc_sgvMaintainer Byron Jaeger<**********************>Repository CRANDate/Publication2017-08-0510:26:17UTCR topics documented:calc_sgv (2)cmp_R2 (3)glmPQL (4)pSym (5)make.partial.C (5)plot.R2 (6)pqlmer (7)print.R2 (8)r2beta (8)r2dt (10)Index12 calc_sgv Compute the standardized generalized variance(SGV)of a blockeddiagonal matrix.DescriptionCompute the standardized generalized variance(SGV)of a blocked diagonal matrix.Usagecalc_sgv(nblocks=NULL,blksizes=NULL,vmat)Argumentsnblocks Number of blocks in the matrix.blksizes vector of block sizesvmat The blocked covariance matrixValueThe SGV of the covariance matrix vmat.Exampleslibrary(Matrix)v1=matrix(c(1,0.5,0.5,1),nrow=2)v2=matrix(c(1,0.2,0.1,0.2,1,0.3,0.1,0.3,1),nrow=3)v3=matrix(c(1,0.1,0.1,0.1,1,0.2,0.1,0.2,1),nrow=3)calc_sgv(nblocks=3,blksizes=c(2,3,3),vmat=Matrix::bdiag(v1,v2,v3))cmp_R23 cmp_R2Compute R2with a speciﬁed C matrixDescriptionCompute R2with a speciﬁed C matrixUsagecmp_R2(c,x,SigHat,beta,method,obsperclust=NULL,nclusts=NULL)Argumentsc Contrast matrix forﬁxed effectsx Fixed effects design matrixSigHat estimated model covariance(matrix or scalar)betaﬁxed effects estimatesmethod the method for computing r2betaobsperclust number of observations per cluster(i.e.subject)nclusts number of clusters(i.e.subjects)ValueA vector with the Wald statistic(ncp),approximate Wald F statistic(F),numerator degrees of free-dom(v1),denominator degrees of freedom(v2),and the speciﬁed r squared value(Rsq)Exampleslibrary(nlme)library(lme4)library(mgcv)lmemod=lme(distance~age*Sex,random=~1|Subject,data=Orthodont)X=model.matrix(lmemod,data=Orthodont)SigHat=extract.lme.cov(lmemod,data=Orthodont)beta=fixef(lmemod)p=length(beta)obsperclust=as.numeric(table(lmemod$data[, Subject ]))nclusts=length(obsperclust)C=cbind(rep(0,p-1),diag(p-1))partial.c=make.partial.C(p-1,p,2)cmp_R2(c=C,x=X,SigHat=SigHat,beta=beta,obsperclust=obsperclust,nclusts=nclusts,method= sgv )cmp_R2(c=partial.c,x=X,SigHat=SigHat,beta=beta,obsperclust=obsperclust,nclusts=nclusts,method= sgv )4glmPQLglmPQL Compute PQL estimates forﬁxed effects from a generalized linearmodel.DescriptionCompute PQL estimates forﬁxed effects from a generalized linear model.UsageglmPQL(glm.mod,niter=20,data=NULL)Argumentsglm.mod a generalized linear modelﬁtted with the glm function.niter maximum number of iterations allowed in the PQL algorithm.data The data used by theﬁtted model.This argument is required for models with special expressions in their formula,such as offset,log,cbind(sucesses,trials),etc.ValueA glmPQL object(i.e.a linear model using pseudo outcomes).Examples#Load the datasets package for example codelibrary(datasets)library(dplyr)#We ll model the number of world changing discoveries per year for the#last100years as a poisson outcome.First,we set up the datadat=data.frame(discoveries)%>%mutate(year=1:length(discoveries))#Fit the GLM with a poisson link functionmod<-glm(discoveries~year+I(year^2),family= poisson ,data=dat)#Find PQL estimates using the original GLMmod.pql=glmPQL(mod)#Note that the PQL model yields a higher R Squared statistic#than the fit of a strictly linear model.This is attributed#to correctly modelling the distribution of outcomes and then#linearizing the model to measure goodness of fit,rather than#simply fitting a linear modelsummary(mod.pql)pSym5summary(linfit<-lm(discoveries~year+I(year^2),data=dat))r2beta(mod.pql)r2beta(linfit)pSym Checks if a matrix is Compound Symmetric.DescriptionChecks if a matrix is Compound Symmetric.UsagepSym(mat,tol=1e-05)Argumentsmat The matrix to be tested.tol a number indicating the smallest acceptable difference between off diagonal val-ues.ValueTrue if the matrix is compound symmetric.Examplesgcmat<-matrix(c(1,0.2,0.1,0.2,1,0.3,0.1,0.3,1),nrow=3)csmat<-matrix(c(1,0.2,0.2,0.2,1,0.2,0.2,0.2,1),nrow=3)pSym(csmat)make.partial.C Generate partial contrast matricesDescriptionGenerate partial contrast matricesUsagemake.partial.C(rows,cols,index)6plot.R2Argumentsrows Number of rows in the contrast matrixcols Number of columns in the contrast matrixindex A number corresponding to the position of theﬁxed effect in the vector ofﬁxed effect parameter estimates.ValueA contrast matrix designed to test theﬁxed effect corresponding to index in the vector ofﬁxedeffects.Examplesmake.partial.C(4,5,2)make.partial.C(4,5,3)make.partial.C(4,5,2:4)plot.R2Visualize standardized effect sizes and model R squaredDescriptionVisualize standardized effect sizes and model R squaredUsage##S3method for class R2plot(x,y=NULL,txtsize=10,maxcov=3,r2labs=NULL,r2mthd="sgv",cor=TRUE,...)Argumentsx An R2object from the r2beta function.y An R2object from the r2beta function.txtsize The text size of the axis labels.maxcov Maximum number of covariates to include in the semi-partial plots.r2labs a character vector containing labels for the models.The labels are printed as subscripts on a covariance model matrix.r2mthd The method used to compute R2cor An argument to be passed to the r2dt function.Only relevant if comparing two R2objects....Arguments to be passed to plotpqlmer7ValueA visual representation of the model and semi-partial R squared from the r2object provided.Exampleslibrary(nlme)library(r2glmm)data(Orthodont)#Linear mixed modellmemod=lme(distance~age*Sex,random=~1|Subject,data=Orthodont)r2=r2beta(model=lmemod,partial=TRUE,method= sgv )plot(x=r2)pqlmer pqlmerDescriptionFit a GLMM model with multivariate normal random effects using Penalized Quasi-Likelihood for mermod objects.Usagepqlmer(formula,family,data,niter=40,verbose=T)Argumentsformula The lme4model formula.family a family function of the error distribution and link function to be used in the model.data the dataframe containing the variables in the model.niter Maximum number of iterations to perform.verbose if TRUE,iterations are printed to console.ValueA pseudo linear mixed model of class"lme".See AlsoglmmPQLExamples#Compare lmer PQL with lme PQLlibrary(MASS)lmePQL=glmmPQL(y~trt+week+I(week>2),random=~1|ID,family=binomial,data=bacteria,verbose=FALSE)merPQL=pqlmer(y~trt+week+I(week>2)+(1|ID),family=binomial,data=bacteria,verbose=FALSE)summary(lmePQL)summary(merPQL)print.R2Print the contents of an R2objectDescriptionPrint the contents of an R2objectUsage##S3method for class R2print(x,...)Argumentsx an object of class R2...other arguments passed to the print function.r2beta r2beta Compute R Squared for Mixed ModelsDescriptionComputes coefﬁcient of determination(R squared)from edwards et al.,2008and the generalized R squared from Jaeger et al.,2016.Currently implemented for linear mixed models with lmer and lme objects.For generalized linear mixed models,only glmmPQL are supported.Usager2beta(model,partial=TRUE,method="sgv",data=NULL)Argumentsmodel aﬁtted mermod,lme,or glmmPQL model.partial if TRUE,semi-partial R squared are calculated for eachﬁxed effect in the mixed model.method Speciﬁes the method of computation for R squared beta:if method=’sgv’then the standardized generalized variance approach is applied.This method is rec-ommended for covariance model selection.if method=’kr’,then the KenwardRoger approach is applied.This option is only available for lme models.ifmethod=’nsj’,then the Nakagawa and Schielzeth approach is applied.This op-tion is available for lmer and lme objects.if method=’lm’,the classical Rsquared from the linear model is computed.This method should only be usedon glm and lm object.data The data used by theﬁtted model.This argument is required for models with special expressions in their formula,such as offset,log,cbind(sucesses,trials),etc.ValueA dataframe containing the model F statistic,numerator and denominator degrees of freedom,non-centrality parameter,and R squared statistic with95If partial=TRUE,then the dataframe also contains partial R squared statistics for allﬁxed effects in the model.ReferencesEdwards,Lloyd J.,et al."An R2statistic forﬁxed effects in the linear mixed model."Statistics in medicine27.29(2008):6137-6157.Nakagawa,Shinichi,and Holger Schielzeth."A general and simple method for obtaining R2from generalized linear mixed effects models."Methods in Ecology and Evolution4.2(2013):133-142.Jaeger,Byron C.,et al.,"An R Squared Statistic for Fixed Effects in the Generalized Linear Mixed Model."Journal of Applied Statistics(2016).Exampleslibrary(nlme)library(lme4)data(Orthodont)#Linear mixed modelsmermod=lmer(distance~age*Sex+(1|Subject),data=Orthodont)lmemod=lme(distance~age*Sex,random=~1|Subject,data=Orthodont)#The Kenward-Roger approachr2beta(mermod,method= kr )#Standardized Generalized Variancer2beta(mermod,method= sgv )r2beta(lmemod,method= sgv )10r2dt#The marginal R squared by Nakagawa and Schielzeth(extended by Johnson)r2beta(mermod,method= nsj )#linear and generalized linear modelslibrary(datasets)dis=data.frame(discoveries)dis$year=1:nrow(dis)lmod=lm(discoveries~year+I(year^2),data=dis)glmod=glm(discoveries~year+I(year^2),family= poisson ,data=dis)#Using an inappropriate link function(normal)leads to#a poor fit relative to the poisson link function.r2beta(lmod)r2beta(glmod)#PQL models#Currently only SGV method is supportedlibrary(MASS)PQL_bac=glmmPQL(y~trt+I(week>2),random=~1|ID,family=binomial,data=bacteria,verbose=FALSE)r2beta(PQL_bac,method= sgv )r2dt R Squared Difference Test(R2DT).Test for a statistically signiﬁcantdifference in generalized explained variance between two candidatemodels.DescriptionR Squared Difference Test(R2DT).Test for a statistically signiﬁcant difference in generalized ex-plained variance between two candidate models.Usager2dt(x,y=NULL,cor=TRUE,fancy=FALSE,onesided=TRUE,clim=95,nsims=2000,mu=NULL)Argumentsx An R2object from the r2beta function.y An R2object from the r2beta function.If y is not speciﬁed,Ho:E[x]=mu is tested(mu is speciﬁed by the user).r2dt11 cor if TRUE,the R squared statistics are assumed to be positively correlated anda simulation based approach is used.If FALSE,the R squared are assumedindependent and the difference of independent beta distributions is used.Thisonly needs to be speciﬁed when two R squared measures are being considered.fancy if TRUE,the output values are rounded and changed to characters.onesided if TRUE,the alternative hypothesis is that one model explains a larger propor-tion of generalized variance.If false,the alternative is that the amount of gener-alized variance explained by the two candidate models is not equal.clim Desired conﬁdence level for interval estimates regarding the difference in gen-eralized explained variance.nsims number of samples to draw when simulating correlated non-central beta random variables.This parameter is only relevant if cor=TRUE.mu Used to test Ho:E[x]=mu.ValueA conﬁdence interval for the difference in R Squared statistics and a p-value corresponding to thenull hypothesis of no difference.Exampleslibrary(nlme)library(lme4)library(r2glmm)data(Orthodont)#Comparing two linear mixed modelsm1=lmer(distance~age*Sex+(1|Subject),Orthodont)m2=lmer(distance~age*Sex+(1+age|Subject),Orthodont)m1r2=r2beta(model=m1,partial=FALSE)m2r2=r2beta(model=m2,partial=FALSE)#Accounting for correlation can make a substantial difference.r2dt(x=m1r2,y=m2r2,cor=TRUE)r2dt(x=m1r2,y=m2r2,cor=FALSE)Indexcalc_sgv,2cmp_R2,3glmmPQL,7,8glmPQL,4pSym,5lme,8,9lmer,8,9make.partial.C,5plot.R2,6pqlmer,7print.R2,8r2beta,8r2dt,1012。

Generalizedlinearmodels：广义线性模型

Generalized linear modelsDouglas BatesNovember01,2010Contents1Deﬁnition1 2Links2 3Estimating parameters5 4Example6 5Model building8 6Conclusions8 7Summary9 1Generalized Linear ModelsGeneralized Linear Models•When using linear models(LMs)we assume that the response being modeled is on a continuous scale.•Sometimes we can bend this assumption a bit if the response is an ordinal response witha moderate to large number of levels.For example,the Scottish secondary school testresults in the mlmRev package are integer values on the scale of1to10but we analyze them on a continuous scale.•However,an LM is not suitable for modeling a binary response,an ordinal response with few levels or a response that represents a count.For these we use generalized linear models(GLMs).•To describe GLMs we return to the representation of the response as an n-dimensional, vector-valued,random variable,Y.Parts of LMs carried over to GLMs•Random variables–Y the response variable•Parameters–β-ﬁxed-eﬀects coeﬃcients–σ-a scale parameter(not always used)•The linear predictor Xβwhere–X is the n×p model matrix forβThe probability model•For GLMs we retain some of the properties of the LM probability modelY∼NµY,σ2IwhereµY=XβSpeciﬁcally–The distribution of Y depends onβonly through the mean,µY=Xβ.–Elements of Y are independent.That is,the distribution of Y is completely speciﬁed by the univariate distributions,Y i,i=1,...,n.–These univariate,distributions all have the same form.They diﬀer only in their means.•GLMs diﬀer from LMs in the form of the univariate distributions and in howµY depends on the linear predictor,Xβ.2Speciﬁc distributions and linksSome choices of univariate distributions•Typical choices of univariate distributions are:–The Bernoulli distribution for binary(0/1)data,which has probability mass func-tionp(y|µ)=µy(1−µ)1−y,0<µ<1,y=0,1–Several independent binary responses can be represented as a binomial response, but only if all the Bernoulli distributions have the same mean.–The Poisson distribution for count(0,1,...)data,which has probability mass func-tionp(y|µ)=e−µµyy!,0<µ,y=0,1,2,...•All of these distributions are completely speciﬁed by the mean.This is diﬀerent from the normal(or Gaussian)distribution,which also requires the scale parameter,σ.The link function,g•When the univariate distributions have constraints onµ,such as0<µ<1(Bernoulli) or0<µ(Poisson),we cannot deﬁne the mean,µY,to be equal to the linear predictor, Xβ,which is unbounded.•We choose an invertible,univariate link function,g,such thatη=g(µ)is unconstrained.The vector-valued link function,g,is deﬁned by applying g component-wise.η=g(µ)whereηi=g(µi),i=1,...,n•We require that g be invertible so thatµ=g−1(η)is deﬁned for−∞<η<∞and is in the appropriate range(0<µ<1for the Bernoulli or0<µfor the Poisson).The vector-valued inverse link,g−1,is deﬁned component-wise.“Canonical”link functions•There are many choices of invertible scalar link functions,g,that we could use for a given set of constraints.•For the Bernoulli and Poisson distributions,however,one link function arises naturally from the deﬁnition of the probability mass function.(The same is true for a few other, related but less frequently used,distributions,such as the gamma distribution.)•To derive the canonical link,we consider the logarithm of the probability mass function (or,for continuous distributions,the probability density function).•For distributions in this“exponential”family,the logarithm of the probability mass or density can be written as a sum of terms,some of which depend on the response,y,only and some of which depend on the mean,µ,only.However,only one term depends on both y andµ,and this term has the form y·g(µ),where g is the canonical link.The canonical link for the Bernoulli distribution•The logarithm of the probability mass function islog(p(y|µ))=log(1−µ)+y logµ1−µ,0<µ<1,y=0,1.•Thus,the canonical link function is the logit linkη=g(µ)=logµ1−µ.•Becauseµ=P[Y=1],the quantityµ/(1−µ)is the odds ratio(in the range(0,∞))and g is the logarithm of the odds ratio,sometimes called“log odds”.•The inverse link isµ=g−1(η)=eη1+eη=11+e−ηPlot of canonical link for the Bernoulli distributionµη=l o g (µ1−µ)−550.00.20.40.60.8 1.0Plot of inverse canonical link for the Bernoulli distributionηµ=11+e x p (−η)0.00.20.40.60.81.0−505Evaluating links for the Bernoulli•As a monotone increasing function the maps (−∞,∞)to (0,1)the inverse link,g −1,is a cumulative distribution function for a continuous random variable.The link,q ,is the corresponding quantile function.•The canonical link is the quantile function for the logistic distribution (R functions qlogis and plogis ).•In the past the"probit"link was sometimes used instead of the logit link.For this the link is the standard normal quantile and the inverse link is the standard normal cumulative,sometimes writtenΦ(z)(R functions qnorm and pnorm).•As described in R’s help page for the binomial functionthe binomial family(allows)the links"logit","probit","cauchit",(cor-responding to logistic,normal and Cauchy CDFs respectively)"log"and"cloglog"(complementary log-log)The canonical link for the Poisson distribution•The logarithm of the probability mass islog(p(y|µ))=log(y!)−µ+y log(µ)•Thus,the canonical link function for the Poisson is the log linkη=g(µ)=log(µ)•The inverse link isµ=g−1(η)=eηThe canonical link related to the variance•For the canonical link function,the derivative of its inverse is the variance of the response.•For the Bernoulli,the canonical link is the logit and the inverse link isµ=g−1(η)= 1/(1+e−η).Thendµdη=e−η(1+e−η)2=11+e−ηe−η1+e−η=µ(1−µ)=Var(Y)•For the Poisson,the canonical link is the log and the inverse link isµ=g−1(η)=eη.Thendµdη=eη=µ=Var(Y)3Estimating parametersEstimating parameters•We determine the maximum likelihood estimates,(mle’s),of the coeﬃcients,β,in the linear predictor and,if used,the scale parameter.•There is no direct algorithm for determining the mle’s in a GLM and we must use iterative algorithms.Fortunately,there is a very eﬀective iterative algorithm which is based on weighted least squares to update parameter estimates.The weights are inversely proportional to the variance of the response,but the variance depends on the mean which, in turn,depends on the parameters.•The IRLS (iteratively reweighted least squares)algorithm ﬁxes the weights,determines the parameter values that minimize the weighted sum of squared residuals,then updates the weights and repeats the process until the weights stabilize.•This algorithm converges very quickly.•The original description of IRLS from McCullagh and Nelder’s book is,I feel,overly complex.4Data description and initial explorationContraception data•One of the data sets in the "mlmRev"package,derived from data ﬁles available on the multilevel modelling web site,is from a fertility survey of women in Bangladesh.•One of the (binary)responses recorded is whether or not the woman currently uses artiﬁcial contraception.•Covariates included the woman’s age (on a centered scale),the number of live children she had,whether she lived in an urban or rural setting,and the district in which she lived.•Instead of plotting such data as points,we use the 0/1response to generate scatterplot smoother curves versus age for the diﬀerent groups.Contraception use versus age by urban and livchCentered ageP r o p o r t i o n 0.00.20.40.60.81.00123+Comments on the data plot•These observational data are unbalanced(some districts have only2observations,some have nearly120).They are not longitudinal(no“time”variable).•Binary responses have low per-observation information content(exactly one bit per ob-servation).Districts with few observations will not contribute strongly to estimates of random eﬀects.•Within-district plots will be too imprecise so we only examine the global eﬀects in plots.•The comparisons on the multilevel modelling site are forﬁts of a model that is linear in age,which is clearly inappropriate.•The form of the curves suggests at least a quadratic in age.•The urban versus rural diﬀerences may be additive.•It appears that the livch factor could be dichotomized into“0”versus“1or more”.Preliminary modelﬁtCall:glm(formula=use~age+I(age^2)+urban+livch,family=binomial,data=Contraception)Deviance Residuals:Min1Q Median3Q Max-1.4738-1.0369-0.6683 1.2401 1.9765Coefficients:Estimate Std.Error z value Pr(>|z|)(Intercept)-0.94995210.1560118-6.0891.14e-09age0.00458370.00890840.5150.607I(age^2)-0.00428650.0007002-6.1229.23e-10urbanY0.76809750.10619167.2334.72e-13livch10.78311280.1569096 4.9916.01e-07livch20.85490400.1783573 4.7931.64e-06livch3+0.80602510.1784817 4.5166.30e-06(Dispersion parameter for binomial family taken to be1)Null deviance:2590.9on1933degrees of freedomResidual deviance:2417.7on1927degrees of freedomAIC:2431.7Number of Fisher Scoring iterations:4Comments on the modelﬁt•There is a highly signiﬁcant quadratic term in age.•The linear term in age is not signiﬁcant but we retain it because the age scale has been centered at an arbitrary value(which,unfortunately,is not provided with the data).•The urban factor is highly signiﬁcant(as indicated by the plot).•Levels of livch greater than0are signiﬁcantly diﬀerent from0but may not be diﬀerent from each other.5Model buildingReduced model with dichotomized livchCall:glm(formula=use~age+I(age^2)+urban+ch,family=binomial,data=Contraception)Deviance Residuals:Min1Q Median3Q Max-1.4544-1.0371-0.6682 1.2378 1.9790Coefficients:Estimate Std.Error z value Pr(>|z|)(Intercept)-0.94293550.1497300-6.2983.02e-10age0.00496750.00758890.6550.513I(age^2)-0.00432750.0006918-6.2553.97e-10urbanY0.76636160.10594637.2334.71e-13chY0.80621720.1422619 5.6671.45e-08(Dispersion parameter for binomial family taken to be1)Null deviance:2590.9on1933degrees of freedomResidual deviance:2417.9on1929degrees of freedomAIC:2427.9Number of Fisher Scoring iterations:4Comparing the modelﬁts•A likelihood ratio test can be used to compare these nested models.>anova(cm2,cm1,test="Chisq")Analysis of Deviance TableModel1:use~age+I(age^2)+urban+chModel2:use~age+I(age^2)+urban+livchResid.Df Resid.Dev Df Deviance P(>|Chi|)119292417.9219272417.720.205250.9025•The large p-value indicates that we would not reject cm2in favor of cm1hence we prefer the more parsimonious cm2.•The plot of the scatterplot smoothers according to live children or none indicates that there may be a diﬀerence in the age pattern between these two groups.6Conclusions from the exampleConclusions from the example•Carefully plotting the data is enormously helpful in formulating the model.•Observational data tend to be unbalanced and have many more covariates than data from a designed experiment.Formulating a model is typically more diﬃcult than in a designed experiment.•A generalized linear model isﬁt with the function glm()which requires a family ar-gument.Typical values are binomial or poisson.By default,the family will use the canonical link.Other links can be speciﬁed as,e.g.binomial(link="cloglog").•We use likelihood-ratio tests in model building.The z-tests provided in the model sum-mary provide a good indication of the results but should be conﬁrmed byﬁtting the reduced model and performing a likelihood-ratio test.A word about overdispersion•In many application areas using“pseudo”distribution families,such as quasibinomial and quasipoisson,is a popular and well-accepted technique for accomodating variability that is apparently larger than would be expected from a binomial or a Poisson distribu-tion.•This amounts to adding an extra parameter,likeσ,the common scale parameter in a LM,to the distribution of the response.•It is possible to form an estimate of such a quantity during the IRLS algorithm but it is an artiﬁcial construct.There is no probability distribution with such a parameter.7SummarySummary•GLMs allow for the distribution of Y to be other than a Gaussian.A Bernoulli(or, more generally,a binomial)distribution is used to model binary or binomial responses.A Poisson distribution is used to model responses that are counts.•The mean depends upon the linear predictor,Xβ,through the inverse link function, g−1.•The mle’s of the parameters are determined through iteratively reweighted least squares (IRLS).•“Wald-type”tests on the coeﬃcients are displayed in the summary but to be more con-ﬁdent of the results we shouldﬁt both models and use a comparative anova with the optional argument,test="Chisq".。

generalized linear model结果解释-概述说明以及解释

generalized linear model结果解释-概述说明以及解释1.引言1.1 概述概述部分的内容可以包括对广义线性模型的简要介绍以及结果解释的重要性。

以下是一种可能的编写方式：在统计学和机器学习领域，广义线性模型（Generalized Linear Model，简称GLM）是一种常用的统计模型，用于建立因变量与自变量之间的关系。

与传统的线性回归模型不同，广义线性模型允许因变量（也称为响应变量）的分布不服从正态分布，从而更适用于处理非正态分布的数据。

广义线性模型的理论基础是广义线性方程（Generalized Linear Equation），它通过引入连接函数（Link Function）和系统误差分布（Error Distribution）的概念，从而使模型能够适应不同类型的数据。

结果解释是广义线性模型分析中的一项重要任务。

通过解释模型的结果，我们可以深入理解自变量与因变量之间的关系，并从中获取有关影响因素的信息。

结果解释能够帮助我们了解自变量的重要性、方向性及其对因变量的影响程度。

通过对结果进行解释，我们可以推断出哪些因素对于观察结果至关重要，从而对问题的本质有更深入的认识。

本文将重点讨论如何解释广义线性模型的结果。

我们将介绍广义线性模型的基本概念和原理，并指出结果解释中需要注意的要点。

此外，我们将提供实际案例和实例分析，以帮助读者更好地理解结果解释的方法和过程。

通过本文的阅读，读者将能够更全面地了解广义线性模型的结果解释，并掌握解释结果的相关技巧和方法。

本文的目的是帮助读者更好地理解和运用广义线性模型，从而提高统计分析和机器学习的能力。

在接下来的章节中，我们将详细介绍广义线性模型及其结果解释的要点，希望读者能够从中受益。

1.2文章结构文章结构部分的内容应该是对整篇文章的结构进行简要介绍和概述。

这个部分通常包括以下内容：文章结构部分的内容：本文共分为引言、正文和结论三个部分。

其中，引言部分主要概述了广义线性模型的背景和重要性，并介绍了文章的目的。

glmmfields 0.1.8 用户手册说明书

Package‘glmmﬁelds’October20,2023Type PackageTitle Generalized Linear Mixed Models with Robust Random Fields for Spatiotemporal ModelingVersion0.1.8Description Implements Bayesian spatial and spatiotemporalmodels that optionally allow for extreme spatial deviations throughtime.'glmmﬁelds'uses a predictive process approach with randomﬁelds implemented through a multivariate-t distribution instead ofthe usual multivariate normal.Sampling is conducted with'Stan'.References:Anderson and Ward(2019)<doi:10.1002/ecy.2403>. License GPL(>=3)URL https:///seananderson/glmmfieldsBugReports https:///seananderson/glmmfields/issues Depends methods,R(>=3.4.0),Rcpp(>=0.12.18)Imports assertthat,broom,broom.mixed,cluster,dplyr(>=0.8.0), forcats,ggplot2(>=2.2.0),loo(>=2.0.0),mvtnorm,nlme,RcppParallel(>=5.0.1),reshape2,rstan(>=2.26.0),rstantools(>=2.1.1),tibbleSuggests bayesplot,coda,knitr,parallel,rmarkdown,testthat,viridisLinkingTo BH(>=1.66.0),Rcpp(>=0.12.8),RcppEigen(>=0.3.3.3.0), RcppParallel(>=5.0.1),rstan(>=2.26.0),StanHeaders(>=2.26.0)VignetteBuilder knitrEncoding UTF-8RoxygenNote7.2.3SystemRequirements GNU makeNeedsCompilation yesBiarch true12glmmﬁelds-package Author Sean C.Anderson[aut,cre],Eric J.Ward[aut],Trustees of Columbia University[cph]Maintainer Sean C.Anderson<********************>Repository CRANDate/Publication2023-10-2017:50:02UTCR topics documented:glmmﬁelds-package (2)format_data (3)glmmﬁelds (4)lognormal (7)loo.glmmﬁelds (8)nbinom2 (8)plot.glmmﬁelds (9)predict (10)sim_glmmﬁelds (12)stan_pars (13)student_t (14)tidy (14)Index16 glmmfields-package The’glmmﬁelds’package.DescriptionImplements Bayesian spatial and spatiotemporal models that optionally allow for extreme spatial deviations through time.’glmmﬁelds’uses a predictive process approach with randomﬁelds im-plemented through a multivariate-t distribution instead of the usual multivariate normal.Sampling is conducted with’Stan’.ReferencesStan Development Team(2018).RStan:the R interface to Stan.R package version2.18.2.format_data3 format_data Format data forﬁtting a glmmﬁelds modelDescriptionFormat data forﬁtting a glmmﬁelds modelUsageformat_data(data,y,X,time,lon="lon",lat="lat",station=NULL,nknots=25L,covariance=c("squared-exponential","exponential","matern"),fixed_intercept=FALSE,cluster=c("pam","kmeans"))Argumentsdata A data frame to be formattedy A numeric vector of the responseX A matrix of the predictorstime A character object giving the name of the time columnlon A character object giving the name of the longitude columnlat A character object giving the name of the latitude columnstation A numeric vector giving the integer ID of the stationnknots The number of knotscovariance The type of covariance functionfixed_interceptShould the intercept beﬁxed?cluster The type of clustering algorithm used to determine the not locations."pam"= pam.kmeans is faster for large datasets.glmmfields Fit a spatiotemporal randomﬁelds GLMMDescriptionFit a spatiotemporal randomﬁelds model that optionally uses the MVT distribution instead of a MVN distribution to allow for spatial extremes through time.It is also possible toﬁt a spatial randomﬁelds model without a time component.Usageglmmfields(formula,data,lon,lat,time=NULL,nknots=15L,prior_gp_theta=half_t(3,0,5),prior_gp_sigma=half_t(3,0,5),prior_sigma=half_t(3,0,5),prior_rw_sigma=half_t(3,0,5),prior_intercept=student_t(3,0,10),prior_beta=student_t(3,0,3),prior_phi=student_t(1000,0,0.5),fixed_df_value=1000,fixed_phi_value=0,estimate_df=FALSE,estimate_ar=FALSE,family=gaussian(link="identity"),binomial_N=NULL,covariance=c("squared-exponential","exponential","matern"),matern_kappa=0.5,algorithm=c("sampling","meanfield"),year_re=FALSE,nb_lower_truncation=0,control=list(adapt_delta=0.9),save_log_lik=FALSE,df_lower_bound=2,cluster=c("pam","kmeans"),offset=NULL,...)Argumentsformula The model formula.data A data frame.lon A character object giving the name of the longitude column.lat A character object giving the name of the latitude column.time A character object giving the name of the time column.Leave as NULL toﬁt a spatial GLMM without a time element.nknots The number of knots to use in the predictive process model.Smaller values will be faster but may not adequately represent the shape of the spatial pattern. prior_gp_theta The prior on the Gaussian Process scale parameter.Must be declared with half_t().Here,and throughout,priors that are normal or half-normal can beimplemented by setting theﬁrst parameter in the half-t or student-t distributionto a large value.E.g.something greater than100.prior_gp_sigma The prior on the Gaussian Process eta parameter.Must be declared with half_t(). prior_sigma The prior on the observation process scale parameter.Must be declared with half_t().This acts as a substitute for the scale parameter in whatever obser-vation distribution is being used.E.g.the CV for the Gamma or the dispersionparameter for the negative binomial.prior_rw_sigma The prior on the standard deviation parameter of the random walk process(if speciﬁed).Must be declared with half_t().prior_interceptThe prior on the intercept parameter.Must be declared with student_t(). prior_beta The prior on the slope parameters(if any).Must be declared with student_t(). prior_phi The prior on the AR parameter.Must be declared with student_t().fixed_df_value Theﬁxed value for the student-t degrees of freedom parameter if the degrees of freedom parameter isﬁxed in the MVT.If the degrees of freedom parameteris estimated then this argument is ignored.Must be1or greater.Very largevalues(e.g.the default value)approximate the normal distribution.If the valueis>=1000then a true MVN distribution will beﬁt.fixed_phi_valueTheﬁxed value for temporal autoregressive parameter,between randomﬁeldsat time(t)and time(t-1).If the phi parameter is estimated then this argument isignored.estimate_df Logical:should the degrees of freedom parameter be estimated?estimate_ar Logical:should the AR(autoregressive)parameter be estimated?Here,this refers to a autoregressive process in the evolution of the spatialﬁeld throughtime.family Family object describing the observation model.Note that only one link is implemented for each distribution.Gamma,negative binomial(speciﬁed vianbinom2()as nbinom2(link="log"),and Poisson must have a log link.Bi-nomial must have a logit link.Also implemented is the lognormal(speciﬁedvia lognormal()as lognormal(link="log").Besides the negative binomialand lognormal,other families are speciﬁed as shown in family.binomial_N A character object giving the optional name of the column containing Binomial sample size.Leave as NULL toﬁt a spatial GLMM with sample sizes(N)=1,equivalent to bernoulli model.covariance The covariance function of the Gaussian Process.One of"squared-exponential", "exponential",or"matern".matern_kappa Optional parameter for the Matern covariance function.Optional values are1.5 or2.5.Values of0.5are equivalent to exponential.algorithm Character object describing whether the model should beﬁt with full NUTS MCMC or via the variational inference mean-ﬁeld approach.See rstan::vb().Note that the variational inference approach should not be trusted forﬁnal infer-ence and is much more likely to give incorrect inference than MCMC.year_re Logical:estimate a random walk for the time variable?If TRUE,then noﬁxed effects(B coefﬁcients)will be estimated.In this case,prior_intercept willbe used as the prior for the initial value in time.nb_lower_truncationFor NB2only:lower truncation value. E.g.0for no truncation,1for1andall values above.Note that estimation is likely to be considerably slower withlower truncation because the sampling is not vectorized.Also note that thelog likelihood values returned for estimating quantities like LOOIC will not becorrect if lower truncation is implemented.control List to pass to rstan::sampling().For example,increase adapt_delta if there are warnings about divergent transitions:control=list(adapt_delta=0.99).By default,glmmﬁelds sets adapt_delta=0.9.save_log_lik Logical:should the log likelihood for each data point be saved so that informa-tion criteria such as LOOIC or W AIC can be calculated?Defaults to FALSE sothat the size of model objects is smaller.df_lower_bound The lower bound on the degrees of freedom parameter.Values that are too low,e.g.below2or3,it might affect chain convergence.Defaults to2.cluster The type of clustering algorithm used to determine the knot locations."pam"= cluster::pam().The"kmeans"algorithm will be faster on larger datasets.offset An optional offset vector....Any other arguments to pass to rstan::sampling().DetailsNote that there is no guarantee that the default priors are reasonable for your data.Also,there is no guarantee the default priors will remain the same in future versions.Therefore it is important that you specify any priors that are used in your model,even if they replicate the defaults in the package.It is particularly important that you consider that prior on gp_theta since it depends on the distance between your location points.You may need to scale your coordinate units so they are on a ballpark range of1-10by,say,dividing the coordinates(say in UTMs)by several order of magnitude.Examples#Spatiotemporal example:set.seed(1)s<-sim_glmmfields(n_draws=12,n_knots=12,gp_theta=1.5,gp_sigma=0.2,sd_obs=0.2)lognormal7 print(s$plot)#options(mc.cores=parallel::detectCores())#for parallel processing#should use4or more chains for real model fitsm<-glmmfields(y~0,time="time",lat="lat",lon="lon",data=s$dat,nknots=12,iter=1000,chains=2,seed=1)#Spatial example(with covariates)from the vignette and customizing#some priors:set.seed(1)N<-100#number of data pointstemperature<-rnorm(N,0,1)#simulated temperature dataX<-cbind(1,temperature)#design matrixs<-sim_glmmfields(n_draws=1,gp_theta=1.2,n_data_points=N,gp_sigma=0.3,sd_obs=0.1,n_knots=12,obs_error="gamma",covariance="squared-exponential",X=X,B=c(0.5,0.2))#B represents our intercept and sloped<-s$datd$temperature<-temperaturelibrary(ggplot2)ggplot(s$dat,aes(lon,lat,colour=y))+viridis::scale_colour_viridis()+geom_point(size=3)m_spatial<-glmmfields(y~temperature,data=d,family=Gamma(link="log"),lat="lat",lon="lon",nknots=12,iter=2000,chains=2,prior_beta=student_t(100,0,1),prior_intercept=student_t(100,0,5),control=list(adapt_delta=0.95))lognormal Lognormal familyDescriptionLognormal familyUsagelognormal(link="log")Argumentslink The link(must be log)Exampleslognormal()8nbinom2 loo.glmmfields Return LOO information criteriaDescriptionExtract the LOOIC(leave-one-out information criterion)using loo::loo().Usage##S3method for class glmmfieldsloo(x,...)Argumentsx Output from glmmfields().Must beﬁt with save_log_lik=TRUE,which is not the default....Arguments for loo::relative_eff()and loo::loo.array().Examplesset.seed(1)s<-sim_glmmfields(n_draws=12,n_knots=12,gp_theta=1.5,gp_sigma=0.2,sd_obs=0.2)#options(mc.cores=parallel::detectCores())#for parallel processing#save_log_lik defaults to FALSE to save space but is needed for loo():m<-glmmfields(y~0,time="time",lat="lat",lon="lon",data=s$dat,nknots=12,iter=1000,chains=4,seed=1,save_log_lik=TRUE)loo(m)nbinom2Negative binomial familyDescriptionThis is the NB2parameterization where the variance scales quadratically with the mean.Usagenbinom2(link="log")plot.glmmﬁelds9Argumentslink The link(must be log)Examplesnbinom2()plot.glmmfields Plot predictions from an glmmﬁelds modelDescriptionPlot predictions from an glmmﬁelds modelUsage##S3method for class glmmfieldsplot(x,type=c("prediction","spatial-residual","residual-vs-fitted"),link=TRUE,...)Argumentsx An object returned by glmmfieldstype Type of plotlink Logical:should the plots be made on the link scale or on the natural scale?...Other arguments passed to predict.glmmfieldsExamples#Spatiotemporal example:set.seed(1)s<-sim_glmmfields(n_draws=12,n_knots=12,gp_theta=1.5,gp_sigma=0.2,sd_obs=0.1)#options(mc.cores=parallel::detectCores())#for parallel processingm<-glmmfields(y~0,time="time",lat="lat",lon="lon",data=s$dat,nknots=12,iter=600,chains=1)x<-plot(m,type="prediction")xx+ggplot2::scale_color_gradient2()plot(m,type="spatial-residual")plot(m,type="residual-vs-fitted")10predict predict Predict from a glmmﬁelds modelDescriptionThese functions extract posterior draws or credible intervals.The helper functions are named to match those in the rstanarm package and call the function predict()with appropriate argument values.Usage##S3method for class glmmfieldspredictive_interval(object,...)##S3method for class glmmfieldsposterior_linpred(object,...)##S3method for class glmmfieldsposterior_predict(object,...)##S3method for class glmmfieldspredict(object,newdata=NULL,estimate_method=c("median","mean"),conf_level=0.95,interval=c("confidence","prediction"),type=c("link","response"),return_mcmc=FALSE,offset=NULL,iter="all",...)Argumentsobject An object returned by glmmfields()....Ignored currentlynewdata Optionally,a data frame to predict onestimate_methodMethod for computing point estimate("mean"or"median") conf_level Probability level for the credible intervals.interval Type of interval calculation.Same as for stats::predict.lm().type Whether the predictions are returned on"link"scale or"response"scale(Same as for stats::predict.glm()).predict11return_mcmc Logical.Should the full MCMC draws be returned for the predictions?offset Optional offset vector to be used in prediction.iter Number of MCMC iterations to draw.Defaults to all.Exampleslibrary(ggplot2)#simulate:set.seed(1)s<-sim_glmmfields(n_draws=12,n_knots=12,gp_theta=2.5,gp_sigma=0.2,sd_obs=0.1)#fit:#options(mc.cores=parallel::detectCores())#for parallel processingm<-glmmfields(y~0,data=s$dat,time="time",lat="lat",lon="lon",nknots=12,iter=800,chains=1)#Predictions:#Link scale credible intervals:p<-predict(m,type="link",interval="confidence")head(p)#Prediction intervals on new observations(include observation error):p<-predictive_interval(m)head(p)#Posterior prediction draws:p<-posterior_predict(m,iter=100)dim(p)#rows are iterations and columns are data elements#Draws from the linear predictor(not in link space):p<-posterior_linpred(m,iter=100)dim(p)#rows are iterations and columns are data elements#Use the tidy method to extract parameter estimates as a data frame:head(tidy(m,conf.int=TRUE,conf.method="HPDinterval"))#Make predictions on a fine-scale spatial grid:pred_grid<-expand.grid(lat=seq(min(s$dat$lat),max(s$dat$lat),length.out=25),lon=seq(min(s$dat$lon),max(s$dat$lon),length.out=25),time=unique(s$dat$time))pred_grid$prediction<-predict(m,newdata=pred_grid,type="response",iter=100,12sim_glmmﬁelds estimate_method="median",offset=rep(0,nrow(pred_grid)))$estimateggplot(pred_grid,aes(lon,lat,fill=prediction))+facet_wrap(~time)+geom_raster()+scale_fill_gradient2()sim_glmmfields Simulate a randomﬁeld with a MVT distributionDescriptionSimulate a randomﬁeld with a MVT distributionUsagesim_glmmfields(n_knots=15,n_draws=10,gp_theta=0.5,gp_sigma=0.2,mvt=TRUE,df=1e+06,seed=NULL,n_data_points=100,sd_obs=0.1,covariance=c("squared-exponential","exponential","matern"),matern_kappa=0.5,obs_error=c("normal","gamma","poisson","nb2","binomial","lognormal"),B=c(0),phi=0,X=rep(1,n_draws*n_data_points),g=data.frame(lon=runif(n_data_points,0,10),lat=runif(n_data_points,0,10)) )Argumentsn_knots The number of knotsn_draws The number of draws(for example,the number of years)gp_theta The Gaussian Process scale parametergp_sigma The Gaussian Process variance parametermvt Logical:MVT?(vs.MVN)df The degrees of freedom parameter for the MVT distributionseed The random seed valuestan_pars13n_data_points The number of data points per drawsd_obs The observation process scale parametercovariance The covariance function of the Gaussian process("squared-exponential","expo-nential","matern")matern_kappa The optional matern parameter.Can be1.5or2.5.Values of0.5equivalent to exponential model.obs_error The observation error distributionB A vector of parameters.Theﬁrst element is the interceptphi The auto regressive parameter on the mean of the randomﬁeld knotsX The model matrixg Grid of pointsExampless<-sim_glmmfields(n_draws=12,n_knots=12,gp_theta=1.5,gp_sigma=0.2,sd_obs=0.2)names(s)stan_pars Return a vector of parametersDescriptionReturn a vector of parametersUsagestan_pars(obs_error,estimate_df=TRUE,est_temporalRE=FALSE,estimate_ar=FALSE,fixed_intercept=FALSE,save_log_lik=FALSE)Argumentsobs_error The observation error distributionestimate_df Logical indicating whether the degrees of freedom parameter should be esti-matedest_temporalRE Logical:estimate a random walk for the time variable?estimate_ar Logical indicating whether the ar parameter should be estimated14tidy fixed_interceptShould the intercept beﬁxed?save_log_lik Logical:should the log likelihood for each data point be saved so that informa-tion criteria such as LOOIC or W AIC can be calculated?Defaults to FALSE sothat the size of model objects is smaller.student_t Student-t and half-t priorsDescriptionStudent-t and half-t priors.Note that this can be used to represent an effectively normal distribution prior by setting theﬁrst argument(the degrees of freedom parameter)to a large value(roughly50 or above).Usagestudent_t(df=3,location=0,scale=1)half_t(df=3,location=0,scale=1)Argumentsdf Degrees of freedom parameterlocation Location parameterscale Scale parameterExamplesstudent_t(3,0,1)half_t(3,0,1)tidy Tidy model outputDescriptionTidy model outputUsagetidy(x,...)##S3method for class glmmfieldstidy(x,...)tidy15Argumentsx Output from glmmfields()...Other argumentsIndexcluster::pam(),6family,5format_data,3glmmfields,4,9glmmfields(),8,10,15glmmfields-package,2half_t(student_t),14half_t(),5lognormal,7lognormal(),5loo(loo.glmmfields),8loo.glmmfields,8loo::loo(),8loo::loo.array(),8loo::relative_eff(),8nbinom2,8nbinom2(),5pam,3plot.glmmfields,9posterior_linpred(predict),10posterior_predict(predict),10predict,10predict.glmmfields,9predictive_interval(predict),10 rstan::sampling(),6rstan::vb(),6sim_glmmfields,12stan_pars,13stats::predict.glm(),10stats::predict.lm(),10student_t,14student_t(),5tidy,1416。

05 - Generalized Linear Models 广义线性模型

Lecture 5 Generalized Linear Models
Outline
Introduction Theory of the Generalized Linear Models
Logistic Regression
Poisson and Negative Binomial Regression
Introduction
Review: How to compare treatments?
Usually an endpoint is compared across treatment groups,
while controlling for important predictors
• Example: Control for baseline measurements • Predictors could be continuous or categorical
Relapse
Remission
10 | Basic Statistics in Clinical Trials | Generalized Linear Models | All Rights Reserved
Time
Introduction
RRMS study example – MRI scans
- T1 lesion - Combined unique active lesion (CUAL)
11 | Basic Statistics in Clinical Trials | Generalized Linear Models | All Rights Reserved
Introduction
Introduction

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 modiﬁcations 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 inefﬁcient 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 deﬁciencies of a linear model.In sec-tion3,I will present the generalized linear mixed model.In section4,I will present the estimation equations for theﬁxed 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 theﬁxed 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 difﬁcult 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 deﬁciencies 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.Aﬁnal option is to ignore the deﬁciencies 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.Speciﬁcally we have a set of estimation procedures which are based on a linear mixed model and manipulate the data to make itﬁt 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 extraﬂexibility 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 ﬁxed 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,theﬁxed 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),speciﬁes 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 whichﬁxed 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 theﬁxed 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 theﬁxed 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η.Theﬁrst 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 ﬁt 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 NEFFﬁxed 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 theﬁrst 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.Theﬁrst 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 modiﬁcations for a GLMM we will brieﬂy 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.2Modiﬁcations 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 difﬁculties associated with generalized linear mixed models are related to either estimates going to inﬁnity or divide by zero problems.Consider the uni-variate logit model.The calculations involve 1i i .Provided 0<µ<1this quantity is well deﬁned.If µapproaches either zero or one,then the quantity 1µi (1−µi )approaches inﬁnity.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 asﬁxed 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 inﬁnity.Asηi→±∞,W ii→0andhry i→n(y i−µi)In the limit W ii and hry i are both well deﬁned.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 sufﬁciently 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,theﬁxed effect which was approaching inﬁnity will be set to zero.This substitution can result in interesting convergence problems.A solution would be to“ﬁx”the offensiveﬁxed 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 ﬂexible way to model production traits which do not satisfy the assumptions of a linear mixed model.This ﬂexibility allows the researcher to focus more on selecting an appropriate model as opposed to ﬁnding manip-ulations to make the data ﬁt a restricted class of models.The ﬂexibility 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,modiﬁcations 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.Wolﬁnger.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基础知识 (1)

很全的sas基础知识SAS里面的PROC一览The ACECLUS Procedure :聚类的协方差矩阵近似估计（approximate covariance estimation for clustering）The ANOVA Procedure ：方差分析The BOXPLOT Procedure ：箱形图The CALIS Procedure ：结构方程模型The CANCORR Procedure ：典型相关分析The CANDISC Procedure ：主成分分析和典型相关分析The CATMOD Procedure ：类别分析The CLUSTER Procedure ：聚类分析，包括11种（average linkage, the centroid method, complete linkage, density linkage (including Wong’s hybrid and th-nearest-neighbor methods), maximum likelihood for mixtures of spherical multivariate normal distributions with equal variances but possibly unequal mixing proportions, the flexible-beta method, McQuitty’s similarity analysis, the median method, single linkage, two-stage density linkage, and Ward’s minimum-variance method，机器翻译为：平均联动，重心法，完全连锁，密度连接（包括Wong混合模型，最近邻的方法），最大的可能性，McQuitty的相似性分析，中位数法，单联动，两阶段密度联动，Ward最小方差法）。

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

广义线性混合效应模型试验设计若干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.。

glmm.hp 0.1-0 软件文档说明书

Package‘glmm.hp’May23,2023Type PackageTitle Hierarchical Partitioning of Marginal R2for GeneralizedMixed-Effect ModelsVersion0.1-0Date2023-5-23Depends R(>=3.4.0),MuMIn,ggplot2,veganImports lme4Maintainer Jiangshan Lai<************.cn>Description Conducts hierarchical partitioning to calculate individual contributions of each predic-tor(ﬁxed effects)towards marginal R2for generalized linear mixed-effect model(includ-ing lm,glm and glmm)based on output of r.squaredGLMM()in'MuMIn',applying the algo-rithm of Lai J.,Zou Y.,Zhang S.,Zhang X.,Mao L.(2022)glmm.hp:an R package for comput-ing individual effect of predictors in generalized linear mixed models.Journal of Plant Ecol-ogy,15(6)1302-1307<doi:10.1093/jpe/rtac096>.License GPLEncoding UTF-8URL https:///laijiangshan/glmm.hpRoxygenNote7.1.1NeedsCompilation noAuthor Jiangshan Lai[aut,cre](<https:///0000-0002-0279-8816>), Kim Nimon[aut]Repository CRANDate/Publication2023-05-2302:50:02UTCR topics documented:glmm.hp (2)plot.glmmhp (3)Index512glmm.hpglmm.hp Hierarchical Partitioning of Marginal R2for Generalized Mixed-Effect ModelsDescriptionHierarchical Partitioning of Marginal R2for Generalized Mixed-Effect ModelsUsageglmm.hp(mod,type="adjR2",commonality=FALSE)Argumentsmod Fitted lme4,nlme,glmmTMB,glm or lm model objects.type The type of R-square of lm,either"R2"or"adjR2",in which"R2"is unadjustedR-square and"adjR2"is adjusted R-square,the default is"adjR2".The adjustedR-square is calculated using Ezekiel’s formula(Ezekiel1930)for lm.commonality Logical;If TRUE,the result of commonality analysis(2^N-1fractions for Npredictors)is shown,the default is FALSE.DetailsThis function conducts hierarchical partitioning to calculate the individual contributions of each predictor towards total(marginal)R2for Generalized Linear Mixed-effect Model(including lm,glm and glmm).The marginal R2is the output of r.squaredGLMM in MuMIn package for glm and glmm.Valuer.squaredGLMM The R2for the full model.hierarchical.partitioningA matrix containing individual effects and percentage of individual effects to-wards total(marginal)R2for each predictor.Author(s)Jiangshan Lai<************.cn>Kim Nimon<*******************>References•Lai J.,Zou Y.,Zhang S.,Zhang X.,Mao L.(2022)glmm.hp:an R package for computing indi-vidual effect of predictors in generalized linear mixed models.Journal of Plant Ecology,15(6):1302-1307<DOI:10.1093/jpe/rtac096>•Lai J.,Zou Y.,Zhang J.,Peres-Neto P.(2022)Generalizing hierarchical and variation partition-ing in multiple regression and canonical analyses using the rdacca.hp R package.Methods in Ecology and Evolution,13(4):782-788<DOI:10.1111/2041-210X.13800>•Chevan,A.&Sutherland,M.(1991).Hierarchical partitioning.American Statistician,45, 90-96.doi:10.1080/00031305.1991.10475776•Nimon,K.,Oswald,F.L.&Roberts,J.K.(2013).Yhat:Interpreting regression effects.R package version2.0.0.•Nakagawa,S.,&Schielzeth,H.(2013).A general and simple method for obtaining R2from generalized linear mixed-effects models.Methods in Ecology and Evolution,4(2),133-142.•Nakagawa,S.,Johnson,P.C.,&Schielzeth,H.(2017).The coefﬁcient of determination R2 and intra-class correlation coefﬁcient from generalized linear mixed-effects models revisited and expanded.Journal of the Royal Society Interface,14(134),20170213.•Ezekiel,M.(1930)Methods of Correlational Analysis.Wiley,New York.Exampleslibrary(MuMIn)library(lme4)mod1<-lmer(Sepal.Length~Petal.Length+Petal.Width+(1|Species),data=iris)r.squaredGLMM(mod1)glmm.hp(mod1)plot(glmm.hp(mod1))mod2<-glm(Sepal.Length~Petal.Length+Petal.Width,data=iris)r.squaredGLMM(mod2)glmm.hp(mod2)plot(glmm.hp(mod2))mod3<-lm(Sepal.Length~Petal.Length+Petal.Width,data=iris)glmm.hp(mod3,type="R2")glmm.hp(mod3,commonality=TRUE)plot.glmmhp Plot for a glmm.hp objectDescriptionPlot for a glmm.hp objectUsage##S3method for class glmmhpplot(x,plot.perc=FALSE,color=NULL,n=1,...)Argumentsx A glmm.hp object.plot.perc Logical;if TRUE,the bar plot(based on ggplot2package)of the percentage to individual effects of variables or groups towards total explained variation,thedefault is FALSE to show plot with original individual effects.color Color of variables.n Integer;which marginal R2in output of r.squaredGLMM to plot....unusedValuea ggplot objectAuthor(s)Jiangshan Lai<************.cn>Exampleslibrary(MuMIn)library(lme4)mod1<-lmer(Sepal.Length~Petal.Length+Petal.Width+(1|Species),data=iris) plot(glmm.hp(mod1))mod3<-lm(Sepal.Length~Petal.Length+Petal.Width,data=iris)plot(glmm.hp(mod3,type="R2"))plot(glmm.hp(mod3,commonality=TRUE))Indexglmm.hp,2,3,4plot.glmmhp,35。

glmmTMB包用户指南说明书

Getting started with the glmmTMB packageBen BolkerOctober7,20231Introduction/quick startglmmTMB is an R package built on the Template Model Builder automatic differentiation engine,for fitting generalized linear mixed models and exten-sions.(Not-yet-implemented features are denoted like this)response distributions:Gaussian,binomial,beta-binomial,Poisson, negative binomial(NB1and NB2parameterizations),Conway-Maxwell-Poisson,generalized Poisson,Gamma,Beta,Tweedie;as well as zero-truncated Poisson,Conway-Maxwell-Poisson,generalized Poisson,and negative binomial;Student t.See?family glmmTMB for a current list including details of parameterizations.link functions:log,logit,probit,complementary log-log,inverse,iden-tityzero-inflation with fixed and random-effects components;hurdle models via truncated Poisson/NBsingle or multiple(nested or crossed)random effectsoffsetsfixed-effects models for dispersiondiagonal,compound-symmetric,or unstructured random effects variance-covariance matrices;first-order autoregressive(AR1)variance struc-tures1In order to use glmmTMB effectively you should already be reasonably fa-miliar with generalized linear mixed models(GLMMs),which in turn requires familiarity with(i)generalized linear models(e.g.the special cases of logis-tic,binomial,and Poisson regression)and(ii)‘modern’mixed models(those working via maximization of the marginal likelihood rather than by manipu-lating sums of squares).Bolker et al.(2009)and Bolker(2015)are reasonable starting points in this area(especially geared to biologists and less-technical readers),as are Zuur et al.(2009),Millar(2011),and Zuur et al.(2013).In order to fit a model in glmmTMB you need to:specify a model for the conditional effects,in the standard R(Wilkinson-Rogers)formula notation(see?formula or Section11.1of the Intro-duction to R.Formulae can also include offsets.specify a model for the random effects,in the notation that is commonto the nlme and lme4packages.Random effects are specified as x|g,where x is an effect and g is a grouping factor(which must be a factorvariable,or a nesting of/interaction among factor variables).For ex-ample,the formula would be1|block for a random-intercept model ortime|block for a model with random variation in slopes through timeacross groups specified by block.A model of nested random effects(block within site)could be1|site/block if block labels are reusedacross multiple sites,or(1|site)+(1|block)if the nesting structureis explicit in the data and each level of block only occurs within onesite.A model of crossed random effects(block and year)would be(1|block)+(1|year).choose the error distribution by specifying the family(family argu-ment).In general,you can specify the function(binomial,gaussian,poisson,Gamma from base R,or one of the options listed at family glmmTMB [nbinom2,beta family(),betabinomial,...])).choose the error distribution by specifying the family(family argu-ment).You can choose among–Distributions defined in base R and documented in?family,ex-cept the inverse Gaussian and quasi-families–Distributions defined in?glmmTMB::family glmmTMB2optionally specify a zero-inflation model(via the ziformula argument) with fixed and/or random effectsoptionally specify a dispersion model with fixed effectsThis document was generated using4.3.1and package versions:##bbmle glmmTMB## 1.0.25 1.1.8The current citation for glmmTMB is:Brooks ME,Kristensen K,van Benthem KJ,Magnusson A,BergCW,Nielsen A,Skaug HJ,Maechler M,Bolker BM(2017).“glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflatedGeneralized Linear Mixed Modeling.”The R Journal,9(2),378–400.doi:10.32614/RJ-2017-066.2Preliminaries:packages and dataLoad required packages:library("glmmTMB")library("bbmle")##for AICtablibrary("ggplot2")##cosmetictheme_set(theme_bw()+theme(panel.spacing=grid::unit(0,"lines")))The data,taken from Zuur et al.(2009)and ultimately from Roulin and Bersier(2007),quantify the number of negotiations among owlets(owl chicks)in different nests prior to the arrival of a provisioning parent as a function of food treatment(deprived or satiated),the sex of the parent,and arrival time.The total number of calls from the nest is recorded,along with the total brood size,which is used as an offset to allow the use of a Poisson response.Since the same nests are measured repeatedly,the nest is used as a random effect.The model can be expressed as a zero-inflated generalized linear mixed model(ZIGLMM).3Various small manipulations of the data set:(1)reorder nests by mean negotiations per chick,for plotting purposes;(2)add log brood size variable(for offset);(3)rename response variable and abbreviate one of the input variables.Owls<-transform(Owls,Nest=reorder(Nest,NegPerChick),NCalls=SiblingNegotiation,FT=FoodTreatment)(If you were really using this data set you should start with summary(Owls)to explore the data set.)We should explore the data before we start to build models,e.g.by plotting it in various ways,but this vignette is about glmmTMB,not aboutdata visualization...Now fit some models:The basic glmmTMB fit—a zero-inflated Poisson model with a single zero-inflation parameter applying to all observations(ziformula~1).(Excludingzero-inflation is glmmTMB’s default:to exclude it explicitly,use ziformula~0.)fit_zipoisson<-glmmTMB(NCalls~(FT+ArrivalTime)*SexParent+offset(log(BroodSize))+(1|Nest),data=Owls,ziformula=~1,family=poisson)summary(fit_zipoisson)We can also try a standard zero-inflated negative binomial model;the default is the“NB2”parameterization(variance=µ(1+µ/k):Hardinand Hilbe(2007)).To use families(Poisson,binomial,Gaussian)that are defined in R,you should specify them as in?glm(as a string referringto the family function,as the family function itself,or as the result ofa call to the family function:i.e.family="poisson",family=poisson, family=poisson(),and family=poisson(link="log")are all allowed andall equivalent(the log link is the default for the Poisson family).Someof the additional families that are not defined in base R(at this point4nbinom2and nbinom1)can be specified using the same format.Other-wise,for families that are implemented in glmmTMB but for which glmmTMB does not provide a function,you should specify the family argument as a list containing(at least)the(character)elements family and link,e.g. family=list(family="nbinom2",link="log").(In order to be able to re-trieve Pearson(variance-scaled)residuals from a fit,you also need to specify a variance component;see?family glmmTMB.)fit_zinbinom<-update(fit_zipoisson,family=nbinom2) Alternatively,we can use an“NB1”fit(variance=ϕµ).fit_zinbinom1<-update(fit_zipoisson,family=nbinom1) Relax the assumption that total number of calls is strictly proportional to brood size(ing log(brood size)as an offset):fit_zinbinom1_bs<-update(fit_zinbinom1,.~(FT+ArrivalTime)*SexParent+BroodSize+(1|Nest))Every change we have made so far improves the fit—changing distri-butions improves it enormously,while changing the role of brood size makes only a modest(-1AIC unit)difference:AICtab(fit_zipoisson,fit_zinbinom,fit_zinbinom1,fit_zinbinom1_bs)2.1Hurdle modelsIn contrast to zero-inflated models,hurdle models treat zero-count and non-zero outcomes as two completely separate categories,rather than treating the zero-count outcomes as a mixture of structural and sampling zeros.glmmTMB includes truncated Poisson and negative binomial familes and hence can fit hurdle models.5fit_hnbinom1<-update(fit_zinbinom1_bs,ziformula=~.,data=Owls,family=truncated_nbinom1)Then we can use AICtab to compare all the models.AICtab(fit_zipoisson,fit_zinbinom,fit_zinbinom1,fit_zinbinom1_bs,fit_hnbinom1)3Sample timingsTo get a rough idea of glmmTMB’s speed relative to lme4(the most commonly used mixed-model package for R),we try a few standard problems,enlarging the data sets by cloning the original data set(making multiple copies and sticking them together).Figure??shows the results of replicating the Contraception data set (1934observations,60levels in the random effects grouping level)from1 to40times.glmmADMB is sufficiently slow(≈1minute for a single copy of the data)that we didn’t try replicating very much.On average,glmmTMB is about2.3times faster than glmer for this problem.Figure1shows equivalent timings for the InstEval data set,although in this case since the original data set is large(73421observations)we subsample the data set rather than cloning it:in this case,the advantage is reversed and lmer is about5times faster.In general,we expect glmmTMB’s advantages over lme4to be(1)greater flexibility(zero-inflation etc.);(2)greater speed for GLMMs,especially those with large number of“top-level”parameters(fixed effects plus random-effects variance-covariance parameters).In contrast,lme4should be faster for LMMs(for maximum speed,you may want to check the MixedModels.jl package for Julia);lme4is more mature and at present has a wider variety of diagnostic checks and methods for using model results,including downstream packages.6Figure1:Timing for fitting subsets of the InstEval data set.7ReferencesBolker, B.M.(2015).Linear and generalized linear mixed models.In G.A.Fox,S.Negrete-Yankelevich,and V.J.Sosa(Eds.),Ecological Statis-tics:Contemporary theory and application,Chapter13.Oxford University Press.Bolker,B.M.,M.E.Brooks,C.J.Clark,S.W.Geange,J.R.Poulsen, M.H.H.Stevens,and J.S.White(2009).Generalized linear mixed mod-els:a practical guide for ecology and evolution.Trends in Ecology& Evolution24,127–135.Hardin,J.W.and J.Hilbe(2007,February).Generalized linear models and extensions.Stata Press.Millar,R.B.(2011,July).Maximum Likelihood Estimation and Inference: With Examples in R,SAS and ADMB.John Wiley&Sons.Roulin,A.and L.Bersier(2007).Nestling barn owls beg more intensely in the presence of their mother than in the presence of their father.Animal Behaviour74,1099–1106.Zuur,A.F.,J.M.Hilbe,and E.N.Leno(2013,May).A Beginner’s Guide to GLM and GLMM with R:A Frequentist and Bayesian Perspective for Ecologists.Highland Statistics Ltd.Zuur,A.F.,E.N.Ieno,N.J.Walker,A.A.Saveliev,and G.M.Smith(2009, March).Mixed Effects Models and Extensions in Ecology with R(1ed.). Springer.8。

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

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
glmm
Department of Biostatistics University of Copenhagen
Problem: The Likelihood
The likelihood of a GLMM model will involve an integral, which cannot in general be computed explicitly. Following Molenberghs and Verbeke, this has contribution has the form
glmm
Department of Biostatistics University of Copenhagen
The PQL approach
A Taylor expansion of Yij = h(xij β + zij bi ) + εij yields (in vector notation) ˆ ˆ b Vi−1 (Yi − µi ) + Xi β + Zi ˆi ≈ Xi β + Zi bi + ε∗ ˆ i which gives a straightforward updating scheme, considering the LHS as pseudo data. This is known as penalized quasi likelihood because it obtains from optimizing a quasi-likelihood (involving only 1st and 2nd derivatives) with a penalty term on the random effects.
glmm
Department of Biostatistics University of Copenhagen
R Functions for GLMM
lmer used to be able to handle all three approaches, but currently cannot do AGQ Model speciﬁcation is quite trivial: Just add a family argument Older PQL code is glmmPQL in the MASS package (based on lme)
Overview of GLMM
Generalized Linear Models Mixed effects extensions Approximations Examples
glmm
Department of Biostatistics University of Copenhagen
Generalized Linear Models
glmm
Department of Biostatistics University of Copenhagen
Example 2: Pressure sensitivity of teeth
(Misc. data manipulation skipped)
d$response <- c(.5,1,0)[d$hilo] d$type <- gl(2,240,labels=c("Implant","Control")) d$dif <- with(d,test-ref) library(lme4) fit <- lmer(response~dif*type + (dif-1 |ID) + (dif-1|ID:type), data=d,family=quasibinomial) ff <- plogis(fitted(fit)) plot(ff,d$response-ff) lines(smooth.spline(ff, d$response-ff)) plot(ff,(d$response-ff)^2) lines(smooth.spline(ff, (d$response-ff)^2))
glmm
Department of Biostatistics University of Copenhagen
Adaptive Gauss Quadrature
This approximates the integral by evaluating it at a set of quadrature point and obtaining an approximation of the form wi f (xi ) Optimal weights and xi have been tabulated for expansions in terms of Hermite polynomials. The “adaptive” refers to scaling the integrand using the Hessian at the optimum point, like in the Laplace method. Drawback: Complexity increases exponentially with dimension of random effects vector.
Statistical distribution (exponential) family Link function transforming mean to linear scale Deviance Examples; Binomial, Poisson, Gaussian (σ known — in principle) Canonical link functions: logit, log, identity Fit using glm in R
glmm
Department of Biostatistics University of Copenhagen
The Laplace Method
A trick from Bayes theory: Likelihoods often look like Gaussian densities, which have known integrals The log-likelihood is then close to parabolic in shape Approximate using value at peak and Hessian (comes automatically from optimizers) Apply to the joint likelihood of (yi , bi )
glmm
Department of Biostatistics University of Copenhagen
Generalized Linear Mixed Models
Basic idea: Just add random effects on the linear scale Prototypically: Logistic regression with random effects log(p/(1 − p)) = Xi β + Zi bi (Notice that there is no error term because we are modeling probabilities) Typically, a normal distribution is assumed for the random effects. Notice that the interpretation can be quite different from a GEE analysis.
Generalized Linear Mixed Models
Peter Dalgaard
Department of Biostatistics University of CoR, January 2006
glmm
Department of Biostatistics University of Copenhagen