Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Not able to run simex with model fit by nls() #10

Open
gjmzajac opened this issue Oct 25, 2022 · 0 comments
Open

Not able to run simex with model fit by nls() #10

gjmzajac opened this issue Oct 25, 2022 · 0 comments

Comments

@gjmzajac
Copy link

The simex documentation claims that it can be used to correct for measurement error in models fit by nls() but I encounter an error when I try to do so. I am attempting to use nls() to fit the same model and data as in the help page for the simex function.

> library(simex)
> 
> ## Seed
> set.seed(49494)
> 
> ## simulating the measurement error standard deviations
> sd_me <- 0.3
> sd_me2 <- 0.4
> temp <- runif(100, min = 0, max = 0.6)
> sd_me_het1 <- sort(temp)
> temp2 <- rnorm(100, sd = 0.1)
> sd_me_het2 <- abs(sd_me_het1 + temp2)
> 
> ## simulating the independent variables x (real and with measurement error):
> 
> x_real <- rnorm(100)
> x_real2 <- rpois(100, lambda = 2)
> x_real3 <- -4*x_real + runif(100, min = -10, max = 10)  # correlated to x_real
> 
> x_measured <- x_real + sd_me * rnorm(100)
> x_measured2 <- x_real2 + sd_me2 * rnorm(100)
> x_het1 <- x_real + sd_me_het1 * rnorm(100)
> x_het2 <- x_real3 + sd_me_het2 * rnorm(100)
> 
> ## calculating dependent variable y:
> y <- x_real + rnorm(100, sd = 0.05)
> y2 <- x_real + 2*x_real2 + rnorm(100, sd = 0.08)
> y3 <- x_real + 2*x_real3 + rnorm(100, sd = 0.08)
> 
> ### one variable with homoscedastic measurement error
> (model_real <- lm(y ~ x_real))

Call:
lm(formula = y ~ x_real)

Coefficients:
(Intercept)       x_real  
  -0.004752     0.994554  

> 
> (model_naiv <- lm(y ~ x_measured, x = TRUE))

Call:
lm(formula = y ~ x_measured, x = TRUE)

Coefficients:
(Intercept)   x_measured  
   -0.01217      0.90327  

> 
> (model_simex <- simex(model_naiv, SIMEXvariable = "x_measured", measurement.error = sd_me))

Naive model:
lm(formula = y ~ x_measured, x = TRUE)

SIMEX-Variables: x_measured
Number of Simulations: 100

Coefficients:
(Intercept)   x_measured  
   -0.01528      0.97007  

> 
> #for nls
> (nls_naiv = nls(y ~ intercept + beta * x_measured, start = list("intercept" = -0.01217, "beta" = 0.90327), model=T))
Nonlinear regression model
  model: y ~ intercept + beta * x_measured
   data: parent.frame()
intercept      beta 
 -0.01217   0.90327 
 residual sum-of-squares: 8.326

Number of iterations to convergence: 1 
Achieved convergence tolerance: 2.973e-08
> (nls_simex <- simex(nls_naiv, SIMEXvariable = "x_measured", measurement.error = sd_me, asymptotic = F))
Error in dimnames(x) <- dn : 
  length of 'dimnames' [2] not equal to array extent
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant