RE: [R] Repeated measures
If you fit the model as suggested below and store it in an object, say model.lme, then you can update the model as follows to include the continuous AR1 structure model2.lme-update(model.lme, correlation=corCAR1(form=~time|cow) ) To compare the two models, you might use the LRT as: anova(model.lme,model2.lme) Harold -Original Message- From: [EMAIL PROTECTED] on behalf of Spencer Graves Sent: Sun 6/6/2004 12:21 PM To: Alex Bach Cc: [EMAIL PROTECTED] Subject: Re: [R] Repeated measures 1. You didn't say which manual you were reading on lme, if you have not consulted Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer), I suggest you do so. The issues are discussed in greater depth with many examples. I found this book well worth the time and money I invested in it. 2. Have you considered the following: lme(milk yield = cow + treatment + time + treatment*time, random = ~time|cow ... ) This will NOT have the ARCH(1) R RCORR error structure; lme can handle certain types of autocorrelated error structures, but I don't remember the details at the moment. Pinheiro and Bates discuss some capabilities of this nature. hope this helps. spencer graves Alex Bach wrote: Dear R-gurus, I am pretty much new on R. I am trying to to do a repeated analysis of a linear mixed model with R, and I consistently fail... The problem is: Cow is the random factor, treatment is the fixed factor. The dependent variable is milk yield, which is measured several times (repeatedly over time), thus there is another variable which is time (i.e. week). The model would be something like this: milk yield = cow + treatment + time + treatment*time With time as a repeated measure. Would some one be kind enough to guide on how could I set up the statement in R? I imagine I have to use LME but I have not been smart enough to figure out how to do it by just reading its manual. Also, with SAS there is a nice contrast, called SLICE that allows you to test when in time there is difference between treatments. I do not know if there is something like this is R. Thank you very much! PD. For SAS users, what I am using in SAS to perform this analysis (with an Autoregressive covariance structure) the program would read like this: proc MIXED covtest; CLASS cow treat time; MODEL yield= treat time treat*time; REPEATED time/SUB=cow(treat) TYPE=ARH(1) R RCORR; RANDOM cow; LSMEANS treat time treat*time/SLICE = time; RUN; Thank you very much, Sincerely, Alex __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] signifikanz?
You need cor.test(x,y). However, I do not think you mean significant at the .95 or .99 level, do you? -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Martin Klaffenboeck Sent: Tuesday, June 01, 2004 9:15 AM To: r-help Subject: [R] signifikanz? Hello, When I use: cor(x, y) I get a correlations coefficient. But how can I see if it is signifikant on a .95 or .99 level? Thanks, Martin __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Merging nlme output
Dear list: I am trying to merge two files together from output I get based on the coef() command. Here is what I am running into. I have two simple linear mixed models mod1.lme-lme(math~year, data=sample, random=~year|childid/schoolid) mod2.lme-lme(math~year, data=sample, random=~year|childid) I then call the coefficients and store them in the following objects using mod1.coef-coef(mod1.lme, level=2) mod2.coef-coef(mod2.lme) The problem is that the IDs are different for the two dataframes and I cannot seem to figure out how to make the merge happen nicely. For example, this is output from mod1.lme for the first 5 students mod1.coef[1:5,] (Intercept) year 2020/273026452 0.03923346 0.9575926 2020/273030991 0.91772318 1.0773731 2020/273059461 0.43865139 1.0111789 2020/278058841 1.11292376 1.0526057 2020/292017571 1.09661340 1.0774692 As opposed to out for the first 5 cases in mod2.coef mod2.coef (Intercept) year 101480302-1.13483082 0.7023644 173559292 0.52818576 0.8639216 174743401-1.18038537 0.6824612 174755092 0.04760782 0.7402835 179986251-2.42848405 0.5677182 You can see that the ID variable for the first output includes the grouping variable schoolid followed by a / and then the childid. I am interested in exploring the differences in the growth trajectories between the two models and would like to merge the two datasets to be able to do so. What do I need to do to somehow eliminate the schoolid and / from mod1.lme so that I can use the merge() command? Thanks, Harold C. Doran Director of Assessment One Massachusetts Avenue, NW · Suite 700 Washington, DC 20001-1431 202.336.7075 [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Help with Plotting Function
Dear List: I cannot seem to find a way to plot my data correctly. I have a small data frame with 6 total variables (x_1 ... x_6). I am trying to plot x_1 against x_2 and x_3. I have tried plot(x_2, x_1) #obviously works fine plot(x_3, x_1, add=TRUE) # Does not work. I keep getting error messages. I would also like to add ablines to this plot. I have experimented with a number of other plotting functions and I cannot seem to get this to work. The data are student achievement data. I am trying to plot percentile ranks against scale scores for different grade levels. When plotted as such, they look like logistic curves. I am trying to show graphically the distance along x a student must grow simply to remain at the same percentile rank, y. Thanks. Harold C. Doran One Massachusetts Avenue, NW · Suite 700 Washington, DC 20001-1431 202.336.7075 [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Fatal Error
Dear List: When trying to open 1.9.0 this morning, I have the following error: Fatal Error: Unable to restore saved data in .Rdata I am using Windows 2000. The program then quits. Do I need to reinstall? Harold C. Doran One Massachusetts Avenue, NW · Suite 700 Washington, DC 20001-1431 202.336.7075 [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Fatal Error
Thank you. Locating and deleting the bad .Rdata file did the trick and I can now work in R. Thanks. -Original Message- From: Duncan Murdoch [mailto:[EMAIL PROTECTED] Sent: Monday, May 17, 2004 8:45 AM To: Harold Doran Cc: [EMAIL PROTECTED] Subject: Re: [R] Fatal Error On Mon, 17 May 2004 08:30:57 -0400, Harold Doran [EMAIL PROTECTED] wrote: Dear List: When trying to open 1.9.0 this morning, I have the following error: Fatal Error: Unable to restore saved data in .Rdata I am using Windows 2000. The program then quits. Do I need to reinstall? It should be sufficient to delete the bad .Rdata file. It's normally stored in the directory that's current when R starts; if you're using Windows, you can look in the shortcut to find what the starting directory is. If you can make this error reproducible, we'd like to hear about it. Duncan Murdoch __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] 2 lme questions
There are two way to accomplish this in nlme. First try using the summary() command, which will produce all variance components and estimates for the fixed effects. Also, try the following to extract the point estimates and approximate CIs for the variance comonents. intervals(model.lme, which=var) Harold -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Steve Roberts Sent: Monday, April 05, 2004 3:32 PM To: [EMAIL PROTECTED] Cc: Steve Roberts Subject: [R] 2 lme questions Greetings, 1) Is there a nice way of extracting the variance estimates from an lme fit? They don't seem to be part of the lme object. 2) In a series of simulations, I am finding that with ML fitting one of my random effect variances is sometimes being estimated as essentially zero with massive CI instead of the finite value it should have, whilst using REML I get the expected value. I guess it is a numerical/optimisation problem but don't know enough about the lme fitting algorithm to know which tollerance/scale parameter to mess about with. Any suggestions where to start? Thanks, Steve. [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] general mixed model statement
This would be accomplished via the nlme library for mixed linear models. See Pinhiero and Bates (2000) for all relevant documentation and data applications. Harold -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Sent: Thursday, April 01, 2004 10:21 AM To: [EMAIL PROTECTED] Cc: [EMAIL PROTECTED] Subject: [R] general mixed model statement I am trying to analyze a general mixed model, but am confused about the model statement to use. The structure of the problem is the following: y = X b + Z u + e where X and Z are known incidence matrices for fixed and random effects respectively, b is a vector of fixed effects to be estimated, u is a vector of random effects, and e is a vector of error terms. Additionally, the covariance matrices of the random effects and error terms are given by G = g A, and R = r E where A and E are known matrices and g and r must be estimated. Presumably there is a way to specify in a model statement the known covariance structure for the random and error terms, but I cannot find the documentation I need. Please either clarify for me how to write this out or point me to the appropriate documentation that outlines the full range of model statements that are possible. Thank you very much for your help. Cheers, Brook __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] rate of change
Type help(deriv). However, it may be easier to compute the derivative by hand for a simple expression. HCD -Original Message- From: Fred J. [mailto:[EMAIL PROTECTED] Sent: Monday, March 15, 2004 10:20 PM To: r help Subject: [R] rate of change Hello I am wondering, how do I find if R has a certain funciton to do a given task. do I just type help.search(rate). I am just trying to find a function to calculate the rate of change for a variable. I could come up with one if there isn't any allready builtin. thanks __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] still spss
It appears you may have left off the file extension, experiencial.sav -Original Message- From: Margarida Júlia Rodrigues Igreja [mailto:[EMAIL PROTECTED] Sent: Tuesday, March 04, 1997 2:24 AM To: [EMAIL PROTECTED] Subject: [R] still spss hi again, i still cannot open the file in spss :( i type: library(foreign) read.spss(H:\\Desktop\\bd1\\experiencia1) and the error comes: Error in read.spss(H:\\Desktop\\bd1\\experiencia1) : unable to open file can you help me? margarida,portugal __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] how to delete a matrix column
You can use dataset$columnname=NULL HCD -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 -Original Message- From: Aimin Yan [mailto:[EMAIL PROTECTED] Sent: Tuesday, March 02, 2004 11:59 AM To: [EMAIL PROTECTED] Subject: [R] how to delete a matrix column Hello, I am new to R, How to delete a matrix column. Thanks, __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Area between CDFs
Thanks. I have been able to create the following simple function to examine the
vertical gap between two CDFs at a value along the x-axis that I specify. For example,
I create the ECDFs:
male.ecdf-ecdf(egmale$math)
female.ecdf-ecdf(egfemale$math)
I then define the following function:
dif.cdf-function(x){return(abs(female.ecdf(x)-male.ecdf(x)))}
Now, I can use the function to measure the gap at values along the x-axis (i.e., gap =
F(x)-G(x). Also, the CDFs do not cross at any point):
dif.cdf(0)
which returns a value that is the size of the gap at a specific score between males
and females.
What I would like to be able to do is measure the vertical gap at each point along the
x-axis and then plot the gap. This would illustrate for how large differences in
student achievement are at different score values into a nice visual display.
The brute force way seems to use the function above for each score value. However,
this is, of course, inefficient. Any ideas on how I might be able to create a function
that would be more efficient?
Many thanks,
Harold
--
Harold C. Doran
Director of Research and Evaluation
New American Schools
675 N. Washington Street, Suite 220
Alexandria, Virginia 22314
703.647.1628
-Original Message-
From: Thomas Lumley [mailto:[EMAIL PROTECTED]
Sent: Thursday, February 19, 2004 10:49 AM
To: Samuelson, Frank*
Cc: [EMAIL PROTECTED]
Subject: RE: [R] Area between CDFs
On Wed, 18 Feb 2004, Samuelson, Frank* wrote:
You may not want to integrate cdfs. They're already probabilities. :)
Nice analytic statistics exist for just the maximum distance between
the cdfs, for example.
And for the area between cdfs, which is perhaps better known as the
difference in means.
-thomas
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Area between CDFs
Dear List: I am trying to find the area between two ECDFs. I am examining the gap in performance between two groups, males and females on a student achievement test in math, which is a continuous metric. I start by creating a subset of the dataframe male-subset(datafile, female=Male) female-subset(datafile, female=Female) I then plot the two CDFs via plot.ecdf(male$math) plot.ecdf(female$math, add=TRUE) This produces the visual display that reveals a gap in performance. What I would like to do is learn to perform the integration between the two ECDFs to examine the size of this gap. I would also like to try and examine the horizontal distance between the two CDFs via another visual display. In other words, the distance between, say the 50th percentile, from each CDF (or, the distance along the x-axis from cdf1 to cdf2 at each percentile. Ideally, I would like to plot this horizontal gap at each percentile. Secondly, I would like to try and measure and plot the vertical gap, i.e., the distance along the y-axis from cdf 1 to cdf2 at each value along the x-axis. I am not sure if I first need to smooth the ECDFs before performing these operations. Any help would be appreciated. I hope this makes sense. Harold -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Sweave/LaTeX Problem with EPS PDF
Yes, they were lattice and your suggestion did the trick. Many thanks! -Original Message- From: Jason Turner [mailto:[EMAIL PROTECTED] Sent: Sun 12/21/2003 12:30 AM To: Harold Doran Cc: [EMAIL PROTECTED] Subject: Re: [R] Sweave/LaTeX Problem with EPS PDF Harold Doran wrote: [... Sweave use...] However, when I examine the pdf or eps files created, there is nothing there. When I view the EPS using Ghostview, the file is empty, but there appears to be a bounding box surrounding nothing. When I open the pdf graphic, there is nothing there either. I have tried creating both dvis and pdf files. Again, text works perfectly, but graphics do not work. Are they lattice graphics? They need to be wrapped in a print() statement, like print(xyplot(...)) Without that, they don't produce output to eps or pdf. Cheers Jason -- Indigo Industrial Controls Ltd. http://www.indigoindustrial.co.nz 64-21-343-545 [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Sweave/LaTeX Problem with EPS PDF
Dear List:
I am unsure if my problem is with Sweave or LaTeX. Anyhow, I am using the MikTeX
distribution and TexnicCenter.
I can easily create Sweave files and all goes well until I try to incorporate
graphics. I use the same code as found in the examples found in the users manual.
In R, the graphics I want are created as Sweave is creating the .tex file. When I
examine the .tex file created by Sweave, it includes the includegraphics{} statement
needed.
When I run LaTeX on the .tex file, everything works except that the graphics I want
are not displayed.
However, when I examine the pdf or eps files created, there is nothing there. When I
view the EPS using Ghostview, the file is empty, but there appears to be a bounding
box surrounding nothing. When I open the pdf graphic, there is nothing there either. I
have tried creating both dvis and pdf files. Again, text works perfectly, but graphics
do not work.
Does anyone have any suggestions?
Many thanks,
Harold
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] varFixed
Dear List: Earlier this week I posted a question and received no response, and I continue to struggle with my model. My original question is pasted below. I am using lme and want to fix the variance of the within group residual at 1 (e~n(0,1). I think the varFixed function should be used to accomplish this, but I am struggling to figure out how to do this. Can anyone offer suggestions on how this might be accomplished? Thanks, I would appreciate any suggestions. Harold Dear List: I am trying to figure out how to incorporate measurement error in an longitudinal educational data set using lme to create a true score model. As a by-product of the procedures used to scale educational tests, one can obtain a person-specific measurement error associated with each score, or a conditional standard error. For example, a score of 200 would have measurement error specific to that score that would be different than, say, a score of 250. I have been rather successful in figuring out how to rescale the necessary components to create this true score model. This simply requires that the response variable, the intercept, and any other variables in the design matrix be multiplied by the reciprocal of the standard error of measurement for the associated score. There may be a better way to do this, but I manually create a vector of 1s for all observations and multiply this vector by 1/sem. This is the new intercept. I also multiply any other predictors in the design matrix by the same value. In the R code, I remove the intercept included by default (-1) and include the newly created intercept (which is no longer a constant) as well as the new response variable and rescaled predictors. However, I am confused regarding the within-group error term. Fitting this model requires that the variance be fixed at 1: e ~ n(0,1). Is it possible to constrain the variance for this model as such? I would appreciate any comments or suggestions regarding this model. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] error constraints in lme
Dear List: I am trying to figure out how to incorporate measurement error in an longitudinal educational data set using lme to create a true score model. As a by-product of the procedures used to scale educational tests, one can obtain a person-specific measurement error associated with each score, or a conditional standard error. For example, a score of 200 would have measurement error specific to that score that would be different than, say, a score of 250. I have been rather successful in figuring out how to rescale the necessary components to create this true score model. This simply requires that the response variable, the intercept, and any other variables in the design matrix be multiplied by the reciprocal of the standard error of measurement for the associated score. There may be a better way to do this, but I manually create a vector of 1s for all observations and multiply this vector by 1/sem. This is the new intercept. I also multiply any other predictors in the design matrix by the same value. In the R code, I remove the intercept included by default (-1) and include the newly created intercept (which is no longer a constant) as well as the new response variable and rescaled predictors. However, I am confused regarding the within-group error term. Fitting this model requires that the variance be fixed at 1: e ~ n(0,1). Is it possible to constrain the variance for this model as such? I would appreciate any comments or suggestions regarding this model. HCD -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] conf int mixed effects
I am very curious about this. If a particular growth model is specified to reflect
repeated observations on individual i in unit j, such as:
y_{tij} = [B_{00} + B_{01}*(TIME)]+[u_{00}+u_{01}*(TIME)+ e_{tij}]
where Bs are the fixed effects and the u's are the random effects.
The growth of individual i is then computed as:
B_{01} + u_{01}
Is it appropriate to compute a confidence interval around this growth rate? I so, how
might this be accomplished? Based on Doug's comments below, it would seem that only a
CI can be formulated for the fixed portion of the model.
I would appreciate any clarification.
HCD
--
Harold C. Doran
Director of Research and Evaluation
New American Schools
675 N. Washington Street, Suite 220
Alexandria, Virginia 22314
703.647.1628
-Original Message-
From: Douglas Bates [mailto:[EMAIL PROTECTED]
Sent: Thursday, November 13, 2003 10:11 AM
To: Joerg Schaber
Cc: [EMAIL PROTECTED]
Subject: Re: [R] conf int mixed effects
Joerg Schaber [EMAIL PROTECTED] writes:
I have a linear mixed-effects model object and want to extract the 95%
confidence intervals for the fixed and random effects, respectively. I
found the function intervals() for confidence intervals for the fixed
effects but no corresponding function for the random effects. Does it
exist or do I have to calculate the confidence intervals for the
random effects myself?
You have to calculate them yourself, partly because it is not clear
what such an interval should be. Technically, the random effects are
not parameters and defining a confidence interval on a random
variable that is part of the model is, at the very least, awkward.
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Selecting a random sample for lmList()
Dear List: I have a data set with over 7000 students with about 4 observations over time per student. I want to examine the within-group fits of a random sample of this group as it takes forever to compute and draw all 7000 regressions. Here is the code I have used so far. group-groupedData(math~year|childid, data=scores) group.list-lmList(group) plot(augPred(group.list)) How might I select a random sample of say 100 students so that I can visually examine their trajectories? Thank you -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Doubly Multivariate LME
Dear R: I am trying to fit a doubly multivariate LME (DM) where I have two response variables measured on two occasions per person. Specifically, reading and math scores measured at the beginning and ending of a school year. The response variables have a correlation of r = .85. The response variables in the data matrix are stacked in a vector with a dummy code flagging each outcome and with time variables for each outcome. The model was fit by removing the overall intercept, but creating fixed effects and random effects for each using the following: mult2.lme-lme(fixed=score~-1+read+math+time.m+time.r, data=mult.samp, random=~-1+read+math+time.r+time.m|childid) This worked and seemed to produce reasonable estimates. I then ran a model using only a single outcome (reading) and found that the estimates are very similar, so I am relatively confident in the results of the model. Now, I have a couple of questions regarding the DM LME: 1) If I wanted to explore model assumptions, i.e., homoscedasticity and dependence among the residuals, how might I do this. For example, how might I specify an AR1 structure? I have explored the assumptions in the single response model, but am having trouble now. 2) I presume it is safe to explore the intercepts and slopes using lmList () one variable at a time. Is this correct? 3) The AIC statistic for the single response model is half the size of the DM model. It is my understanding that this statistic is more appropriate than LL test in this scenario. Is this correct? If so, is there a reason that may explain the larger AIC in the DM model? Or, is this indicating that the single response model is a better fit? 4) Last, is there any other issue that I am missing? I may not have asked the question, but would appreciate any suggestions related to the model. Regards, -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Downloading LME4?
Dear R: Am I having trouble downloading the LME4 library. I am using Windows and am using ver 1.7 I have tried the following: 1) Install package from CRAN, but LME4 is not listed 2) Downloaded LME4 from http://cran.us.r-project.org/, however, I cannot open the file when I try install from local drive. I get the following error: Error in file(file, r) : unable to open connection In addition: Warning messages: 1: error 1 in extracting from zip file 2: cannot open file `lme4_0.3-7.tar.gz/DESCRIPTION' Any help is appreciated. -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Off Topic: Good reference for sample size calculations
Jacob Cohen's book Statistical Power Analysis for the Behavioral Sciences is one.
--
Harold C. Doran
Director of Research and Evaluation
New American Schools
675 N. Washington Street, Suite 220
Alexandria, Virginia 22314
703.647.1628
http://www.edperform.net
-Original Message-
From: Warnes, Gregory R [mailto:[EMAIL PROTECTED]
Sent: Wednesday, September 10, 2003 8:37 AM
To: '[EMAIL PROTECTED]'
Subject: [R] Off Topic: Good reference for sample size calculations
Hi All,
This is off topic, but we're drawing a blank here..
In a presentation I'll be giving next week, I want to include a reference
to a good general text on computing sample sizes for standard experiments.
Can anyone recommend a good book to use for this purpose?
Thanks,
-Greg
LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}}
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Levene test of homogeneity of variance
I believe it is in the Rcmdr package, which requires the car library to be loaded. You can also perform an ANOVA using the absolute value of the deviations from each respective group mean, which is what Levene's Test does. -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net -Original Message- From: Haynes, Maurice (NIH/NICHD) [mailto:[EMAIL PROTECTED] Sent: Wednesday, August 13, 2003 10:40 AM To: 'r-help' Subject: [R] Levene test of homogeneity of variance Has the Levene test of homogeneity of variance been implemented in any library in R? Thanks, Maurice Haynes National Institute of Child Health and Human Development Child and Family Research Section 6705 Rockledge Drive Bethesda, MD 20892 Voice: 301-496-8180 Fax: 301-496-2766 E-Mail: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Conditional Statements for Graphing
Dear List I have math test scores for male and female students where gender is a dummy code (female =1). I also have a variety of other demographic variables. However to begin, I want to create a very simple stripchart where female math scores are a blue circle and male scores are a red triangle. I am having difficulty using conditional statements to accomplish this. Thank you. -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] how to test whether two slopes are sign. different?
Dear Stoet This can be handled well in using a mixed-effects model, library (nmle). You can use the lmList option to check whether the slopes differ across populations. -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net -Original Message- From: Gijsbert Stoet [mailto:[EMAIL PROTECTED] Sent: Sunday, July 20, 2003 10:51 PM To: [EMAIL PROTECTED] Subject: [R] how to test whether two slopes are sign. different? Hi, suppose I do want to test whether the slopes (e.g. determined with lsfit) of two different population are significantly different, how do I test this (in R). Say for example, I found out what the slope between age and number of books read per year is for two different populations of subjects (e.g. 25 man and 25 woman), say using lsfit. How can I tell whether the slopes are different in R. (And how would I do it for regression coefficients?) Thanks a lot for your help. __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] Maximum Likelihood Estimation and Optimisation
Well, lm() produces an OLS solution, which are also MLE solutions for the fixed effects. I think this is an easy way, although maybe not the best. BHHH is a numerical approximation that can be used when a closed form solution is not available. It is less sophisticated than Newton-Raphson. Is this helpful? -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 -Original Message- From: Fohr, Marc [AM] [mailto:[EMAIL PROTECTED] Sent: Thursday, July 10, 2003 10:17 AM To: [EMAIL PROTECTED] Subject: [R] Maximum Likelihood Estimation and Optimisation Hello, I want to calculate a maximum likelihood funktion in R in order to solve for the parameters of an estimator. Is there an easy way to do this in R? How do I get the parameters and the value of the maximum likelihood funktion. More, I want to specify the algorithm of the optimisation above: BHHH (Berndt Hall Hall Hausman). Is this possible? Thanks a lot for your help and best regards Marc - Marc Fohr, CFA Equity Portfolio Manager First Private Investment Management Neue Mainzer Strasse 75 D-60311 Frankfurt/Main Phone: ++49 - 69 - 2607 5424 Fax: ++49 - 69 - 2607 5440 Email: [EMAIL PROTECTED] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
RE: [R] NLME Fitted Values
I would like to be able to add the random effect for the slope to the fitted value for the slope at each level (i.e, for each school and for each student). The fitted values at the population level should be the fixed effects (B_00 and B_01). So each school should have the same fitted values, but a unique random effect. Adding the two together results in the mean adjusted growth rate for each school. But when I use fitted() and set level=0, I get results that are not the fixed effects. I am unsure what the fitted values are that are produced. At the child level, each child has a unique random effect, but each child within the same school should share the same fitted value. But the fitted values for children at different schools should be different. Adding these should result in the mean adjusted growth rate for each child. Although I get unique random effects for each child using ranef(), I do not get the fitted values that I think I should get wth fitted(). However, when I use level.1-data.frame(coef(eg.model1), level=1) this gives me what HLM calls the fitted values for each child. So, I think I can add these to the random effects at the same level to get the mean growth rate for each child. So, I am a little unclear why fitted() is not producing results that are similar to HLM fitted values, but coef() produces the values that HLM calls the level 2 fitted values. So, the R and HLM random effects are convergent at all levels. The HLM fitted values at the population level are the fixed effects, but this is not the case when I use fitted () at level =0 in R. However, the fitted values at the observation level in R are the same fitted values at the observation level in HLM. In sum, HLM and R random effects are exactly the same at all levels. HLM and R fitted values are exactly the same at the level of observation. But, the fitted vaues at other levels are very different. Am I using fitted incorrectly and should instead be using coef() to accomplish my goal? -Original Message- From: Douglas Bates [mailto:[EMAIL PROTECTED] Sent: Tue 7/8/2003 5:52 PM To: Harold Doran Cc: [EMAIL PROTECTED] Subject: Re: [R] NLME Fitted Values Harold Doran [EMAIL PROTECTED] writes: Dear List: I am having difficulties with the fitted values at different levels of a multilevel model. My data set is a series of student test scores over time with a total of 7,280 observations, 1,720 students nested witin 60 schools. The data set is not balanced. The model was fit using eg.model.1-lme(math~year, random=~year|schoolid/childid, data=single). When I call the random effects at all levels using EB.1-data.frame(ranef(eg.model1, level=1)) and EB.2-data.frame(ranef(eg.model1, level=2)), I get the shrinkage estimators that I expect. That is, I get 2 random effects for each child (1 intercept and 1 slope) and 2 for each school (1 intercept and 1 slope). However, when I call the fitted values using: fitted-data.frame(fitted(eg.model1, level=0:2)), I get 7,280 fitted values at the level of observation. This makes sense (one for each observed score). However, I also get 7,280 fitted values at the child and at the school level. This does not seem correct to me. I should only have, I think, 60 fitted values at the school level (actually, 1 intercept and 1 slope for each of 60 schools) and 1,720 fitted values at the child level (again, 1 intercept and one for the slope for each child). I think you are confusing the random effects with the fitted values. The fitted values will depend on the fixed-effects and, when level 0, on the random effects. Because the year term, which is associated with the fixed effects, can change within school and within child we always return one fitted value for each observation. Why am I always getting 7,280 fitted values? There is some redundancy but we cannot determine the redundancy without knowing the exact form of both the fixed effects and the random effects. In your case where year is the only term in the fixed effects when level = 0 the fitted values for the same year should be the same. When level = 1 the fitted values for the same year and same school but different children should be the same. When level = 2 potentially all the fitted values will be different. However, it could be that a fixed effects term would have a unique value at each row and then even the level = 0 fitted values could all be distinct. I have tried
RE: [R] within group variance of the coeficients in LME
I think using the metrics provided (standard deviation units) one can assess the degree of variability and determine whether differences are large enough to be considered practically significant. I suppose Dr. Bates and others can comment further, but making decisions solely on the basis of statistical significance is not always the most reasonable method for making research decisions, especially when it is largely a function of sample size. In my HLM to R experience, I have found the lmList function to provide the best answer to your question, although it is not a quantitative test like the chi-square in HLM. You can fit a linear model for each subject, use the intervals function, and plot the intercepts and the slopes for all units with 95% CIs. You can visually see whether the intercepts and slopes vary across units. I think this is a very strong method for examining variability. However, you can eliminate the random slope (or intercept, or both) and compare the models using the ANOVA command. This will give you an empirical test of the overall model fit. I was reading in Kirk (Experimental Design) last night and was reminded that examining descriptives, raw data, and visual displays are too often overlooked and dismissed, but are extremely important. Best, HCD -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net -Original Message- From: Andrej Kveder [mailto:[EMAIL PROTECTED] Sent: Wednesday, July 02, 2003 5:41 AM To: Douglas Bates; J.R. Lockwood Cc: Harold Doran; R-Help Subject: RE: [R] within group variance of the coeficients in LME Firstly let me thank all for your answers and suggestions. I would like to followup on your comments. Currently I'm implementing multilevel models within a simulation run. So I was looking for an estimate weather a covariate varies across second level units or not. I was using HLM before so I knew about the test implemented in the package and tried to reporoduce it in R. However I think I can follow your logic about not using that. I would appreciate your reflection on the following. I need a quantitative figure to evaluate weather the covariate varies across second level units in the process of simulation. Of course I will be running thousands of them and would need to program the condition in code. In one of the previous questions to the group dr. Bates suggested to use the CI estmates, however he warned me about their very conservative nature (I got the same tip from the book). I thought about using the lower bound of the CI as an estimate with the rule if above 0 then the covariate varies. Would that be a sound think to do? Do you have any other suggestions? I would really appreciate the feedback. thanks Andrej -Original Message- From: Douglas Bates [mailto:[EMAIL PROTECTED] Behalf Of Douglas Bates Sent: Monday, June 30, 2003 6:09 PM To: J.R. Lockwood Cc: Harold Doran; R-Help; Andrej Kveder Subject: Re: [R] within group variance of the coeficients in LME J.R. Lockwood [EMAIL PROTECTED] writes: Dear listers, I can't find the variance or se of the coefficients in a multilevel model using lme. The component of an lme() object called apVar provides the estimated asymptotic covariance matrix of a particular transformation of the variance components. Dr. Bates can correct me if I'm wrong but I believe it is the matrix logarithm of Cholesky decomposition of the covariance matrix of the random effects. I believe the details are in the book by Pinheiro and Bates. Once you know the transformation you can use the apVar elements to get estimated asympotic standard errors for your variance components estimates using the delta method. J.R. Lockwood 412-683-2300 x4941 [EMAIL PROTECTED] http://www.rand.org/methodology/stat/members/lockwood/ First, thanks to those who answered the question. I have been away from my email for about a week and am just now catching up on the r-help list. As I understand the original question from Andrej he wants to obtain the standard errors for coefficients in the fixed effects part of the model. Those are calculated in the summary method for lme objects and returned as the component called 'tTable'. Try library(nlme) example(lme) summary(fm2)$tTable to see the raw values. Other software for fitting mixed-effects models, such as SAS PROC MIXED and HLM, return standard errors along with the estimates of the variances and covariances of the random effects. We don't return standard errors of estimated variances because we don't think they are useful. A standard error for a parameter estimate is most useful when the distribution of the estimator is approximately symmetric, and these are not. Instead we feel that the variances and covariances should be converted to an unconstrained scale, and preferably a scale for which
RE: [R] within group variance of the coeficients in LME
lme does not produce standard errors for the variance components like HLM does. It does produce SEs for the fixed effects, however, along with t-statistics and p-values, just like HLM. Use the summary() command to see these. When you do this, you will get the AIC, BIC, and loglik values. Just below this output will be the variance components for the random effects. But, the level 2 variance components are reported as standard deviations and SEs do not accompany these random effects. In lme, the residual is the within-group error, which is the sigma-squared in HLM. In terms of lme, the plot(intervals) can be used to assess variability rather than the chi-square in HLM. -Original Message- From: Andrej Kveder [mailto:[EMAIL PROTECTED] Sent: Wed 6/25/2003 5:24 PM To: R-Help Cc: Subject: [R] within group variance of the coeficients in LME Dear listers, I can't find the variance or se of the coefficients in a multilevel model using lme. I want to calculate a Chi square test statistics for the variability of the coefficients across levels. I have a simple 2-level problem, where I want to check weather a certain covariate varies across level 2 units. Pinheiro Bates suggest just looking at the intervals or doing a rather conservative ANOVA test in their book. I have also consultet Raudenbush Bryk on the subject and they suggest using a Chi sqare statistics. It is defined as follows: SUM by j( (beta_hat_qj - y_hat_q0 - sum(y_hat_qs*w_sj))^2/V_hat_qqj) beta are the within 2-level coffecients - got them with the coef() y is a fixed effect or grand mean the sum is for accounting of the level 2 predictors, that I don't presently have, but will the problem is V_hat_qqj which are the variances of the coefficients. I can't get to them. Does anybody have an idea how to get to them? I would really appreciate any suggestion. Andrej _ Andrej Kveder, M.A. researcher Institute of Medical Sciences SRS SASA; Novi trg 2, SI-1000 Ljubljana, Slovenia phone: +386 1 47 06 440 fax: +386 1 42 61 493 [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] R Commander
I am trying to import a file using R Commander. It was working a few days ago, but now I get the following message when I try to import from SPSS. Any thoughts? Error in parse(file, n, text, prompt) : parse error -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] read.spss
I have loaded the foreign package and am still having problems with an import. I get a message that reads, unable to open file. Whe I try different files I get the same message. Here is the code I used. Am I missing something? I am using 1.7 and have also tried this in 1.6 with the same problem. hsb-read.spss(C:\HLM504_Student\Examples\AppendxA\HSB1.SAV, use.value.labels=TRUE, to.data.frame=FALSE, max.value.labels=Inf) -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
[R] Sparse Matrix
I am learning about sparse matrices and wonder if R can create them from a full matrix. Can anyone tell me how I might be able to accomplish this. -- Harold C. Doran Director of Research and Evaluation New American Schools 675 N. Washington Street, Suite 220 Alexandria, Virginia 22314 703.647.1628 http://www.edperform.net/ [[alternate HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
