Philippe Massicotte - 10 months ago 57

R Question

I am using the

`lme4`

`lmer()`

`SAS`

`PARMS`

In the following example, the estimated variances are:

`c(0.00000, 0.03716, 0.00000, 0.02306)`

I would like to

`c(0.09902947, 0.02460464, 0.05848691, 0.06093686)`

so there are not estimated.

`> summary(mod1)`

Linear mixed model fit by maximum likelihood ['lmerMod']

Formula: log_cumcover_mod ~ (1 | kildestationsnavn) + (1 | year) + (1 |

kildestationsnavn:year) + (1 | proevetager)

Data: res

AIC BIC logLik deviance df.resid

109.9 122.9 -48.9 97.9 59

Scaled residuals:

Min 1Q Median 3Q Max

-2.1056 -0.6831 0.2094 0.8204 1.7574

Random effects:

Groups Name Variance Std.Dev.

kildestationsnavn:year (Intercept) 0.00000 0.0000

kildestationsnavn (Intercept) 0.03716 0.1928

proevetager (Intercept) 0.00000 0.0000

year (Intercept) 0.02306 0.1518

Residual 0.23975 0.4896

Number of obs: 65, groups:

kildestationsnavn:year, 6; kildestationsnavn, 3; proevetager, 2; year, 2

Fixed effects:

Estimate Std. Error t value

(Intercept) 4.9379 0.1672 29.54

Answer Source

This is possible, if a little hacky. Here's a reproducible example:

Fit the original model:

```
library(lme4)
set.seed(101)
ss <- sleepstudy[sample(nrow(sleepstudy),size=round(0.9*nrow(sleepstudy))),]
m1 <- lmer(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
fixef(m1)
## (Intercept) Days
## 251.55172 10.37874
```

Recover the deviance (in this case REML-criterion) function:

```
dd <- as.function(m1)
```

I'm going to set the standard deviations to zero so that I have something to compare with, i.e. the coefficients of the regular linear model. (The parameter vector for `dd`

is a vector containing the column-wise, lower-triangular, concatenated Cholesky factors for the random effects terms in the model. Luckily, if all you have are scalar/intercept-only random effects (e.g. `(1|x)`

), then these correspond to the random-effects standard deviations.)

```
(ff <- dd(c(0,0))) ## new REML: 1704.708
environment(dd)$pp$beta(1) ## new parameters
## [1] 251.11920 10.56979
```

Matches:

```
coef(lm(Reaction~Days,ss))
## (Intercept) Days
## 251.11920 10.56979
```

If you want to construct a new `merMod`

object you can do it as follows ...

```
opt <- list(par=c(0,0),fval=ff,conv=0)
lmod <- lFormula(Reaction~Days+(1|Subject)+(0+Days|Subject),ss)
m1X <- mkMerMod(environment(dd), opt, lmod$reTrms, fr = lmod$fr,
mc = quote(hacked_lmer()))
```

If you want to set the variances to (700,30), then follow the same procedure but start with `dd(sqrt(700),sqrt(30))`

. Calling the deviance function sets the variances internally *and* re-estimates the fixed-effect parameters internally ...

By the way, it's not generally recommended to use random effects when the grouping variables have a very small number (e.g. <5 or 6) levels: see here ...