Bagaimana cara menguji pengaruh variabel pengelompokan dengan model non-linear?

15

Saya punya pertanyaan tentang penggunaan variabel pengelompokan dalam model non-linear. Karena fungsi nls () tidak memungkinkan untuk variabel faktor, saya telah berjuang untuk mencari tahu apakah seseorang dapat menguji efek dari faktor pada model fit. Saya telah memasukkan contoh di bawah ini di mana saya ingin mencocokkan model pertumbuhan "musiman von Bertalanffy" dengan perlakuan pertumbuhan yang berbeda (paling umum diterapkan pada pertumbuhan ikan). Saya ingin menguji efek danau tempat ikan tumbuh serta makanan yang diberikan (hanya contoh buatan). Saya akrab dengan solusi untuk masalah ini - menerapkan F-test membandingkan model yang cocok untuk data yang dikumpulkan vs cocok terpisah seperti yang diuraikan oleh Chen et al. (1992) (ARSS - "Analisis jumlah residu kuadrat"). Dengan kata lain, untuk contoh di bawah ini,

masukkan deskripsi gambar di sini

Saya membayangkan ada cara sederhana untuk melakukan ini di R menggunakan nlme (), tapi saya mengalami masalah. Pertama-tama, dengan menggunakan variabel pengelompokan, derajat kebebasan lebih tinggi daripada yang saya dapatkan dengan pemasangan model terpisah saya. Kedua, saya tidak dapat membuat variabel pengelompokan sarang - Saya tidak melihat di mana masalah saya. Setiap bantuan menggunakan nlme atau metode lain sangat dihargai. Di bawah ini adalah kode untuk contoh buatan saya:

###seasonalized von Bertalanffy growth model
soVBGF <- function(S.inf, k, age, age.0, age.s, c){
    S.inf * (1-exp(-k*((age-age.0)+(c*sin(2*pi*(age-age.s))/2*pi)-(c*sin(2*pi*(age.0-age.s))/2*pi))))
}

###Make artificial data
food <- c("corn", "corn", "wheat", "wheat")
lake <- c("king", "queen", "king", "queen")

#cornking, cornqueen, wheatking, wheatqueen
S.inf <- c(140, 140, 130, 130)
k <- c(0.5, 0.6, 0.8, 0.9)
age.0 <- c(-0.1, -0.05, -0.12, -0.052)
age.s <- c(0.5, 0.5, 0.5, 0.5)
cs <- c(0.05, 0.1, 0.05, 0.1)

PARS <- data.frame(food=food, lake=lake, S.inf=S.inf, k=k, age.0=age.0, age.s=age.s, c=cs)

#make data
set.seed(3)
db <- c()
PCH <- NaN*seq(4)
COL <- NaN*seq(4)
for(i in seq(4)){
    age <- runif(min=0.2, max=5, 100)
    age <- age[order(age)]
    size <- soVBGF(PARS$S.inf[i], PARS$k[i], age, PARS$age.0[i], PARS$age.s[i], PARS$c[i]) + rnorm(length(age), sd=3)
	PCH[i] <- c(1,2)[which(levels(PARS$food) == PARS$food[i])]
	COL[i] <- c(2,3)[which(levels(PARS$lake) == PARS$lake[i])]
	db <- rbind(db, data.frame(age=age, size=size, food=PARS$food[i], lake=PARS$lake[i], pch=PCH[i], col=COL[i]))
}

#visualize data
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)


###fit growth model
library(nlme)

starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)

#fit to pooled data ("small model")
fit0 <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values
)
summary(fit0)

#fit to each lake separatly ("large model")
fit.king <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values,
  subset=db$lake=="king"
)
summary(fit.king)

fit.queen <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values,
  subset=db$lake=="queen"
)
summary(fit.queen)


#analysis of residual sum of squares (F-test)
resid.small <- resid(fit0)
resid.big <- c(resid(fit.king),resid(fit.queen))
df.small <- summary(fit0)$df
df.big <- summary(fit.king)$df+summary(fit.queen)$df

F.value <- ((sum(resid.small^2)-sum(resid.big^2))/(df.big[1]-df.small[1])) / (sum(resid.big^2)/(df.big[2]))
P.value <- pf(F.value , (df.big[1]-df.small[1]), df.big[2], lower.tail = FALSE)
F.value; P.value


