Multilevel Logistic Regression: Respiritory



Multilevel Logistic Regression: RespiritoryCarolyn J. Anderson3/20/2020Table of ContentsTOC \o "1-3" \h \z \uSetup PAGEREF _Toc35598541 \h 1Logistic regression PAGEREF _Toc35598542 \h 2Generalized Estimating Equations PAGEREF _Toc35598543 \h 3Random effects logistic regression via glmer PAGEREF _Toc35598544 \h 4LaPlace PAGEREF _Toc35598545 \h 5MLE: quass quadrature PAGEREF _Toc35598546 \h 5The R in this document reproduces the results in the lecture on multilevel logistic regression. The data are from Skrondal & Rabe-Hesketh (2004) which were also analyzed by Zeger & Karim (1991), Diggle et al.?(2002), but originally from Sommer et al.?(1983).SetupThe package that we’ll uselibrary(lme4)library(lmerTest)library(gee) # for generalized estimating equationslibrary(MASS) # For confidence limits for parametersAnd read in the dataresp <- read.table(file="D:/Dropbox/edps587/lectures/10 logreg/respitory.txt", header=TRUE)head(resp)## id resp cons age xero cosine sine female height stunted time age1## 1 121013 0 1 31 0 -1 0 0 -3 0 1 31## 2 121013 0 1 34 0 0 -1 0 -3 0 2 31## 3 121013 0 1 37 0 1 0 0 -2 0 3 31## 4 121013 0 1 40 0 0 1 0 -2 0 4 31## 5 121013 1 1 43 0 -1 0 0 -2 0 5 31## 6 121013 0 1 46 0 0 -1 0 -3 0 6 31## season time2## 1 2 1## 2 3 4## 3 4 9## 4 1 16## 5 2 25## 6 3 36Logistic regressionThis is the example shown on page 26, except that in notes I give Wald# glm = Generalized Linear modelmodel1 <- glm(resp ~ 1 + age, data=resp, family=binomial)summary(model1)## ## Call:## glm(formula = resp ~ 1 + age, family = binomial, data = resp)## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -0.6204 -0.4988 -0.3940 -0.3135 2.6238 ## ## Coefficients:## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.343573 0.105297 -22.257 < 2e-16 ***## age -0.024795 0.005558 -4.461 8.15e-06 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## (Dispersion parameter for binomial family taken to be 1)## ## Null deviance: 721.45 on 1199 degrees of freedom## Residual deviance: 699.69 on 1198 degrees of freedom## AIC: 703.69## ## Number of Fisher Scoring iterations: 5Note that z^2 above equals Wald statistics in the lecture notes. These are compared to N(0,1) to get the p-value.The following puts 95% confidence levels on parameters (i.e., the coefficients or beta’s)round(confint(model1, level = 0.95),digits=2)## Waiting for profiling to be done...## 2.5 % 97.5 %## (Intercept) -2.56 -2.14## age -0.04 -0.01For interpetation of the resultsodds <- exp(coefficients(model1))round(odds,digits=2)## (Intercept) age ## 0.10 0.98odds2 <- 1/oddsround(odds2,digits=2) ## (Intercept) age ## 10.42 1.03Out of curiosity and as a check on whether variance function of mean (i.e., binomial distribution OK).model1b <- glm(resp ~ 1 + age, data=resp, family=quasibinomial)summary(model1b)## ## Call:## glm(formula = resp ~ 1 + age, family = quasibinomial, data = resp)## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -0.6204 -0.4988 -0.3940 -0.3135 2.6238 ## ## Coefficients:## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -2.343573 0.103890 -22.558 < 2e-16 ***## age -0.024795 0.005484 -4.522 6.75e-06 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## (Dispersion parameter for quasibinomial family taken to be 0.9734542)## ## Null deviance: 721.45 on 1199 degrees of freedom## Residual deviance: 699.69 on 1198 degrees of freedom## AIC: NA## ## Number of Fisher Scoring iterations: 5Since the dispersion paramters is is close to 1, binomial is probably just fine for these data.Generalized Estimating EquationsThis yields a population averge modelmodel.gee <- gee(resp ~1 + age, id, data=resp,family=binomial,corstr="exchangeable")## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27## running glm to get initial regression estimate## (Intercept) age ## -2.34357294 -0.02479512summary(model.gee)## ## GEE: GENERALIZED LINEAR MODELS FOR DEPENDENT DATA## gee S-function, version 4.13 modified 98/01/27 (1998) ## ## Model:## Link: Logit ## Variance to Mean Relation: Binomial ## Correlation Structure: Exchangeable ## ## Call:## gee(formula = resp ~ 1 + age, id = id, data = resp, family = binomial, ## corstr = "exchangeable")## ## Summary of Residuals:## Min 1Q Median 3Q Max ## -0.17380518 -0.11711642 -0.07547880 -0.04895995 0.96704630 ## ## ## Coefficients:## Estimate Naive S.E. Naive z Robust S.E. Robust z## (Intercept) -2.33553443 0.112748751 -20.714504 0.113371780 -20.600668## age -0.02426996 0.005886584 -4.122928 0.005143234 -4.718814## ## Estimated Scale Parameter: 0.9661839## Number of Iterations: 2## ## Working Correlation## [,1] [,2] [,3] [,4] [,5] [,6]## [1,] 1.00000000 0.04999101 0.04999101 0.04999101 0.04999101 0.04999101## [2,] 0.04999101 1.00000000 0.04999101 0.04999101 0.04999101 0.04999101## [3,] 0.04999101 0.04999101 1.00000000 0.04999101 0.04999101 0.04999101## [4,] 0.04999101 0.04999101 0.04999101 1.00000000 0.04999101 0.04999101## [5,] 0.04999101 0.04999101 0.04999101 0.04999101 1.00000000 0.04999101## [6,] 0.04999101 0.04999101 0.04999101 0.04999101 0.04999101 1.00000000Random effects logistic regression via glmerThe command glmer is an option in the lme4 package.LaPlaceNow for the multilevel model or random effects model using the default estimation method (i.e., LaPlace)mod1.laplace <- glmer(resp ~ 1 + age + (1 | id), data=resp, family=binomial)summary(mod1.laplace)## Generalized linear mixed model fit by maximum likelihood (Laplace## Approximation) [glmerMod]## Family: binomial ( logit )## Formula: resp ~ 1 + age + (1 | id)## Data: resp## ## AIC BIC logLik deviance df.resid ## 696.7 711.9 -345.3 690.7 1197 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -0.8165 -0.2981 -0.2309 -0.1770 4.8868 ## ## Random effects:## Groups Name Variance Std.Dev.## id (Intercept) 0.8949 0.946 ## Number of obs: 1200, groups: id, 275## ## Fixed effects:## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.676724 0.184338 -14.521 < 2e-16 ***## age -0.026678 0.006735 -3.961 7.47e-05 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Correlation of Fixed Effects:## (Intr)## age 0.163MLE: quass quadratureIf you want MLE, glmer uses Guass quadrature — MLE gold standard. glmer will use adaptive guass quadrature and the number of quadrature points is indicated by nAGQ. These results are reported ~ page 50,mod1.quad <- glmer(resp ~ 1 + age + (1 | id),data=resp, family=binomial, nAGQ=10 )summary(mod1.quad)## Generalized linear mixed model fit by maximum likelihood (Adaptive## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]## Family: binomial ( logit )## Formula: resp ~ 1 + age + (1 | id)## Data: resp## ## AIC BIC logLik deviance df.resid ## 698.6 713.8 -346.3 692.6 1197 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -0.7688 -0.3095 -0.2392 -0.1847 5.0157 ## ## Random effects:## Groups Name Variance Std.Dev.## id (Intercept) 0.7273 0.8528 ## Number of obs: 1200, groups: id, 275## ## Fixed effects:## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.612985 0.172299 -15.165 < 2e-16 ***## age -0.026760 0.006592 -4.059 4.92e-05 ***## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## ## Correlation of Fixed Effects:## (Intr)## age 0.195Computing the confidence intervals take a bit of time. The following gives 95% profile confidence intervals for all parametersround(confint(mod1.quad),digits=4)## Computing profile confidence intervals ...## 2.5 % 97.5 %## .sig01 0.3993 1.2727## (Intercept) -2.9887 -2.3059## age -0.0403 -0.0142For interpreations, we use the oddsodds <- exp(fixef(mod1.quad))round(odds,digits=2)## (Intercept) age ## 0.07 0.97 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download