[R] nlme fixed effects specification
roger koenker
roger at ysidro.econ.uiuc.edu
Sat May 5 17:21:08 CEST 2007
Ivo,
I don't know whether you ever got a proper answer to this question.
Here is a kludgy one -- someone else can probably provide
a more elegant version of my diffid function.
What you want to do is sweep out the mean deviations from both y
and x based on the factor fe and then estimate the simple y on x
linear model.
I have an old function that was originally designed to do panel data
models that looks like this:
"diffid" <- function(h, id)
{
if(is.vector(h))
h <- matrix(h, ncol = 1)
Ph <- unique(id)
Ph <- cbind(Ph, table(id))
for(i in 1:ncol(h))
Ph <- cbind(Ph, tapply(h[, i], id, mean))
is <- tapply(id, id)
Ph <- Ph[is, - (1:2)]
h - Ph
}
With this you can do:
set.seed(1);
fe = as.factor( as.integer( runif(100)*10 ) ); y=rnorm(100); x=rnorm
(100);
summary(lm(diffid(y,fe) ~ diffid(x,fe)))
HTH,
Roger
On May 4, 2007, at 3:08 PM, ivo welch wrote:
> hi doug: yikes. could I have done better? Oh dear. I tried to make
> my example clearer half-way through, but made it worse. I meant
>
> set.seed(1);
> fe = as.factor( as.integer( runif(100)*10 ) ); y=rnorm(100); x=rnorm
> (100);
> print(summary(lm( y ~ x + fe)))
> <deleted>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 0.1128 0.3680 0.31 0.76
> x 0.0232 0.0960 0.24 0.81
> fe1 -0.6628 0.5467 -1.21 0.23
> <deleted more fe's>
> Residual standard error: 0.949 on 89 degrees of freedom
> Multiple R-Squared: 0.0838, Adjusted R-squared: -0.0192
> F-statistic: 0.814 on 10 and 89 DF, p-value: 0.616
>
> I really am interested only in this linear specification, the
> coefficient on x (0.0232) and the R^2 of 8.38% (adjusted -1.92%). If
> I did not have so much data in my real application, I would never have
> to look at nlme or nlme4. I really only want to be able to run this
> specification through lm with far more observations (100,000) and
> groups (10,000), and be done with my problem.
>
> now, with a few IQ points more, I would have looked at the lme
> function instead of the nlme function in library(nlme). [then
> again, I could understand stats a lot better with a few more IQ
> points.] I am reading the lme description now, but I still don't
> understand how to specify that I want to have dummies in my
> specification, plus the x variable, and that's it. I think I am not
> understanding the integration of fixed and random effects in the same
> R functions.
>
> thanks for pointing me at your lme4 library. on linux, version
> 2.5.0, I did
> R CMD INSTALL matrix*.tar.gz
> R CMD INSTALL lme4*.tar.gz
> and it installed painlessly. (I guess R install packages don't have
> knowledge of what they rely on; lme4 requires matrix, which the docs
> state, but having gotten this wrong, I didn't get an error. no big
> deal. I guess I am too used to automatic resolution of dependencies
> from linux installers these days that I did not expect this.)
>
> I now tried your specification:
>
>> library(lme4)
> Loading required package: Matrix
> Loading required package: lattice
>> lmer(y~x+(1|fe))
> Linear mixed-effects model fit by REML
> Formula: y ~ x + (1 | fe)
> AIC BIC logLik MLdeviance REMLdeviance
> 282 290 -138 270 276
> Random effects:
> Groups Name Variance Std.Dev.
> fe (Intercept) 0.000000000445 0.0000211
> Residual 0.889548532468 0.9431588
> number of obs: 100, groups: fe, 10
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) -0.0188 0.0943 -0.199
> x 0.0528 0.0904 0.585
>
> Correlation of Fixed Effects:
> (Intr)
> x -0.022
> Warning messages:
> 1: Estimated variance for factor 'fe' is effectively zero
> in: `LMEoptimize<-`(`*tmp*`, value = list(maxIter = 200L, tolerance =
> 0.0000000149011611938477,
> 2: $ operator not defined for this S4 class, returning NULL in: x
> $symbolic.cor
>
> Without being a statistician, I can still determine that this is not
> the model I would like to work with. The coefficient is 0.0528, not
> 0.0232. (I am also not sure why I am getting these warning messages
> on my system, either, but I don't think it matters.)
>
> is there a simple way to get the equivalent specification for my smple
> model, using lmer or lme, which does not choke on huge data sets?
>
> regards,
>
> /ivo
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list