Joshua, this is amazingly simple and works!!! Thank you so much!

Jason


 
Dr. Iasonas Lamprianou
Department of Social and Political Sciences
University of Cyprus


>________________________________
> From: Joshua Wiley <jwiley.psych@gmail.com>
>To: Iasonas Lamprianou <lamprianou@yahoo.com> 
>Cc: "r-sig-mixed-models@r-project.org" <r-sig-mixed-models@r-project.org> 
>Sent: Monday, 1 April 2013, 0:19
>Subject: Re: [R-sig-ME] lme heteroscedasticity
> 
>Dear Iasonas,
>
>Consider this simple example:
>
>mtcars <- mtcars
>mtcars$weights <- sample(0:4, size = nrow(mtcars), replace=TRUE)
>lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~weights), data = mtcars)
>
>which gives many errors.  But...
>
>lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+1)),
>data = mtcars)
>
>and even
>
>lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+.001)),
>data = mtcars)
>
>but not
>
>lme(mpg ~ 1, random = ~ 1 | cyl, weights = varFixed(~(weights+0)),
>data = mtcars)
>
>
>Anyway, the short of it is that having 0 weights does not play well
>with lme.  They do not have to be far away from zero, but it does not
>like them to be exactly zero.
>
>In your case, you have count of days absent, with most students having
>zero.  If that is the baseline, and then there is more variability in
>the students who are absent at least some days, since the weights are
>a multiplier essentially of the residual variance, consider adding 1,
>which then would make the residual variance for no absent, and then it
>would increase based on number of days absent.
>
>Does that make sense?
>
>Cheers,
>
>Joshua
>
>
>
>On Sun, Mar 31, 2013 at 1:07 AM, Iasonas Lamprianou
><lamprianou@yahoo.com> wrote:
>>
>>
>> Dear all, a few days ago I sent an email to this list and got some responses (I thank you all). This time, I am trying to fit an heteroscedastic model  on lme:
>>
>> model0<-lme(score ~ 1  + overallabs +PARENTALEDUCATIONALLEVEL + PARENTALOCCUPATIONALLEVEL + GENDER +  SUBJECT *  ETHNICITY   +  SCHOOL2 ,random=~1|id, weights=varFixed(~overallabs),  data=galatia.comb_small.vert)
>>
>> overallabs is the number of absences from school (which is a heavily skewed variable; most students have zero absences) and this causes some heteroscedasticity in the model. id is the id of the students (each one has scores on 4 subjects which is modelled as fixed effects)
>>
>>
>> However, when I try this model, I get the result:
>>
>> Error in MEestimate(lmeSt, grps) :
>>   NA/NaN/Inf in foreign function call (arg 1)
>>
>> Not sure why this is the case. Anyone any ideas? The internet cannot help me in this case....
>>
>> I am using the laterst nlme version on linux mint.
>>
>> Thanks
>>
>>
>> Dr. Iasonas Lamprianou
>> Department of Social and Political Sciences
>> University of Cyprus
>>
>>
>>>________________________________
>>> From: "r-sig-mixed-models-request@r-project.org" <r-sig-mixed-models-request@r-project.org>
>>>To: r-sig-mixed-models@r-project.org
>>>Sent: Friday, 29 March 2013, 8:23
>>>Subject: R-sig-mixed-models Digest, Vol 75, Issue 45
>>>
>>>Send R-sig-mixed-models mailing list submissions to
>>>    r-sig-mixed-models@r-project.org
>>>
>>>To subscribe or unsubscribe via the World Wide Web, visit
>>>    https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>or, via email, send a message with subject or body 'help' to
>>>    r-sig-mixed-models-request@r-project.org
>>>
>>>You can reach the person managing the list at
>>>    r-sig-mixed-models-owner@r-project.org
>>>
>>>When replying, please edit your Subject line so it is more specific
>>>than "Re: Contents of R-sig-mixed-models digest..."
>>>
>>>
>>>Today's Topics:
>>>
>>>   1. how to impose same variance for two random effects with fixed
>>>      correlation patterns? (Gabriel Baud-Bovy)
>>>   2. Re: glmmADMB troubles [diagnosed + worked around] (Ross Boylan)
>>>   3. glmmADMB truncated distributions (Ross Boylan)
>>>   4. Re: glmmADMB truncated distributions (Ben Bolker)
>>>   5. confidence interval for the SD of random effects
>>>      (Iasonas Lamprianou)
>>>   6. variance of fixed effects (Iasonas Lamprianou)
>>>
>>>
>>>----------------------------------------------------------------------
>>>
>>>Message: 1
>>>Date: Thu, 28 Mar 2013 21:14:47 +0100
>>>From: Gabriel Baud-Bovy <baud-bovy.gabriel@hsr.it>
>>>To: r-sig-mixed-models@r-project.org
>>>Subject: [R-sig-ME] how to impose same variance for two random effects
>>>    with fixed correlation patterns?
>>>Message-ID: <5154A4B7.2070702@hsr.it>
>>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>>
>>>Dear all,
>>>
>>>Pinheiro & Bates (1996) give an example of use of pdBlocked to represent
>>>a two-level mixed-effects model as a single-level model (p. 162):
>>>
>>># two-level model
>>>lme(yield ~ nitro, data = Oats,   random = list(Block=pdIdent(~1),
>>>Variety=pdIdent(~1)) )
>>>
>>># single-level model
>>>lme(yield ~ nitro, data = Oats,
>>>      random = list(Block=pdBlocked(list(pdIdent(~1),pdIdent(~Variety-1))))
>>>
>>>This yields a diagonal covariance matrix with variance sigma_1^2 in the
>>>first block
>>>and variance sigma_2^2 in the second block.
>>>
>>>I would like to have the SAME variance in both blocks (this would be a way
>>>to impose that the two variances in the corresponding two-level model
>>>are equal).
>>>
>>>I have looked at pdConstruct.pdBlocked to see whether it was
>>>possible to impose this restriction but the code is daunting and
>>>I would appreciate any help about how to proceed.
>>>
>>>Best,
>>>
>>>Gabriel
>>>
>>>P.S. Let me know if my question (or the previous one about
>>>how to make the covariance matrix for random effects depend
>>>on some characteristics of the grouping factor) is more
>>>appropriate at R-help .
>>>
>>>
>>>
>>>--
>>>---------------------------------------------------------------------
>>>Gabriel Baud-Bovy               tel.: (+39) 02 2643 4839 (office)
>>>UHSR University                       (+39) 02 2643 3429 (laboratory)
>>>via Olgettina, 58                     (+39) 02 2643 4891 (secretary)
>>>20132 Milan, Italy               fax: (+39) 02 2643 4892
>>>
>>>
>>>
>>>------------------------------
>>>
>>>Message: 2
>>>Date: Thu, 28 Mar 2013 15:22:35 -0700
>>>From: Ross Boylan <ross@biostat.ucsf.edu>
>>>To: Ben Bolker <bbolker@gmail.com>
>>>Cc: r-sig-mixed-models@r-project.org
>>>Subject: Re: [R-sig-ME] glmmADMB troubles [diagnosed + worked around]
>>>Message-ID: <5154C2AB.5000502@biostat.ucsf.edu>
>>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>>
>>>I stepped through the code and think I found the problem.  Since I'm
>>>running under emacs/ESS the value of the environment variable SHELL is
>>>funny:
>>>
>>>[context: run_bin does think it's in windows and invokes
>>>  shell(cmd, invisible = TRUE, intern = !verbose)]
>>>
>>>The relevant branch of shell() is (note it uses shell as a variable
>>>inside the function)
>>>     if (missing(shell)) {
>>>         shell <- Sys.getenv("R_SHELL")
>>>         if (!nzchar(shell))
>>>             shell <- Sys.getenv("SHELL")
>>>         if (!nzchar(shell))
>>>             shell <- Sys.getenv("COMSPEC")
>>>R_SHELL is empty
>>>SHELL is  "C:/Program Files/GNU Emacs 24.2/bin/cmdproxy.exe"
>>>I believe that is where the "C:/Program" was coming from in the error
>>>messages.
>>>COMSPEC is  "C:\\Windows\\system32\\cmd.exe".
>>>
>>>I set R_SHELL to the value of COMSPEC and now the command runs.
>>>
>>>Ross
>>>
>>>P.S. shell() in trun invokes system(), which is why you were seeing it
>>>in the traceback.
>>>
>>>On 3/6/2013 6:37 PM, Ross Boylan wrote:
>>>> On 3/2/2013 3:04 PM, Ben Bolker wrote:
>>>>> On 13-02-28 01:15 PM, Ross Boylan wrote:
>>>>>> On 2/27/2013 8:05 PM, Ben Bolker wrote:
>>>>>>> Ross Boylan <ross@...> writes:
>>>>>    [snip: not including version without NAs removed, since we already
>>>>> know what the issue is there]
>>>>>
>>>>>>> r <- glmmadmb(sexActs~(1|id), sexpartner[!is.na(sexpartner$sexActs),],
>>>>>>> debug=TRUE)
>>>>>> platform: windows 32
>>>>>> executable name: glmmadmb.exe
>>>>>> bin_loc:
>>>>>> c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe
>>>>>>
>>>>>>
>>>>>> using temp directory
>>>>>> C:\Users\rdboylan\AppData\Local\Temp\Rtmpy2JsMY\glmmADMB17085c19f4d
>>>>>> creating temp directory
>>>>>> changed working directory to
>>>>>> C:/Users/rdboylan/AppData/Local/Temp/Rtmpy2JsMY/glmmADMB17085c19f4d
>>>>>> Command line:
>>>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>>>
>>>>>> -maxfn 500 -maxph 5 -noinit -shess
>>>>>> Error in system(cmd, intern = intern, wait = wait | intern,
>>>>>> show.output.on.console = wait,  :
>>>>>>    'C:/Program' not found
>>>>>    I'm a little bit baffled here.  What happens if you use save.dir to
>>>>> save the input files to a temporary directory and run
>>>>>
>>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>>
>>>>> -maxfn 500 -maxph 5 -noinit -shess
>>>>>
>>>>> from the command line?
>>>>
>>>>> r <- glmmadmb(sexActs~(1|id),
>>>>> sexpartner[!is.na(sexpartner$sexActs),], debug=TRUE,
>>>>> save.dir="I:/LAMOC/Ross/")
>>>> platform: windows 32
>>>> executable name: glmmadmb.exe
>>>> bin_loc:
>>>> c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe
>>>> changed working directory to I:/LAMOC/Ross
>>>> Command line:
>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>> -maxfn 500 -maxph 5 -noinit -shess
>>>> Error in system(cmd, intern = intern, wait = wait | intern,
>>>> show.output.on.console = wait,  :
>>>>   'C:/Program' not found
>>>>
>>>> Then from a Windows Command Prompt (not cygwin)
>>>> H:\>
>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32
>>>> /glmmadmb.exe" -maxfn 500 -maxph 5 -noinit -shess
>>>> Error trying to open data input file glmmadmb.dat
>>>>  Error trying to read in model data
>>>>  This is usual caused by a missing DAT file
>>>> H:\>I:
>>>>
>>>> I:\>cd LAMOC/Ross
>>>>
>>>> I:\LAMOC\Ross>
>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/b
>>>> /windows32/glmmadmb.exe" -maxfn 500 -maxph 5 -noinit -shess
>>>>
>>>> Initial statistics: 1 variables; iteration 0; function evaluation 0;
>>>> phase 1
>>>> Function value   8.9003579e+04; maximum gradient component mag
>>>> -5.3609e+02
>>>> Var   Value    Gradient   |Var   Value    Gradient   |Var Value
>>>> Gradient
>>>>
>>>>   1  0.00000 -5.36095e+02 |
>>>>
>>>>  - final statistics:
>>>> 1 variables; iteration 8; function evaluation 14
>>>> Function value   5.0576e+04; maximum gradient component mag 6.2393e-08
>>>> Exit code = 1;  converg criter   1.0000e-04
>>>> Var   Value    Gradient   |Var   Value    Gradient   |Var Value
>>>> Gradient
>>>>
>>>>   1 112.5711  6.23928e-08 |
>>>> etc
>>>>
>>>> So it seems to work, provided I start in the save.dir. Note that is
>>>> the directory that R is running in.
>>>> glmmadmb.exe is still running as I hit send.
>>>>
>>>>>
>>>>>   What is the result of .Platform (and .Platform$OS in particular)
>>>>
>>>>> .Platform
>>>> $OS.type
>>>> [1] "windows"
>>>>
>>>> $file.sep
>>>> [1] "/"
>>>>
>>>> $dynlib.ext
>>>> [1] ".dll"
>>>>
>>>> $GUI
>>>> [1] "RTerm"
>>>>
>>>> $endian
>>>> [1] "little"
>>>>
>>>> $pkgType
>>>> [1] "win.binary"
>>>>
>>>> $path.sep
>>>> [1] ";"
>>>>
>>>> $r_arch
>>>> [1] "i386"
>>>>
>>>> Ross
>>>>
>>>>>
>>>>>   It looks conceivably like R is misdiagnosing your system as *not*
>>>>> being
>>>>> windows, as that's the only way system() should be running.  Are you
>>>>> running under Cygwin (you say it's installed below) ...  ?
>>>>>
>>>>>> changed working directory to i:/LAMOC/Ross
>>>>>> removed temp directory
>>>>>> C:\Users\rdboylan\AppData\Local\Temp\Rtmpy2JsMY\glmmADMB17085c19f4d
>>>>>>> glmmADMB:::get_bin_loc()
>>>>>> $bin_loc
>>>>>> [1]
>>>>>> "c:/Users/rdboylan/Documents/R/R-2.15.2/site-library/glmmADMB/bin/windows32/glmmadmb.exe"
>>>>>>
>>>>>>
>>>>>>
>>>>>> $platform
>>>>>> [1] "windows"
>>>>>>
>>>>>> P.S. about lme4; I don't have a build environment and so trying the
>>>>>> github version will not be my first move.
>>>>>> Although perhaps lack of a build environment is why the second version
>>>>>> is failing.  I do have cygwin installed, althoughI would not expect
>>>>>> R to
>>>>>> know how to find it.
>>>>>    lme4 should be installable from lme4.r-forge.r-project.org/repos now,
>>>>> as stated in a message earlier today.
>>>>>
>>>>
>>>
>>>
>>>
>>>------------------------------
>>>
>>>Message: 3
>>>Date: Thu, 28 Mar 2013 15:39:07 -0700
>>>From: Ross Boylan <ross@biostat.ucsf.edu>
>>>To: r-sig-mixed-models@r-project.org
>>>Subject: [R-sig-ME] glmmADMB truncated distributions
>>>Message-ID: <5154C68B.5090807@biostat.ucsf.edu>
>>>Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>>
>>>How are they truncated?  I expected there might be a way to set the
>>>truncation threshhold, but I don't see it in the docs.
>>>Ross Boylan
>>>
>>>
>>>
>>>------------------------------
>>>
>>>Message: 4
>>>Date: Fri, 29 Mar 2013 03:39:53 +0000 (UTC)
>>>From: Ben Bolker <bbolker@gmail.com>
>>>To: r-sig-mixed-models@r-project.org
>>>Subject: Re: [R-sig-ME] glmmADMB truncated distributions
>>>Message-ID: <loom.20130329T043015-425@post.gmane.org>
>>>Content-Type: text/plain; charset=us-ascii
>>>
>>>Ross Boylan <ross@...> writes:
>>>
>>>>
>>>> How are they truncated?  I expected there might be a way to set the
>>>> truncation threshhold, but I don't see it in the docs.
>>>> Ross Boylan
>>>>
>>>
>>>  They're truncated at zero.
>>>
>>>  One of the problems with developing a flexible system is that everyone
>>>comes up with new uses for it: the original intent of the truncated
>>>distributions was just to make hurdle models possible.
>>>
>>>  It would be possible but not completely trivial to adapt glmmADMB
>>>to handle other truncation points -- you would need a straightforward
>>>way to compute the cumulative distribution function of the response up
>>>to the truncation point. Probably the best thing for people who want
>>>to extend glmmADMB in these various ways is to learn how to use the
>>>full version of AD Model Builder ...
>>>
>>>  Ben Bolker
>>>
>>>
>>>
>>>------------------------------
>>>
>>>Message: 5
>>>Date: Thu, 28 Mar 2013 23:19:24 -0700 (PDT)
>>>From: Iasonas Lamprianou <lamprianou@yahoo.com>
>>>To: "r-sig-mixed-models@r-project.org"
>>>    <r-sig-mixed-models@r-project.org>
>>>Subject: [R-sig-ME] confidence interval for the SD of random effects
>>>Message-ID:
>>>    <1364537964.28361.YahooMailNeo@web160104.mail.bf1.yahoo.com>
>>>Content-Type: text/plain
>>>
>>>
>>>
>>>Dear all, I am running a very simple mixed effects model using lmer. I get the following results:
>>>
>>>Random effects: Groups   Name        Variance Std.Dev. id       (Intercept) 7.4223   2.7244   Residual             1.9767   1.4060
>>>Number of obs: 8080, groups: id, 2020
>>>
>>>
>>>Apparently, a very good question is whethe0r the estimate of the parameters of the random effect are indeed bigger than zero. So, I used the bootstrap method described by Faraway where the distribution of likelihood rations is built between a model with and withot the random effect. I get a probbility p<0.0001, so the random effect estimates are significant. I also confirmed these findings using the library("RLRsim").
>>>
>>>
>>>So far so good, but one of the reviewers asked me to give 95% confidence intervals. I spent three hours on the internet searching for good solutions and decided that I should use the following:
>>>
>>>sm1 <- mcmcsamp(mlm.1.1g, n=1000, saveb = TRUE)
>>>lme4::HPDinterval(sm1)
>>>
>>>
>>>However, the 95% lower and upper bounds do not indclude the estimate!!!! See below:
>>>
>>>> r$sigma  lower    upper
>>>[1,] 1.977582 2.058084
>>>attr(,"Probability")
>>>[1] 0.95
>>>
>>>How is this possible?
>>>
>>>Also, I assume that the ST shows the ratio of [random effect variance]/[total variance]. In this case, my ICC is 7.42/(7.42+1.97)=0.79 but see the foloowing:
>>>
>>>> r$ST  lower     upper
>>>[1,] 0.6905202 0.7342866
>>>attr(,"Probability")
>>>[1] 0.95
>>>My empirical ICC is not included in the confidence interval. Apparently, I am missing something but after spending a lot of time, I need to ask you offer me some of your experience.
>>>Thank you for your time
>>>
>>>Jason
>>>
>>>
>>>
>>>
>>>Dr. Iasonas Lamprianou
>>>Department of Social and Political Sciences
>>>University of Cyprus
>>>
>>>
>>>
>>>
>>>>
>>>    [[alternative HTML version deleted]]
>>>
>>>
>>>
>>>------------------------------
>>>
>>>Message: 6
>>>Date: Thu, 28 Mar 2013 23:23:05 -0700 (PDT)
>>>From: Iasonas Lamprianou <lamprianou@yahoo.com>
>>>To: "r-sig-mixed-models@r-project.org"
>>>    <r-sig-mixed-models@r-project.org>
>>>Subject: [R-sig-ME] variance of fixed effects
>>>Message-ID:
>>>    <1364538185.83291.YahooMailNeo@web160102.mail.bf1.yahoo.com>
>>>Content-Type: text/plain
>>>
>>>
>>>Note: I am sending this again because I got an email saying that it was not delivered. In case it was delivered initially, I apologise for the double-posting.
>>>
>>>
>>>
>>>Dear all, I am running a very simple mixed effects model using lmer. I get the following results:
>>>
>>>Random
>>>effects: Groups   Name        Variance Std.Dev. id       (Intercept)
>>>7.4223   2.7244   Residual             1.9767   1.4060
>>>Number of obs: 8080, groups: id, 2020
>>>
>>>
>>>Apparently,
>>>a very good question is whethe0r the estimate of the parameters of the
>>>random effect are indeed bigger than zero. So, I used the bootstrap
>>>method described by Faraway where the distribution of likelihood rations
>>>is built between a model with and withot the random effect. I get a
>>>probbility p<0.0001, so the random effect estimates are significant. I
>>>also confirmed these findings using the library("RLRsim").
>>>
>>>
>>>So
>>>far so good, but one of the reviewers asked me to give 95% confidence
>>>intervals. I spent three hours on the internet searching for good
>>>solutions and decided that I should use the following:
>>>
>>>sm1 <- mcmcsamp(mlm.1.1g, n=1000, saveb = TRUE)
>>>lme4::HPDinterval(sm1)
>>>
>>>
>>>However, the 95% lower and upper bounds do not indclude the estimate!!!! See below:
>>>
>>>> r$sigma  lower    upper
>>>[1,] 1.977582 2.058084
>>>attr(,"Probability")
>>>[1] 0.95
>>>
>>>How is this possible?
>>>
>>>Also,
>>>I assume that the ST shows the ratio of [random effect variance]/[total
>>>variance]. In this case, my ICC is 7.42/(7.42+1.97)=0.79 but see the
>>>foloowing:
>>>
>>>> r$ST  lower     upper
>>>[1,] 0.6905202 0.7342866
>>>attr(,"Probability")
>>>[1] 0.95
>>>My
>>>empirical ICC is not included in the confidence interval. Apparently, I
>>>am missing something but after spending a lot of time, I need to ask
>>>you offer me some of your experience.
>>>Thank you for your time
>>>
>>>Jason
>>>
>>>
>>>
>>>Dr. Iasonas Lamprianou
>>>Department of Social and Political Sciences
>>>University of Cyprus
>>>    [[alternative HTML version deleted]]
>>>
>>>
>>>
>>>------------------------------
>>>
>>>_______________________________________________
>>>R-sig-mixed-models mailing list
>>>R-sig-mixed-models@r-project.org
>>>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>>
>>>End of R-sig-mixed-models Digest, Vol 75, Issue 45
>>>**************************************************
>>>
>>>
>>>
>>         [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
>
>
>-- 
>Joshua Wiley
>Ph.D. Student, Health Psychology
>University of California, Los Angeles
>http://joshuawiley.com/
>Senior Analyst - Elkhart Group Ltd.
>http://elkhartgroup.com
>
>
>
	[[alternative HTML version deleted]]