###plot models
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)
legend("topleft", legend=c("soVGBF pooled", "soVGBF king", "soVGBF queen"), col=c(1,2,3), lwd=2)

#plot "small" model (pooled data)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100))
pred <- predict(fit0, tmp)
lines(tmp$age, pred, col=1, lwd=2)

#plot "large" model (seperate fits)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="king")
pred <- predict(fit.king, tmp)
lines(tmp$age, pred, col=2, lwd=2)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="queen")
pred <- predict(fit.queen, tmp)
lines(tmp$age, pred, col=3, lwd=2)



###Can this be done in one step using a grouping variable?
#with "lake" as grouping variable
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit1 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake,
  start=starting.values
)
summary(fit1)

#similar residuals to the seperatly fitted models
sum(resid(fit.king)^2+resid(fit.queen)^2)
sum(resid(fit1)^2)

#but different degrees of freedom? (10 vs. 21?)
summary(fit.king)$df+summary(fit.queen)$df
AIC(fit1, fit0)


###I would also like to nest my grouping factors. This doesn't work...
#with "lake" and "food" as grouping variables
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit2 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake/food,
  start=starting.values
)

Referensi: Chen, Y., Jackson, DA dan Harvey, HH, 1992. Perbandingan fungsi von Bertalanffy dan polinomial dalam memodelkan data pertumbuhan ikan. 49, 6: 1228-1235.

r mixed-model nls Marc di dalam kotak
sumber

Jawaban:

6

X1,...,XpYf

Y=f(X1,...,Xp)+ε

ε∼N(0,σ2)fBmBL1L0

Model yang tidak terstruktur jelas merupakan submodel dari model yang bertingkat, sehingga uji rasio kemungkinan sesuai untuk melihat apakah model yang lebih besar bernilai kompleksitas yang ditambahkan - statistik uji adalah

λ=2(L1−L0)

λχ2mp−p=p(m−1)pχ2

Makro
sumber
Apakah Anda menyarankan agar sesuai dengan model yang terpisah, jumlah kemungkinan log dari setiap L1 = SUM (LL_i, i dari 1 hingga m) dan kemudian melanjutkan dengan kemungkinan itu? Juga, apakah L0 model dengan prediktor kategorik dalam pertanyaan dimasukkan (dengan variabel dummy m-1 misalnya)?
B_Miner
Ya, saya menyarankan itu. L.0BB
Terima kasih atas saran Anda Makro. Ini tampaknya sesuai dengan apa yang telah saya lakukan - walaupun Anda menyarankan perbandingan kemungkinan daripada uji-F. Dalam contoh saya, uji-F juga membandingkan residu kecocokan tunggal dengan penjumlahan residual dari beberapa kecocokan yang diterapkan pada setiap tingkat prediktor kategori. Saya kira saya bertanya-tanya apakah seseorang dapat melakukan ini dalam model campuran dalam satu langkah daripada mencocokkan beberapa model. Juga, akankah strategi semacam itu memungkinkan pengujian faktor bersarang?
Marc di dalam kotak
Saya tidak berpikir Anda akan dapat menyesuaikan beberapa model untuk membandingkan model. Juga, ya, uji rasio kemungkinan dapat digunakan untuk menguji faktor-faktor bersarang.
Makro
2

Saya menemukan bahwa adalah mungkin untuk kode variabel variabel dengan nls (), hanya dengan mengalikan vektor benar / salah ke dalam persamaan Anda. Contoh:

# null model (no difference between groups; all have the same coefficients)
nls.null <- nls(formula = percent_on_cells ~ vmax*(Time/(Time+km)),
            data = mehg,
            start = list(vmax = 0.6, km = 10))

# alternative model (each group has different coefficients)
nls.alt <- nls(formula = percent_on_cells ~ 
              as.numeric(DOC==0)*(vmax1)*(Time/(Time+(km1))) 
            + as.numeric(DOC==1)*(vmax2)*(Time/(Time+(km2)))
            + as.numeric(DOC==10)*(vmax3)*(Time/(Time+(km3)))
            + as.numeric(DOC==100)*(vmax4)*(Time/(Time+(km4))),
            data = mehg, 
            start = list(vmax1=0.63, km1=3.6, 
                         vmax2=0.64, km2=3.6, 
                         vmax3=0.50, km3=3.2,
                         vmax4= 0.40, km4=9.7))
housetyrell
sumber