http://blog.gmane.org/gmane.comp.lang.r.lme4.devel hourly 1 1901-01-01T00:00+00:00 http://gmane.org/img/gmane-25t.png http://gmane.org http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8264 <pre>Hi all, I am trying to get predictions for individual level observations for an lme object. Yet, since I am get different types of errors for different trials, it seems to me I am missing something. My model is the following: model &lt;- lme(log(child_mortality) ~ as.factor(cluster)*time + my.new.time.one.transition.low.and.middle + ttd + maternal_educ+ log(IHME_id_gdppc) + hiv_prev-1, merged0,na.action=na.omit,method="ML",weights=varPower(form = ~ time), random= ~ time| country.x, correlation=corAR1(form = ~ time), control=lmeControl(msMaxIter = 200, msVerbose = TRUE)) It runs fine and the results make sense. Now my attempts to get predictions. For example: data.frame(time=c(10,10,10,10),country.x=c("Poland","Brazil","Argentina","France"), + my.new.time.one.transition.low.and.middle=c(1,1,1,0),ttd=c(0,0,0,0), + maternal_educ=c(10,10,10,10),IHME_id_gdppc=c(log(5000),log(8000),log(8000),log(15000)), + hiv_prev=c(.005,.005,.005,.005), + cluster=c("One Transition, Middle Income","One Transition, Middle Income", + "One Transition, Middle Income","Democracy, High Income")) Error in X %*% fixef(object) : non-conformable arguments If I exclude, say, France, and only include countries in which cluster="One Transition, Middle Income" in my toy example, then I get a different error predict(model,test.pred,level=0) Error in `contrasts&lt;-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : contrasts can be applied only to factors with 2 or more levels My object is to calculate counterfactual scenarios (predictions and intervals for them) for different subsets of my data. Could you please help? Many thanks, Antonio Pedro. [[alternative HTML version deleted]] </pre> Antonio P. Ramos 2012-05-25T23:09:52 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8258 <pre>Dear list, The simplest random intercept model that I try to fit with the newest version of lme4 (0.999902344-0) throws an "Object 'multResp' not found" error on my Mac (Mac OSX 10.7.3, German R 2.15.0). Here is a reproducible example: install.packages(c("minqa", "Rcpp")) install.packages("lme4", repos="http://lme4.r-forge.r-project.org/repos") library(lme4) my_data &lt;- read.csv("http://dl.dropbox.com/u/5384027/my_data.csv") my_data$group &lt;- as.factor(my_data$group) # the data set consists of 97 four-person teams # with two variables x and y observed on the individual # level in each team head(my_data) # X group y x # 1 1 1 2.914286 -8.170984 # 2 2 1 2.746269 -13.504318 # 3 3 1 3.171429 -1.504318 # 4 4 1 2.978723 -4.837651 # 5 5 2 2.928571 -8.170984 # 6 6 2 2.987013 8.495682 mlmodel1_ri &lt;- lmer(y ~ x + (1 | group), data = my_data) # Fehler in lmer(y ~ x + (1 | group), data = my_data) : # Objekt 'multResp' nicht gefunden Does anyone know how to fix this? Greetings, Bertolt </pre> Bertolt Meyer 2012-05-24T15:30:07 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8254 <pre> Dear R-sig-mixed-models subscribers, I am analizing survival (ALIVE) of one year-old seedlings of over a hundred species as a function of the total number of seedlings (TOTABN) and the proportion of conspecifics seedlings (CONSPp) in their neighborhood (sampling quadrat). The data comes from monitoring seedling dynamics trough time (YEAR) in a series of forest plots (PLOT). A crucial part of the analysis is to understand the variation in responses between species (SPECIES). For now, I have one continuous species-level variable, abundance across the landscape (LANDSPABN), that I would like to use to explain the variation between species. After reading a number of printed and online references on advanced statistical models it seems to me that my problem could be properly analyzed using hierarchical bayesian modelling but I am just beggining to understand lmer (and R) and I was wondering if using such tools I can also perform a meaningfull analysis.My first approach in lmer was to specify a model with varying intercept and slope for the effect of conspecifics, i.e. (1+CONSPp|SPECIES) and to extract from the results the conditional modes (using terminology in Prof. Bates draft book) of the random effects (and their variances) to model them separately as a function other continuous predictors (using simple weighted linear regressions). The full model of this approach is: M1&lt;-lmer(ALIVE~CONSPp+TOTABN+(1+CONSPp|SPECIES)+(1|PLOT)+(1|YEAR), data=oneyrseedl,family=binomial). After further reading I found an example from Gelman and Hill 2007 in which he includes "group level predictors". Specifically in their example, a predictor X1 appears both as fixed (y~ X1) and random effect (1+ X1 | grouping factor) with the grouping factor predictor (X2) modelled as an interaction with X1. So their full model is y ~ X1 + X2 + X1:X2+ (1 + X1 | grouping factor). I thought I may be able to specify a similar model to directly incorporate my continuous species-level variable (species abundance across the landscape) to explain the variation in the effect of conspecifics on seedling survival. One crucial difference between my case and Gelman and Hill´s case is that my continuous species-level variable has only one value per species (my grouping factor) while their X2 has a series of values per grouping factor. I thus think that if include the continuous species-level variable in the model I do not need to include a random slope for the effects of conspecifics (with SPECIES as a grouping factor) as I feel it would control for all the variation that could be more meaningfully explained by my species-level variable (or any other species-level variable I may come up later). Thus, the model I am thinking about would be: M2&lt;-lmer(ALIVE~CONSPp+TOTABN+ LANDSPABN+ CONSPp:LANDSPABN + (1|PLOT)+(1|YEAR), data=oneyrseedl,family=binomial). I am doubting if a random intercept by SPECIES is still needed (i.e. (1|fSPECIES).I hope I have exposed my problem clearly. I would highly appreciate any insight from the mixed model experts.Kind regards,Edwin [[alternative HTML version deleted]] </pre> Edwin Lebrija Trejos 2012-05-23T03:35:21 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8247 <pre>I know how to obtain the fixed effects from a mer object: or, But how can I obtain the t-statistic values for the fixed effects? Thanks, Gang </pre> Gang Chen 2012-05-22T17:56:16 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8237 <pre>Hello all, The following (i.e. below the hash marks) is a reposting of something originally posted on the R forum. I was told by a helpful person there to redirect my question to this forum. My question concerns what the correct expression is to specify two random factors for the lme function. So far I've tried: random = ~ 1|C + D random = ~ 1 + (1|C) + (1|D) random = ~ 1 + 1|C + 1|D but all of these resulted in errors (mostly of the type " Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups"). To be more specific, "C" refers to a subject factor, while "D" refers to a repetition factor (i.e. multiple evaluations of the same stimulus by the same subject). In this case, would the correct way to code this be: random = ~ 1|C/D ? Additionally, I would be grateful if someone could tell me whether I have to code the interactions in the "random" argument in the lme function, and if so how. Thank you very much in advance for the help. Best regards, Daisuke ############################################################################ Hello, I have a question regarding the syntax of the lme function in the nlme package. What I'm trying to do is to calculate an estimate of R^2 based on the likelihood ratio test. For this calculation, I need to determine the maximum log-likelihood of the intercept-only model and the model of interest (with the desired factors and interactions). My model has four independent factors (i.e. A, B, C and D). A and B are fixed factors, whilst C and D are random factors. Presently, I'm trying to code the lme function for the full ANOVA model, but am unsure how to code both random factors. Additionally, I'm unsure whether I coded the interactions correctly (i.e. I'm unsure whether the interactions which contain random factors also need to be specified in the "random" argument). I've skimmed through the Pinheiro and Bates book, "Mixed-Effects Models in S and S-PLUS", but could not find an answer. My R code for the lme function which includes all the main factors and interactions (2- and 3-way) so far is: lme(Y ~A + B + C + D + A : B + A : C + A : D + B : C + B : D + C : D + A : B : C + A : B : D + A : C : D + B : C : D, data = myData, random = ~ 1|C, method = "ML") Can someone help me with how to code the random factors and interactions which contain these random factors please? I apologise in advance, that I am self-taught in statistics, and just started using R. Hence I hope my question made sense. Thank you very much. Daisuke Koya Confused PhD student [[alternative HTML version deleted]] </pre> Daisuke Koya 2012-05-22T13:16:25 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8234 <pre>Hi I have a repeated measures dataset with unequal intervals (spacing) - the intervals are Day 6, Day 13 and Day 15. Anyone knows how to handle this type of data? Thanks Ahmad [[alternative HTML version deleted]] </pre> Ahmad Rabiee 2012-05-22T04:51:12 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8229 <pre>Sorry, plots are two large. Please read this message. Dear list, I'm running a lmer (package lme4) with a binomial error distribution and "bird" as the random effect (160 observations of 25 birds). The response variable "ars1" is coded as 0, 1. The fixed effect "sit" is numerical and "dive" is categorical (0, 1). What puzzles me a little is that the a variance and sd of the random effect is ZERO. Same question has been posted before and Douglas Bates answer was: "No, an estimate of zero is not suspicious. It is simply an indication that the variability between individuals is not significantly larger than what one would expect from the random variability in the response." While another answer suggested that the model was "wrong": "A zero estimate of a variance possibly indicates the model is wrong." This wrong model seemed to be related to a negative covariation of one of the fixed effects ? My simplified model is: mod6 &lt;- lmer(ars1 ~ sit + dive + (1|bird), data=dat, family=binomial) Generalized linear mixed model fit by the Laplace approximation Formula: ars1 ~ sit + dive + (1 | bird) Data: dat AIC BIC logLik deviance 159.4 171.7 -75.71 151.4 Random effects: Groups Name Variance Std.Dev. bird (Intercept) 0 0 Number of obs: 160, groups: bird, 25 Fixed effects: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.3615 0.4037 -0.895 0.37059 sit 0.7492 0.1448 5.175 2.28e-07 *** dive 1.3076 0.4374 2.990 0.00279 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) sit sit 0.422 dive -0.756 0.005 Based on the summary output (zero variance and sd) I'm inclined to believe that in fact my random effect bird does not account for any of the variance in the model. I.e., that there is no significant variability between birds that I should account for. Further, if I plot *dotplot(ranef(mod6, postVar=TRUE))* to represent the 95% prediction intervals for each of the birds, these are all plotted on zero. The same for: *qqnorm(unlist(ranef(mod6)), main="normal qq-plot, random effects") qqline(unlist(ranef(mod6))) # qq of random effects* All dots are on zero! *QUESTION: Could I be overlooking something or is it justified to run a glm without the random effect bird instead of a lmer?* Thank you! Best regards, Julia dotplot(ranef(mod6, postVar=TRUE)) </pre> Julia Sommerfeld 2012-05-22T01:27:12 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8227 <pre>Dear list, I'm running a lmer (package lme4) with a binomial error distribution and "bird" as the random effect (160 observations of 25 birds). The response variable "ars1" is coded as 0, 1. The fixed effect "sit" is numerical and "dive" is categorical (0, 1). What puzzles me a little is that the a variance and sd of the random effect is ZERO. Same question has been posted before and Douglas Bates answer was: "No, an estimate of zero is not suspicious. It is simply an indication that the variability between individuals is not significantly larger than what one would expect from the random variability in the response." While another answer suggested that the model was "wrong": "A zero estimate of a variance possibly indicates the model is wrong." This wrong model seemed to be related to a negative covariation of one of the fixed effects ? My simplified model is: mod6 &lt;- lmer(ars1 ~ sit + dive + (1|bird), data=dat, family=binomial) Generalized linear mixed model fit by the Laplace approximation Formula: ars1 ~ sit + dive + (1 | bird) Data: dat AIC BIC logLik deviance 159.4 171.7 -75.71 151.4 Random effects: Groups Name Variance Std.Dev. bird (Intercept) 0 0 Number of obs: 160, groups: bird, 25 Fixed effects: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.3615 0.4037 -0.895 0.37059 sit 0.7492 0.1448 5.175 2.28e-07 *** dive 1.3076 0.4374 2.990 0.00279 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) sit sit 0.422 dive -0.756 0.005 Based on the summary output (zero variance and sd) and the two plots below, I'm inclined to believe that in fact my random effect bird does not account for any of the variance in the model. I.e., that there is no significant variability between birds that I should account for. *QUESTION: Could I be overlooking something or is it justified to run a glm without the random effect bird instead of a lmer?* Thank you! Best regards, Julia dotplot(ranef(mod6, postVar=TRUE)) qqnorm(unlist(ranef(mod6)), main="normal qq-plot, random effects") qqline(unlist(ranef(mod6))) # qq of random effects [[alternative HTML version deleted]] </pre> Julia Sommerfeld 2012-05-22T01:07:22 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8221 <pre>Greetings R users, I am trying to match some SAS output of a mixed model. After running PROC MIXED with the covariance structure as AR1 the output below is given. Now when I try to replicate this with lme, I get the correct degrees of freedom and somewhat close values. If I try with gls, I get the exact F value but the incorrect denominator degrees of freedom. Is there some syntax or parameter I can adjust to get lme to have the same F values as gls? That way the correct degrees of freedom would be applied and I would get the replicated model. The data is also included below. Thanks to any who can help with this issue. Regards, Charles Row ID Group Died Event_name var 1 12510 3 YES B -1.05257 2 12510 3 YES S45 -1.00000 3 12510 3 YES FR2 -1.14630 4 12510 3 YES FR8 -1.08831 5 12510 3 YES FR20 -1.03339 6 21510 3 NO B -0.87290 7 21510 3 NO S45 -1.22185 8 21510 3 NO FR2 -1.01592 9 21510 3 NO FR8 -1.06550 10 21510 3 NO PR48 -0.65326 11 22210 3 YES B -1.37059 12 30810 3 NO B -0.92959 13 30810 3 NO S45 -1.26680 14 30810 3 NO FR2 -0.99396 15 30810 3 NO FR8 -1.04769 16 30810 3 NO FR20 -1.04769 17 31510 4 YES B -1.50585 18 31510 4 YES S45 -1.35262 19 31510 4 YES FR2 -1.26520 20 31510 4 YES FR8 -1.11407 21 31510 4 YES FR20 -0.97757 22 40510 2 YES B -1.04721 23 40510 2 YES S45 -1.13371 24 51010 3 NO B -1.21467 25 51010 3 NO S45 -1.18177 26 51010 3 NO FR2 -1.23582 27 51010 3 NO FR8 -1.35262 28 51010 3 NO FR20 -0.69897 29 51010 3 NO PR48 -1.14086 30 51710 4 NO B -1.13371 31 51710 4 NO S45 -1.24336 32 51710 4 NO FR2 -0.85201 33 51710 4 NO FR8 -0.72700 34 51710 4 NO FR20 -0.71220 35 51710 4 NO PR48 -1.10735 36 51809 1 NO B -0.70752 37 51809 1 NO FR2 -0.58120 38 51809 1 NO FR8 -0.71332 39 51809 1 NO FR20 -0.37212 40 61410 4 NO B -1.25493 41 61410 4 NO S45 -1.24949 42 61410 4 NO FR2 -1.18310 43 61410 4 NO FR8 -0.87517 44 61410 4 NO FR20 -0.34679 45 61410 4 NO PR48 -0.76472 46 62110 3 NO B -1.34199 47 62110 3 NO S45 -1.23657 48 62110 3 NO FR2 -1.20412 49 62110 3 NO FR8 -0.72011 50 62110 3 NO FR20 -0.72308 51 62110 3 NO PR48 -0.89688 52 62209 3 NO B -0.74041 53 62209 3 NO S45 -0.99183 54 62209 3 NO FR2 -1.06854 55 62209 3 NO FR8 -0.94386 56 62209 3 NO FR20 -0.38817 57 71210 4 YES B -1.23657 58 71210 4 YES S45 -0.94808 59 71309 4 YES B -1.14630 60 71309 4 YES S45 -1.10568 61 71309 4 YES FR2 -0.94462 62 71309 4 YES FR8 -0.89211 63 71309 4 YES FR20 -0.98970 64 71309 4 YES PR48 -0.83387 65 81009 4 YES B -1.15989 66 81009 4 YES S45 -0.93779 67 81009 4 YES FR2 -0.86710 68 81009 4 YES FR8 -0.69379 69 81610 1 NO B -1.04769 70 81610 1 NO S45 -1.02780 71 81610 1 NO FR2 -0.96098 72 81610 1 NO FR8 -0.91829 73 81610 1 NO FR20 -0.68550 74 81610 1 NO PR48 -0.66374 75 81709 4 YES B -0.65521 76 82410 3 NO B -1.06702 77 82410 3 NO S45 -1.30103 78 82410 3 NO FR2 -0.96617 79 82410 3 NO FR8 -0.84497 80 82410 3 NO FR20 -0.47716 81 82410 3 NO PR48 -0.68131 82 91310 4 NO B -1.30803 83 91310 4 NO S45 -1.07212 84 91310 4 NO FR2 -0.94615 85 91310 4 NO FR8 -1.29414 86 91310 4 NO FR20 -0.81304 87 91310 4 NO PR48 -0.86519 88 91409 4 NO B -1.12321 89 91409 4 NO S45 -1.08831 90 91409 4 NO FR2 -1.11862 91 91409 4 NO FR8 -0.87976 92 91409 4 NO FR20 -0.87517 93 91409 4 NO PR48 -0.68867 94 92109 3 NO B -1.15120 95 92109 3 NO S45 -0.99525 96 92109 3 NO FR2 -0.95429 97 92109 3 NO FR8 -0.54914 98 92109 3 NO FR20 -0.67264 99 92109 3 NO PR48 -0.74256 100 92110 4 YES B -1.05012 101 92110 4 YES S45 -1.14086 102 92110 4 YES FR2 -0.91901 103 92110 4 YES FR8 -0.80521 104 92110 4 YES FR20 -0.88807 105 92710 4 NO B -1.20412 106 92710 4 NO S45 -1.21325 107 92710 4 NO FR2 -1.19518 108 92710 4 NO FR8 -0.92959 109 92710 4 NO FR20 -0.69315 110 92710 4 NO PR48 -0.73779 111 100410 1 NO B -0.92300 112 100410 1 NO S45 -1.07314 113 100410 1 NO FR2 -0.80162 114 100410 1 NO FR8 -0.63283 115 100410 1 NO FR20 -0.73119 116 100410 1 NO PR48 0.10102 117 101209 3 NO B -1.24718 118 101209 3 NO S45 -1.00524 119 101209 3 NO FR2 -1.01592 120 101209 3 NO FR8 -0.85201 121 101209 3 NO FR20 -0.51670 122 101209 3 NO PR48 -0.72239 123 101909 3 NO B -1.13549 124 101909 3 NO S45 -0.92665 125 101909 3 NO FR2 -0.98548 126 101909 3 NO FR8 -1.03152 127 101909 3 NO FR20 -0.71715 128 101909 3 NO PR48 -1.07935 129 102510 2 YES B -1.30103 130 102510 2 YES S45 -0.83387 131 110810 3 YES B -1.21968 132 110810 3 YES S45 -0.68048 133 110810 3 YES FR2 -1.09474 134 110909 2 YES B -1.09259 135 110909 2 YES S45 -0.91293 136 111510 3 NO B -1.39794 137 111510 3 NO S45 -1.42366 138 111510 3 NO FR2 -1.22621 139 111510 3 NO FR8 -0.98885 140 111510 3 NO FR20 -0.66055 141 111510 3 NO PR48 -0.84497 142 111609 1 NO B -0.79942 143 111609 1 NO S45 -0.87517 144 111609 1 NO FR2 -0.74041 145 111609 1 NO FR8 -0.64092 146 111609 1 NO FR20 -0.53835 147 111609 1 NO PR48 -0.68719 148 120610 4 NO B -1.05110 149 120610 4 NO S45 -1.45100 150 120610 4 NO FR2 -1.21183 151 120610 4 NO FR8 -1.04144 152 120610 4 NO FR20 -0.61101 153 120610 4 NO PR48 -0.81702 154 120709 1 NO B -0.84497 155 120709 1 NO S45 -0.99439 156 120709 1 NO FR2 -1.13430 157 120709 1 NO FR20 -0.57187 158 120709 1 NO PR48 -0.80632 159 121310 3 YES B -1.16115 160 121310 3 YES S45 -1.18310 dat=read.table("C:/…/subset.csv",sep=",",header=TRUE, na.strings=".") attach(dat) dat34=dat[Pig_group %in% c("3", "4"),] attach(dat34) dat34=within(dat34, { group=factor(group) Event_name=factor(Event_name) Died=factor(Died) ID=factor(ID) }) attach(dat34) contrasts(dat34$Event_name)=contr.sum(n=6) contrasts(dat34$group)=contr.sum(n=2) contrasts(dat34$Died)=contr.sum(n=2) # What the results should be from SAS AR1 var #Type 3 Tests of Fixed Effects #Effect NumDF DenDF F Value Pr &gt; F #Event_name 5 91 7.43 &lt;.0001 #Died 1 24 0.08 0.7756 #Event_name*Died 5 91 3.09 0.0128 fit.13=lme(var~Event_name*Died, data=liver34, random=~1|ID, corr=corAR1()) anova(fit.13, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 91 1342.7364 &lt;.0001 #Event_name 5 91 8.7143 &lt;.0001 #Died 1 24 0.0527 0.8204 #Event_name:Died 5 91 3.0909 0.0127 ###################################################### ###THIS DOES IT but not lme function (need dfs) fit.16=gls(var~Event_name*Died, data=liver34, corr=corAR1(, ~1|ID)) anova(fit.16, type="marginal", adjustSigma=F) #Denom. DF: 115 # numDF F-value p-value #(Intercept) 1 1496.9032 &lt;.0001 #Event_name 5 7.4304 &lt;.0001 #Died 1 0.0831 0.7736 #Event_name:Died 5 3.0871 0.0119 #Give exact replicated answers 1-pf(7.4304, 5, 91) 1-pf(.0831, 1, 24) 1-pf(3.0871, 5, 91) ####################################################### [[alternative HTML version deleted]] </pre> Charles Determan Jr 2012-05-21T12:45:29 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8214 <pre>Hello, I would like to specify the following model (in lmer() syntax) in MCMCglmm, but seem to be unable to do so. I receive the warning message "some fixed effects are not estimable and have been removed", although I do not believe this should be the case. Here's the lmer() model: fit = lmer(y ~ -1 + cyear:stratum + site + (-1+cyear|site),data=indat) And here's an attempt at fitting this in MCMCglmm(): fitmcmc = MCMCglmm(y ~ -1 + cyear:stratum + site,random=~idh(-1+cyear):site,data=indat) The data structure is relatively simple: stratum and site are categorical (technically site is nested within stratum, so I'd give bonus points if someone can help me specify separate variance components for the site-specific slopes within each stratum), and "cyear" is the continuous predictor for which I would like stratum-specific slopes with random slopes for each site (again, within stratum, if possible). Thanks, Chris Gast ----------------------------- Chris Gast cmgast-Re5JQEeQqe8AvxtiuMwx3w&lt; at &gt;public.gmane.org [[alternative HTML version deleted]] </pre> Chris Gast 2012-05-20T21:12:27 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8213 <pre>Dear list, I am making use of the "clmm" function from the "ordinal" package, for fitting mixed models with ordinal outcomes. When fitting quite a few models, I get the error: Error in if (maxGrad &lt; rho$ctrl$gradTol) { : missing value where TRUE/FALSE needed In addition: Warning message: In sqrt(phi2) : NaNs produced Can anyone help me to understand why this may be happening, and/or whether there's anything obvious I should be doing different? With some seemingly small/arbitrary specification changes--or when taking (some) random subsets of my dataset--I don't get this message. But I'd like to find a way of reducing the frequency with which I do get it. I don't think there's anything unusual about my call: mod1 &lt;- clmm(as.ordered(y) ~ x + (1 | fac1) + (1 | fac2), data=dat, link="logit", threshold="symmetric") I'm running: R version 2.15.0 (2012-03-30) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) Many thanks to anyone who can assist, Malcolm </pre> Malcolm Fairbrother 2012-05-20T15:49:25 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8212 <pre>Dear Mailinglist, I would be very glad to get some assistance with interpreting a tricky model output in LME4. The model is this: Y ~ -1 + Y_index + time + Y_index*x1 + Y_index*x2 + ( - 1 + Y_index | subject) What I am doing with this is modeling 9 items of a questionnaire as a multivariate response variable Y. Y_index is a categorical variable defining the number of item of the questionnaire (1 through 9), and I am checking with interaction effects on x1 and x2 whether these covariates have differential effects on my multivariate response. It is a longitudinal design with 5 measurement points, and I expect that x1 only affects some of my 9 items (the same for x2). The dataset is in long-long format (Y_index * time), so each subject has 9*5 lines. (I found the suggestion for that kind of analysis in Hox, 2010). The for me relevant part of the output looks like this: Fixed effects: Estimate Std. Error t value Y_index1 0.3161592 0.0780922 4.049 Y_index2 0.4685218 0.0775360 6.043 Y_index3 0.9531528 0.0969119 9.835 Y_index4 0.2366093 0.0898923 2.632 Y_index5 0.3055025 0.0955639 3.197 Y_index6 0.2581729 0.0819606 3.150 Y_index7 0.4556287 0.0817002 5.577 Y_index8 0.6027990 0.0691566 8.716 Y_index9 0.8697155 0.0620898 14.007 time 0.5726978 0.0374384 15.297 x1 0.0196260 0.0020225 9.704 x2 -0.0415874 0.0350631 -1.186 Y_index2:x1 0.0023080 0.0018770 1.230 Y_index3:x1 -0.0019870 0.0027166 -0.731 Y_index4:x1 -0.0006285 0.0023784 -0.264 Y_index5:x1 0.0033737 0.0026178 1.289 Y_index6:x1 0.0067164 0.0020428 3.288 Y_index7:x1 -0.0016510 0.0021435 -0.770 Y_index8:x1 -0.0080817 0.0020414 -3.959 Y_index9:x1 -0.0140743 0.0021874 -6.434 Y_index2:x2 0.0358944 0.0325015 1.104 Y_index3:x2 -0.0675878 0.0470604 -1.436 Y_index4:x2 0.0037518 0.0411980 0.091 Y_index5:x2 -0.0456199 0.0453805 -1.005 Y_index6:x2 0.0333067 0.0353716 0.942 Y_index7:x2 0.0443440 0.0371040 1.195 Y_index8:x2 0.0595453 0.0353532 1.684 Y_index9:x2 0.0223958 0.0379038 0.591 What I don't understand is (1) why Y_index1 is missing in my interaction output with x1 and x2 (the lists start with Y_index2), and (2) what the interaction lines exactly mean. Is it a comparison to a baseline? Is it difference from the main effect? I do apologize in case this is a very simple question, but I cannot get my head around it. Thank you --E [[alternative HTML version deleted]] </pre> Eiko Fried 2012-05-20T11:53:00 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8211 <pre>Dear R-listers, I have a couple of questions regrading how to assess the effects of lower order terms when a higher order term is significant using likelihood ratio tests. I am interested to examine how two independent variables (iv1 and iv2) affect subjects' log-transformed reaction time (logRT). A subset of the data can be accessed in the link below if needed: http://www-scf.usc.edu/~hex/testset.txt In a series of responses that Dr. Bates made in this thread: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005399.html, he recommended that non-significant fixed effects be removed from a model before further llkelihood ratio tests. For example, I can first remove *iv2* and compare the model with both iv1 and iv2 with the model without iv2. model.main = lmer(logRT ~ iv1 + iv2 + (1 + iv2 | subject) + (1 | item), Sample, REML = FALSE) model.iv1 = lmer(logRT ~ iv1 + (1 + iv2 | subject) + (1 | item), Sample, REML = FALSE) anova(model.main, model.iv1) #Data: Sample #Models: #model.iv1: logRT ~ iv1 + (1 + iv2 | subject) + (1 | item) #model.main: logRT ~ iv1 + iv2 + (1 + iv2 | subject) + (1 | item) # Df AIC BIC logLik Chisq Chi Df Pr(&gt;Chisq) #model.iv1 7 235.58 266.31 -110.79 #model.main 8 236.31 271.43 -110.15 1.2754 1 0.2588 In this case, there is no significant difference between the two models. Based on Dr. Bates' suggestion, if I would like to further examine the contribution of iv1, I could start from model.iv1 by comparing model.iv1 with a null model, as shown below: model.null = lmer(logRT ~ 1 + (1 + iv2 | subject) + (1 | item), Sample, REML = FALSE) anova(model.iv1, model.null) #Data: Sample #Models: #model.null: logRT ~ 1 + (1 + iv2 | subject) + (1 | item) #model.iv1: logRT ~ iv1 + (1 + iv2 | subject) + (1 | item) # Df AIC BIC logLik Chisq Chi Df Pr(&gt;Chisq) #model.null 6 235.99 262.33 -112.00 #model.iv1 7 235.58 266.31 -110.79 2.4084 1 0.1207 As is shown, iv1 doesn't seem to improve to the fit of the model significantly either. While this procedure is clear, how does one deal with a situation where, say, the full model has a higher order term, such as a two-way interaction, and one would like to know not only how much the interaction contributes to the model fit but also how much the lower order terms contribute to the model. model.inter = lmer(logRT ~ iv1 * iv2 + (1 + iv2 | subject) + (1 | item), Sample, REML = FALSE) anova(model.inter, model.main) #Data: Sample #Models: #model.main: logRT ~ iv1 + iv2 + (1 + iv2 | subject) + (1 | item) #model.inter: logRT ~ iv1 * iv2 + (1 + iv2 | subject) + (1 | item) # Df AIC BIC logLik Chisq Chi Df Pr(&gt;Chisq) #model.main 8 236.31 271.43 -110.15 #model.inter 9 229.10 268.61 -105.55 9.2097 1 0.002407 ** #--- #Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (1). My first question is as follows: For example, in the case above, the interaction is shown as significant. I know that when an interaction is present, main effects are often not interpretable, so probably that makes it unnecessary to even assess main effects? But say, if for whatever reason, I need to assess lower order terms - for example iv2, should I compare the first pair of models below, or the second pair of models below? It seems to me that comparing the 2nd pair makes more sense because the significant interaction probably shouldn't be removed from the model. First pair: model.main: logRT ~ iv1 + iv2 + (1 + iv2 | subject) + (1 | item) model.iv1: logRT ~ iv1 + (1 + iv2 | subject) + (1 | item) Second pair model.main: logRT ~ iv1 + iv2 + iv1:iv2 (1 + iv2 | subject) + (1 | item) model: logRT ~ iv1 + iv1:iv2 + (1 + iv2 | subject) + (1 | item) (2). My second question has to do with random effect specification. As you can see, in all the models I showed above, I have iv2 as the random slope - which I decided based upon comparing models that have the same fixed effects ( iv1 * iv2) but have different nested random effect specifications. While including iv2 makes sense for models that have iv2 as a fixed effect (e.g., model.main, model.inter), its inclusion doesn't seem to make sense in models where iv2 is not included as a fixed effect - for example, the null model and model.iv1. However, it seems necessary to include it when I compare the models below. I wonder if someone can tell me what is the proper way of dealing with random effects in this kind of scenario and the reasons behind. model.main: logRT ~ iv1 + iv2 + (1 + iv2 | subject) + (1 | item) model.iv1: logRT ~ iv1 + (1 + iv2 | subject) + (1 | item) Thank you in advance!! Best, Xiao [[alternative HTML version deleted]] </pre> Xiao He 2012-05-19T14:27:40 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8210 <pre>Dear All I wonder if anyone can help with this issue.?    I have two fixed factors (experimental group and sex) and a random factor (brood). I need the random factor to control for pseudoreplication because each brood provides a?    certain number of siblings, which are pseudoreplicates because share genes/environment. Then I have 8 variables. so I would like to to compare the covariance structures among groups while controlling for the random factor. I was told that R might do that, I mean a manova with mixed effects (manova(lme(....))?    thanks for your help?    David ______________________________________________ David Costantini, PhD http://www.davidcostantini.it NERC Postdoctoral research associate Institute of Biodiversity, Animal Health and Comparative Medicine School of Life Sciences College of Medical, Veterinary and Life Sciences University of Glasgow Graham Kerr Building, room 511 Glasgow G12 8QQ, UK See also my association Ornis italica http://www.ornisitalica.com http://www.birdcam.it ____________________________________________________ [[alternative HTML version deleted]] </pre> David Costantini 2012-05-19T09:16:15 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8209 <pre>Hello mixed modelers, I am having problems with some GEE models I am trying to run using geepack. I have species abundance data for 52 different species in 154 sites over 47 years, and I am trying to extract slope parameter estimates so that I can look at whether these species have increased or decreased in abundance over time, while taking into account the repeated measurements at each site over time. I originally started doing this with mixed models, but have been advised that GEE would be more appropriate for my data as it gives population-averaged responses. However, when I try to run GEE's on my data I get really bizarre parameter estimates for some of my species. As my dataset is huge I unfortunately cannot provide the whole thing, but I have uploaded a subset of the data for one species with a particularly bizarre slope parameter estimate here: http://dl.dropbox.com/u/4481861/Example_for_GEE_one_species.csv The data look like this: Site Year Species Value_Pres Value_Abs 1 1 1961 1 0 2089 2 1 1962 1 0 2120 3 1 1963 1 0 2089 4 1 1964 1 0 2225 5 1 1965 1 0 2197 6 1 1966 1 0 2208 I have been using the following model specification (I have been running a loop to calculate estimates for all 52 species separately, but this is for just one species): speciesA&lt;-orderBy(~Site+Year,data=speciesA) #using the doBy package to order by subject then time speciesA.mod&lt;-geeglm(cbind(Value_Pres,Value_Abs)~I(Year-1961),data=speciesA, family=binomial,id=Site,corstr="ar1") Call: geeglm(formula = cbind(Value_Pres, Value_Abs) ~ I(Year - 1961), family = binomial, data = speciesA, id = Site, corstr = "ar1") Coefficients: Estimate Std.err Wald Pr(&gt;|W|) (Intercept) -2.99e+14 9.10e+11 107705 &lt;2e-16 *** I(Year - 1961) -9.62e+13 3.88e+10 6155147 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Estimated Scale Parameters: Estimate Std.err (Intercept) 6.57e+10 5.62e+30 Correlation: Structure = ar1 Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.98 4.4e+18 Number of clusters: 154 Maximum cluster size: 47 I suspect the problem might have something to do with the correlation structure, as species abundances in subsequent years are often very highly correlated, even if there is substantial change over the 47 years overall. If I use the corstr="independence" command I get parameter estimate that are very similar to those I got using mixed effects models (at least, the slopes for species responses relative to each other are similar). Furthermore, if I use corstr="ar1" but subset my data to every 5 years instead of every year, I get much more reasonable slope estimates for this particular species as well as most of the other species (slope values are very similar to the corstr="independence" value), but a few different species' slopes then get very weird. (By get weird I mean that they have abnormally large positive or negative slopes that don't reflect what's happening in the raw data at all). I would really appreciate some insight into what the problem with my data could be, or, more particularly, how to fix it! My head and my wall would be very grateful! Perhaps I should just give up on GEE's and go back to mixed models?? Thanks very much, Anne [[alternative HTML version deleted]] </pre> Anne Bjorkman 2012-05-18T18:19:28 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8208 <pre>Dear list members, I would be very grateful for any advice on fitting a model for a cross-over design please excuse my lengthy post. I am currently trying to analyze a transient treatment effect (treatment vs. placebo) in a cross-over experiment (2 sessions) in two groups (old, young subjects) with repeated measurements for both treatment conditions (before treatment=baseline, post, post45, post90) in multiple outcome variables separately (2 physiological and 4 behavioral, all of which metric). Points of interest are whether there is an effect of treatment in both/either of the groups in one/more outcome variables and if an inter-correlation between potential changes in physiological and behavioral outcome exists. Problems I am trying to deal with: Although a carry-over of direct treatment effect seems unlikely (sufficient wash-out phase in-between two conditions/sessions), a possible change (improvement) in response rate from 1st to 2nd session irrespective of treatment, i.e. rather "period effect", cannot be ruled out Is it sufficient to fit a random intercept to tackle this problem? OR is it necessary to use baseline adjusted outcome? Is there any OTHER option? Fitting the model (example attached), I am using separate models for each outcome, correlation structure corAR1, random INTERCEPT, and random SLOPE (time point of repeated measurements). I sit necessary to specify the random effects for grouping factors, i.e. for treatment and/or agegroup? ~time point|ID OR |treatment/~time pointID?, and in case of the latter, fitting factor treatment also as fixed? The first time point is baseline, i.e. before onset of treatment effect, and the last time point of measurement is expected to already reflect wash-out of treatment. I read some advice on coding time for repeated measurements but need some help with what might be the most appropriate/accepted way to do it: time as continuous variable (in minutes relative to treatment onset, here: time_point: -15, 20, 45, 90) time point as categorical or ordered? does it make sense to analyze all time points as categorical in a first step, to see whether there is a significant "treatment*time" interaction for specified time points, and then in a second step to analyze slope from baseline to respective time point, i.e. code time as continuous? Is there a reasonable way of quantifying a possible treatment effect in order to analyze relation of potential effects on outcome, i.e. effect on physiological outcome associated with effect on behavioral outcome? I tried to find more detailed advice ( eg. Dı́az-Uriarte, R. (2002). Incorrect analysis of crossover trials in animal behaviour research. Animal Behaviour, 63(4), 815–822. doi:10.1006/anbe.2001.1950) and would appreciate any further recommendations. Thank you very much in advance! kirstin --- data.frame':5760 obs. of 21 variables: $ ID [=subject level] : Factor w/ 20 levels "OMI_01","OMI_02",..: 11 11 11 11 11 11 11 11 11 11 ... $ agegroup : Factor w/ 2 levels "OLD ","YOUNG ": 2 2 2 2 2 2 2 2 2 2 ... $ session : num 1 1 1 1 1 1 1 1 1 1 ... $ stim.cond [=treatment]: Factor w/ 2 levels "sham","tDCS": 1 1 1 1 1 1 1 1 1 1 ... $ timepoint : Factor w/ 4 levels "BL ","P ","P45",..: 1 1 1 1 1 1 1 1 1 1 ... $ time_point : num -15 -15 -15 -15 -15 -15 -15 -15 -15 -15 ... $ outcome : num NA 2.887 0.963 4.006 2.06 ... [....] example for model, BIC comparison favors the corAR1 correlation structure: summary(mR0x_corAR &lt;- lme(outcome~factor(agegroup)*factor(timepoint)*factor(treatment),data=tDCSrest, random=~1|ID, correlation = corAR1(), na.action=na.exclude,method="REML")) Linear mixed-effects model fit by REML Data: tDCSrest AIC BIC logLik 5827.065 5936.875 -2894.533 Random effects: Formula: ~1 | ID (Intercept) Residual StdDev: 0.7819433 0.8284249 Correlation Structure: AR(1) Formula: ~1 | ID Parameter estimate(s): Phi 0.3094027 Fixed effects: outcome ~ factor(agegroup) * factor(timepoint) * factor(stim.cond) Value Std.Error DF t-value p-value (Intercept) 3.656386 0.2625337 2373 13.927301 0.0000 factor(agegroup)young -0.530390 0.3710212 18 -1.429540 0.1700 factor(timepoint)P -0.420957 0.1222664 2373 -3.442950 0.0006 factor(timepoint)P45 -0.120908 0.1265935 2373 -0.955090 0.3396 factor(timepoint)P90 -0.291501 0.1252295 2373 -2.327736 0.0200 factor(stim.cond)sham -0.253259 0.1292979 2373 -1.958726 0.0503 factor(agegroup)young:factor(timepoint)P 0.508492 0.1733311 2373 2.933646 0.0034 factor(agegroup)young:factor(timepoint)P45 0.372538 0.1768388 2373 2.106652 0.0353 factor(agegroup)young:factor(timepoint)P90 0.380923 0.1760245 2373 2.164036 0.0306 factor(agegroup)young:factor(stim.cond)sham 0.284941 0.1792027 2373 1.590047 0.1120 factor(timepoint)P :factor(stim.cond)sham 0.587820 0.1777151 2373 3.307656 0.0010 factor(timepoint)P45:factor(stim.cond)sham 0.050355 0.1824015 2373 0.276066 0.7825 factor(timepoint)P90:factor(stim.cond)sham 0.430209 0.1832581 2373 2.347559 0.0190 factor(agegroup)young:factor(timepoint)P :factor(stim.cond)sham -0.785629 0.2482133 2373 -3.165136 0.0016 factor(agegroup)young:factor(timepoint)P45:factor(stim.cond)sham -0.383033 0.2539655 2373 -1.508210 0.1316 factor(agegroup)young:factor(timepoint)P90:factor(stim.cond)sham -0.478725 0.2549540 2373 -1.877691 0.0605 Correlation: (Intr) fctr() fct()P fc()P45 fc()P90 fct(.) fc():()P factor(agegroup)young -0.708 factor(timepoint)P -0.233 0.165 factor(timepoint)P45 -0.234 0.165 0.503 factor(timepoint)P90 -0.237 0.167 0.489 0.510 factor(stim.cond)sham -0.229 0.162 0.474 0.476 0.480 factor(agegroup)young:factor(timepoint)P 0.164 -0.230 -0.705 -0.355 -0.345 -0.334 factor(agegroup)young:factor(timepoint)P45 0.167 -0.234 -0.360 -0.716 -0.365 -0.341 0.501 factor(agegroup)young:factor(timepoint)P90 0.168 -0.235 -0.348 -0.363 -0.711 -0.342 0.485 factor(agegroup)young:factor(stim.cond)sham 0.165 -0.231 -0.342 -0.344 -0.347 -0.722 0.476 factor(timepoint)P :factor(stim.cond)sham 0.160 -0.113 -0.688 -0.346 -0.337 -0.690 0.485 factor(timepoint)P45:factor(stim.cond)sham 0.162 -0.115 -0.350 -0.695 -0.355 -0.699 0.247 factor(timepoint)P90:factor(stim.cond)sham 0.168 -0.119 -0.347 -0.362 -0.696 -0.708 0.245 factor(agegroup)young:factor(timepoint)P :factor(stim.cond)sham -0.115 0.160 0.493 0.248 0.241 0.494 -0.698 factor(agegroup)young:factor(timepoint)P45:factor(stim.cond)sham -0.117 0.163 0.251 0.499 0.255 0.502 -0.350 factor(agegroup)young:factor(timepoint)P90:factor(stim.cond)sham -0.121 0.169 0.250 0.260 0.501 0.509 -0.347 fc():()P45 fc():()P90 f():(. f()P:( f()P45: f()P90: factor(agegroup)young factor(timepoint)P factor(timepoint)P45 factor(timepoint)P90 factor(stim.cond)sham factor(agegroup)young:factor(timepoint)P factor(agegroup)young:factor(timepoint)P45 factor(agegroup)young:factor(timepoint)P90 0.513 factor(agegroup)young:factor(stim.cond)sham 0.485 0.487 factor(timepoint)P :factor(stim.cond)sham 0.248 0.240 0.498 factor(timepoint)P45:factor(stim.cond)sham 0.497 0.252 0.505 0.508 factor(timepoint)P90:factor(stim.cond)sham 0.259 0.495 0.511 0.496 0.521 factor(agegroup)young:factor(timepoint)P :factor(stim.cond)sham -0.350 -0.338 -0.689 -0.716 -0.364 -0.355 factor(agegroup)young:factor(timepoint)P45:factor(stim.cond)sham -0.697 -0.358 -0.700 -0.365 -0.718 -0.374 factor(agegroup)young:factor(timepoint)P90:factor(stim.cond)sham -0.367 -0.704 -0.711 -0.357 -0.375 -0.719 f():(: f():()P45: factor(agegroup)young factor(timepoint)P factor(timepoint)P45 factor(timepoint)P90 factor(stim.cond)sham factor(agegroup)young:factor(timepoint)P factor(agegroup)young:factor(timepoint)P45 factor(agegroup)young:factor(timepoint)P90 factor(agegroup)young:factor(stim.cond)sham factor(timepoint)P :factor(stim.cond)sham factor(timepoint)P45:factor(stim.cond)sham factor(timepoint)P90:factor(stim.cond)sham factor(agegroup)young:factor(timepoint)P :factor(stim.cond)sham factor(agegroup)young:factor(timepoint)P45:factor(stim.cond)sham 0.506 factor(agegroup)young:factor(timepoint)P90:factor(stim.cond)sham 0.494 0.520 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.613364210 -0.622035919 -0.001779787 0.637882645 4.410916182 Number of Observations: 2407 Number of Groups: 20 anova(mR0x_corAR) numDF denDF F-value p-value (Intercept) 1 2373 353.1295 &lt;.0001 factor(agegroup) 1 18 0.6113 0.4445 factor(timepoint) 3 2373 0.6737 0.5681 factor(stim.cond) 1 2373 1.1819 0.2771 factor(agegroup):factor(timepoint) 3 2373 0.8459 0.4686 factor(agegroup):factor(stim.cond) 1 2373 1.7029 0.1920 factor(timepoint):factor(stim.cond) 3 2373 3.3753 0.0177 factor(agegroup):factor(timepoint):factor(stim.cond) 3 2373 3.4048 0.0170 -- Pflichtangaben gemäß Gesetz über elektronische Handelsregister und Genossenschaftsregister sowie das Unternehmensregister (EHUG): Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg Vorstandsmitglieder: Prof. Dr. Guido Sauter (Vertreter des Vorsitzenden), Dr. Alexander Kirstein, Joachim Prölß, Prof. Dr. Dr. Uwe Koch-Gromus [[alternative HTML version deleted]] </pre> Kirstin-Friederike Heise 2012-05-18T13:13:04 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8206 <pre>Greetings R users, I am trying to match some SAS output of a mixed model. After running PROC MIXED with the covariance structure as AR1 the output below is given. Now when I try to replicate this with lme, I get the correct degrees of freedom and somewhat close values. If I try with gls, I get the exact F value but the incorrect denominator degrees of freedom. Is there some syntax or parameter I can adjust to get lme to have the same F values as gls? That way the correct degrees of freedom would be applied and I would get the replicated model. The data is also included below. Thanks to any who can help with this issue. Regards, Charles Row ID Group Died Event_name var 1 12510 3 YES B -1.05257 2 12510 3 YES S45 -1.00000 3 12510 3 YES FR2 -1.14630 4 12510 3 YES FR8 -1.08831 5 12510 3 YES FR20 -1.03339 6 21510 3 NO B -0.87290 7 21510 3 NO S45 -1.22185 8 21510 3 NO FR2 -1.01592 9 21510 3 NO FR8 -1.06550 10 21510 3 NO PR48 -0.65326 11 22210 3 YES B -1.37059 12 30810 3 NO B -0.92959 13 30810 3 NO S45 -1.26680 14 30810 3 NO FR2 -0.99396 15 30810 3 NO FR8 -1.04769 16 30810 3 NO FR20 -1.04769 17 31510 4 YES B -1.50585 18 31510 4 YES S45 -1.35262 19 31510 4 YES FR2 -1.26520 20 31510 4 YES FR8 -1.11407 21 31510 4 YES FR20 -0.97757 22 40510 2 YES B -1.04721 23 40510 2 YES S45 -1.13371 24 51010 3 NO B -1.21467 25 51010 3 NO S45 -1.18177 26 51010 3 NO FR2 -1.23582 27 51010 3 NO FR8 -1.35262 28 51010 3 NO FR20 -0.69897 29 51010 3 NO PR48 -1.14086 30 51710 4 NO B -1.13371 31 51710 4 NO S45 -1.24336 32 51710 4 NO FR2 -0.85201 33 51710 4 NO FR8 -0.72700 34 51710 4 NO FR20 -0.71220 35 51710 4 NO PR48 -1.10735 36 51809 1 NO B -0.70752 37 51809 1 NO FR2 -0.58120 38 51809 1 NO FR8 -0.71332 39 51809 1 NO FR20 -0.37212 40 61410 4 NO B -1.25493 41 61410 4 NO S45 -1.24949 42 61410 4 NO FR2 -1.18310 43 61410 4 NO FR8 -0.87517 44 61410 4 NO FR20 -0.34679 45 61410 4 NO PR48 -0.76472 46 62110 3 NO B -1.34199 47 62110 3 NO S45 -1.23657 48 62110 3 NO FR2 -1.20412 49 62110 3 NO FR8 -0.72011 50 62110 3 NO FR20 -0.72308 51 62110 3 NO PR48 -0.89688 52 62209 3 NO B -0.74041 53 62209 3 NO S45 -0.99183 54 62209 3 NO FR2 -1.06854 55 62209 3 NO FR8 -0.94386 56 62209 3 NO FR20 -0.38817 57 71210 4 YES B -1.23657 58 71210 4 YES S45 -0.94808 59 71309 4 YES B -1.14630 60 71309 4 YES S45 -1.10568 61 71309 4 YES FR2 -0.94462 62 71309 4 YES FR8 -0.89211 63 71309 4 YES FR20 -0.98970 64 71309 4 YES PR48 -0.83387 65 81009 4 YES B -1.15989 66 81009 4 YES S45 -0.93779 67 81009 4 YES FR2 -0.86710 68 81009 4 YES FR8 -0.69379 69 81610 1 NO B -1.04769 70 81610 1 NO S45 -1.02780 71 81610 1 NO FR2 -0.96098 72 81610 1 NO FR8 -0.91829 73 81610 1 NO FR20 -0.68550 74 81610 1 NO PR48 -0.66374 75 81709 4 YES B -0.65521 76 82410 3 NO B -1.06702 77 82410 3 NO S45 -1.30103 78 82410 3 NO FR2 -0.96617 79 82410 3 NO FR8 -0.84497 80 82410 3 NO FR20 -0.47716 81 82410 3 NO PR48 -0.68131 82 91310 4 NO B -1.30803 83 91310 4 NO S45 -1.07212 84 91310 4 NO FR2 -0.94615 85 91310 4 NO FR8 -1.29414 86 91310 4 NO FR20 -0.81304 87 91310 4 NO PR48 -0.86519 88 91409 4 NO B -1.12321 89 91409 4 NO S45 -1.08831 90 91409 4 NO FR2 -1.11862 91 91409 4 NO FR8 -0.87976 92 91409 4 NO FR20 -0.87517 93 91409 4 NO PR48 -0.68867 94 92109 3 NO B -1.15120 95 92109 3 NO S45 -0.99525 96 92109 3 NO FR2 -0.95429 97 92109 3 NO FR8 -0.54914 98 92109 3 NO FR20 -0.67264 99 92109 3 NO PR48 -0.74256 100 92110 4 YES B -1.05012 101 92110 4 YES S45 -1.14086 102 92110 4 YES FR2 -0.91901 103 92110 4 YES FR8 -0.80521 104 92110 4 YES FR20 -0.88807 105 92710 4 NO B -1.20412 106 92710 4 NO S45 -1.21325 107 92710 4 NO FR2 -1.19518 108 92710 4 NO FR8 -0.92959 109 92710 4 NO FR20 -0.69315 110 92710 4 NO PR48 -0.73779 111 100410 1 NO B -0.92300 112 100410 1 NO S45 -1.07314 113 100410 1 NO FR2 -0.80162 114 100410 1 NO FR8 -0.63283 115 100410 1 NO FR20 -0.73119 116 100410 1 NO PR48 0.10102 117 101209 3 NO B -1.24718 118 101209 3 NO S45 -1.00524 119 101209 3 NO FR2 -1.01592 120 101209 3 NO FR8 -0.85201 121 101209 3 NO FR20 -0.51670 122 101209 3 NO PR48 -0.72239 123 101909 3 NO B -1.13549 124 101909 3 NO S45 -0.92665 125 101909 3 NO FR2 -0.98548 126 101909 3 NO FR8 -1.03152 127 101909 3 NO FR20 -0.71715 128 101909 3 NO PR48 -1.07935 129 102510 2 YES B -1.30103 130 102510 2 YES S45 -0.83387 131 110810 3 YES B -1.21968 132 110810 3 YES S45 -0.68048 133 110810 3 YES FR2 -1.09474 134 110909 2 YES B -1.09259 135 110909 2 YES S45 -0.91293 136 111510 3 NO B -1.39794 137 111510 3 NO S45 -1.42366 138 111510 3 NO FR2 -1.22621 139 111510 3 NO FR8 -0.98885 140 111510 3 NO FR20 -0.66055 141 111510 3 NO PR48 -0.84497 142 111609 1 NO B -0.79942 143 111609 1 NO S45 -0.87517 144 111609 1 NO FR2 -0.74041 145 111609 1 NO FR8 -0.64092 146 111609 1 NO FR20 -0.53835 147 111609 1 NO PR48 -0.68719 148 120610 4 NO B -1.05110 149 120610 4 NO S45 -1.45100 150 120610 4 NO FR2 -1.21183 151 120610 4 NO FR8 -1.04144 152 120610 4 NO FR20 -0.61101 153 120610 4 NO PR48 -0.81702 154 120709 1 NO B -0.84497 155 120709 1 NO S45 -0.99439 156 120709 1 NO FR2 -1.13430 157 120709 1 NO FR20 -0.57187 158 120709 1 NO PR48 -0.80632 159 121310 3 YES B -1.16115 160 121310 3 YES S45 -1.18310 dat=read.table("C:/…/subset.csv",sep=",",header=TRUE, na.strings=".") attach(dat) dat34=dat[Pig_group %in% c("3", "4"),] attach(dat34) dat34=within(dat34, { group=factor(group) Event_name=factor(Event_name) Died=factor(Died) ID=factor(ID) }) attach(dat34) contrasts(dat34$Event_name)=contr.sum(n=6) contrasts(dat34$group)=contr.sum(n=2) contrasts(dat34$Died)=contr.sum(n=2) # What the results should be from SAS AR1 var #Type 3 Tests of Fixed Effects #Effect NumDF DenDF F Value Pr &gt; F #Event_name 5 91 7.43 &lt;.0001 #Died 1 24 0.08 0.7756 #Event_name*Died 5 91 3.09 0.0128 fit.13=lme(var~Event_name*Died, data=liver34, random=~1|ID, corr=corAR1()) anova(fit.13, type="marginal", adjustSigma=F) # numDF denDF F-value p-value #(Intercept) 1 91 1342.7364 &lt;.0001 #Event_name 5 91 8.7143 &lt;.0001 #Died 1 24 0.0527 0.8204 #Event_name:Died 5 91 3.0909 0.0127 ###################################################### ###THIS DOES IT but not lme function (need dfs) fit.16=gls(var~Event_name*Died, data=liver34, corr=corAR1(, ~1|ID)) anova(fit.16, type="marginal", adjustSigma=F) #Denom. DF: 115 # numDF F-value p-value #(Intercept) 1 1496.9032 &lt;.0001 #Event_name 5 7.4304 &lt;.0001 #Died 1 0.0831 0.7736 #Event_name:Died 5 3.0871 0.0119 #Give exact replicated answers 1-pf(7.4304, 5, 91) 1-pf(.0831, 1, 24) 1-pf(3.0871, 5, 91) ####################################################### [[alternative HTML version deleted]] </pre> Charles Determan Jr 2012-05-18T13:21:53 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8205 <pre>Hi folks, This (http://stats.stackexchange.com/questions/24293) recent-ish post to stats.stackexchange.com suggests that when it comes to tools available in R, one cannot fit a non-linear mixed effects model with binomial error. I just wanted to double check that this is an accurate portrayal of the state of the aRt. (In case it matters, I'm looking to follow up the answer to my question here [http://stats.stackexchange.com/questions/22895]) Cheers, Mike -- Mike Lawrence Graduate Student Department of Psychology Dalhousie University Looking to arrange a meeting? Check my public calendar: http://goo.gl/BYH99 ~ Certainty is folly... I think. ~ </pre> Mike Lawrence 2012-05-18T13:02:48 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8202 <pre>Hi All, I am trying to estimate disease heritability ultimately, but to get started, I wanted to model a continuous outcome using lmekin(). I tried something very similar (I think) to example 3 (GAW) of the lmekin vignette from the coxme package http://cran.r-project.org/web/packages/coxme/vignettes/lmekin.pdf However, I get an error about the variance matrix not having dimnames. Although there does not seem to be data provided with the vignette, I believe I should be essentially perfectly replicating the steps. Any ideas on what to do? Below is made up data and an example along with my session info. Thanks! Josh ################################################ require(kinship2) adat &lt;- data.frame(PersonID = 1:12, FatherID = c(NA, NA, NA, 1, 1, 1, 1, 4, 4, 4, 4, 4), MotherID = c(NA, NA, NA, 2, 2, 2, 2, 3, 3, 3, 3, 3), sex = c(0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1)) adat &lt;- do.call("rbind", rep(list(adat), 100)) adat$FamilyID &lt;- rep(1:100, each = 12) set.seed(10) adat$Outcome &lt;- rnorm(nrow(adat), mean = adat$sex + rep(rnorm(100), each = 12)) gadat &lt;- with(adat, pedigree(PersonID, FatherID, MotherID, sex, famid = FamilyID)) ## plot family structure (identical for all) plot(gadat[1]) kamat &lt;- kinship(gadat) require(coxme) (gfit1 &lt;- lmekin(Outcome ~ factor(sex) + (1 | FamilyID), data=adat, varlist=kamat)) # Error in bdsmatrix.reconcile(varlist, bname) : # No dimnames found on a variance matrix Here is my sessionInfo(): R Under development (unstable) (2012-05-08 r59331) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] splines stats graphics grDevices utils datasets methods [8] base other attached packages: [1] coxme_2.2-2 nlme_3.1-103 bdsmatrix_1.3 survival_2.36-14 [5] kinship2_1.3.3 quadprog_1.5-4 Matrix_1.0-6 lattice_0.20-6 loaded via a namespace (and not attached): [1] compiler_2.16.0 grid_2.16.0 tools_2.16.0 </pre> Joshua Wiley 2012-05-18T01:00:11 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8189 <pre>*Hi all, it is counterintuitive to me that adding a random effect to a GLMM should lead to a drop in the likelihood of the fitted model for, after all, one can give the new random effect vanishingly small variance which should lead to the same likelihood as before. But this is what I have encountered. I regret that the output below is not self contained but if needed I could follow-up offline with additional files. The factor "box" is nested inside the factor "tree" and the factor "gap" is crossed with both of these. This is a binomial GLMM with 0/1 response. (1|tree), family = binomial) Generalized linear mixed model fit by the Laplace approximation Formula: fincr ~ icfac + (1 | gap) + (1 | box) + (1 | gap:box) + (1 | tree) AIC BIC logLik deviance 687 732.4 -335.5 671 Random effects: Groups Name Variance Std.Dev. gap:box (Intercept) 0.000000 0.00000 box (Intercept) 0.082784 0.28772 gap (Intercept) 0.122614 0.35016 tree (Intercept) 0.505608 0.71106 Number of obs: 2160, groups: gap:box, 1080; box, 40; gap, 27; tree, 20 Fixed effects: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -3.4348 0.2315 -14.835 &lt; 2e-16 *** icfacmale 1.4535 0.2738 5.308 1.11e-07 *** icfacfem -2.5740 0.7480 -3.441 0.000579 *** icfacmix 0.2935 1.1049 0.266 0.790517 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) icfcml icfcfm icfacmale -0.357 icfacfem -0.114 0.122 icfacmix -0.065 0.066 0.032 Data: Models: f0128bi: fincr ~ icfac + (1 | gap) + (1 | box) + (1 | gap:box) f0128bit: fincr ~ icfac + (1 | gap) + (1 | box) + (1 | gap:box) + (1 | f0128bit: tree) Df AIC BIC logLik Chisq Chi Df Pr(&gt;Chisq) f0128bi 7 615.73 655.48 -300.87 f0128bit 8 687.01 732.43 -335.51 0 1 1 Any suggestions gratefully received. Murray* </pre> Murray Jorgensen 2012-05-16T22:35:02 http://comments.gmane.org/gmane.comp.lang.r.lme4.devel/8187 <pre>I adjusted a couple of models with different packages and their respective roles I think these models are equivalent, ie they have the same fixed and random structures   exp.nlme&lt;- nlme (dL~expT(largo, IncM,k),      data = dL.gdOK,      fixed = IncM +k~ 1,      random= IncM+k~1,      start = list(fixed = c(11.2,0.047)))   exp.nlmer2&lt;-nlmer(dL~expT(largo, IncM,k)~(IncM+k|Codigo), dL.gdOK,     start=c(IncM=11.2,k=0.047))   However, the output shows different values ​​for: estimates of fixed and random effects, variances and standard deviations of random effects, log-likelihoods…   Nonlinear mixed-effects model fit by maximum likelihood   Model: dL ~ expT(largo, IncM, k)   Data: dL.gdOK   Log-likelihood: -2559.485   Fixed: IncM + k ~ 1        IncM           k 11.81793370  0.04770045   Random effects:  Formula: list(IncM ~ 1, k ~ 1)  Level: Codigo  Structure: General positive-definite, Log-Cholesky parametrization          StdDev     Corr IncM     5.17079862 IncM k        0.01189968 0.6 Residual 0.95701497       Number of Observations: 1677 Number of Groups: 118     Nonlinear mixed model fit by the Laplace approximation Formula: dL ~ expT(largo, IncM, k) ~ (IncM + k | Codigo)    Data: dL.gdOK   AIC  BIC logLik deviance  2052 2085  -1020     2040 Random effects:  Groups   Name Variance   Std.Dev. Corr   Codigo   IncM 2.3700e+01 4.868270                 k    1.2203e-04 0.011047 0.543  Residual      9.2676e-01 0.962685       Number of obs: 1677, groups: Codigo, 118   Fixed effects:       Estimate Std. Error t value IncM 12.108597   0.489472   24.74 k     0.049541   0.001274   38.87   Correlation of Fixed Effects:   IncM k 0.600   Could someone explain to me why these differences occur? Could anyone tell me what causes these differences? I have read enough material on nonlinear mixed effects but can not find the answer to my problem. Cheers Gabriela Gabriela Escati Peñaloza Centro Nacional Patagónico(CENPAT). CONICET Bvd. Brown s/nº. (U9120ACV)Pto. Madryn Chubut Argentina Tel: 54-2965/451301/451024/451375/45401 (Int:277) Fax: 54-29657451543 [[alternative HTML version deleted]] </pre> gabriela escati peñaloza 2012-05-16T14:00:59 Search the mailing list at Gmane query http://search.gmane.org/?group=$group=gmane.comp.lang.r.lme4.devel