0

I estimate the following model for estimating the standard errors for random coefficients:

library(nlme)
library(foreign)
Sample_Data <- read.dta (file="Sample_Data.dta")
Sample_Data <- within(Sample_Data, {
industry <- factor(industry)
id <- factor(id)
year <- factor(year)})
Model1 <- lme(MRPK ~ 1+age+empl+net_worth+input_growth+TFPR_shock, random = (~1+net_worth|industry/id), data=Sample_Data,  control=lmeControl(opt='optim'), method="REML")

My sample data on the following link:

https://dl.dropboxusercontent.com/u/99295510/Sample_Data.dta

I am getting the following error:

Error in solve.default(pdMatrix(a, factor = TRUE)) :
Lapack routine dgesv: system is exactly singular: U[1,1] = 0

How could we solve this problem?

mjalam
  • 29
  • 4
  • What are you trying to achieve? The error message says that you can't use the `intervals` function on an lmerMod object. this probably means that the `intervals` package does not support objects produced by `lme4`. – lmo May 03 '16 at 19:37
  • My objective is to estimate the standard errors for random coefficient. I do not know how to estimate the standard error using lme4. I thought I could calculate confidence intervals then I can calculate standard error from the confidence intervals. – mjalam May 03 '16 at 19:51
  • Maybe this answer will help you: [lmer CI](http://stackoverflow.com/questions/11072544/how-to-get-coefficients-and-their-confidence-intervals-in-mixed-effects-models/17329983#17329983). – lmo May 03 '16 at 19:55

1 Answers1

0

At present I don't know of a reliable way to get the Wald standard errors (the only sensible definition I know of) for random effects using lme4. Your alternatives are (1):

variation.model1  <- lmer(sales ~ 1+(1|industry:id)+(1|industry),
    data=data_file, REML=TRUE)
confint(variation.model1, which="theta_")

which will give you 95% profile likelihood confidence intervals (but not standard errors) on the variance parameters (very crudely, you could divide the width of those confidence intervals by 1.96*2=3.92).

And (2), fit the model with lme instead:

library(nlme)
variation.model2  <- lme(sales ~ 1,
        random = (1|industry/id),
        data=data_file, method="REML")

Built-in (no need for intervals package):

intervals(variation.model2)     ## Wald intervals
sqrt(diag(variation.model2$apVar))  ## sd

See this rpubs document for more details about Wald variances in lme4.

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Thanks for your nice feedback. I need to estimate the residual variance of each year for the above model. Could I use the residuals of lme and calculate variance of each year? – mjalam May 08 '16 at 02:18
  • Thanks again. I need to write for the first case confint(variation.model1, parm="theta_"). However, this produces missing values. I got the following errors "There were 50 or more warnings (use warnings() to see the first 50)". For the second case, I also getting errors: " Error in diag(variation.model2$apVar) : invalid 'nrow' value (too large or NA)" "In addition: Warning message: In diag(variation.model2$apVar) : NAs introduced by coercion". What should I do now? Even though I add one more random variable in the model. – mjalam Jun 18 '16 at 12:51
  • Can you please include data and/or add code to your question that will provide me with a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) ? – Ben Bolker Jun 19 '16 at 01:17
  • I have added head of data. – mjalam Jun 19 '16 at 05:57
  • I would like to estimate the following model - – mjalam Jul 17 '16 at 01:13
  • I add the data and model. Please tell me how could I handle this problem. – mjalam Jul 17 '16 at 01:29
  • Do you have time to say something with my model? Thanks! – mjalam Jul 26 '16 at 18:43