Proses menghasilkan data adalah:
Misalkan menjadi urutan dari ke dengan panjang dan menjadi faktor yang sesuai . Ambil semua kemungkinan kombinasi untuk menghitung :
Menggunakan B-spline-Basis (tidak terpusat) untuk untuk setiap tingkat tidak akan layak oleh properti-of-unity-properti (baris jumlah ke 1). Model seperti itu tidak akan dapat diidentifikasi (bahkan tanpa intersep).
Contoh: (Pengaturan: 5 interval simpul dalam (terdistribusi secara merata), B-Spline derajat 2, spline
-fungsi adalah kebiasaan)
# drawing the sequence
n <- 100
x <- seq(-4,4,length.out=n)
z <- seq(-4,4,length.out=n)
d <- as.factor(0:1)
data <- CJ(x=x,z=z,d=d)
set.seed(100)
# setting up the model
data[,y := sin(x+I(d==0)) + sin(x+4*I(d==1)) + I(d==0)*z^2 + 3*I(d==1)*z^2 + rnorm(n,0,1)]
# creating the uncentered B-Spline-Basis for x and z
X <- data[,spline(x,min(x),max(x),5,2,by=d,intercept=FALSE)]
> head(X)
x.1d0 x.2d0 x.3d0 x.4d0 x.5d0 x.6d0 x.7d0 x.1d1 x.2d1 x.3d1 x.4d1 x.5d1 x.6d1 x.7d1
[1,] 0.5 0.5 0 0 0 0 0 0.0 0.0 0 0 0 0 0
[2,] 0.0 0.0 0 0 0 0 0 0.5 0.5 0 0 0 0 0
[3,] 0.5 0.5 0 0 0 0 0 0.0 0.0 0 0 0 0 0
Z <- data[,spline(z,min(z),max(z),5,2,by=d)]
head(Z)
z.1d0 z.2d0 z.3d0 z.4d0 z.5d0 z.6d0 z.7d0 z.1d1 z.2d1 z.3d1 z.4d1 z.5d1 z.6d1
[1,] 0.5000000 0.5000000 0.00000000 0 0 0 0 0.0000000 0.0000000 0.00000000 0 0 0
[2,] 0.0000000 0.0000000 0.00000000 0 0 0 0 0.5000000 0.5000000 0.00000000 0 0 0
[3,] 0.4507703 0.5479543 0.00127538 0 0 0 0 0.0000000 0.0000000 0.00000000 0 0 0
z.7d1
[1,] 0
[2,] 0
[3,] 0
# lm will drop one spline-column for each factor
lm(y ~ -1+X+Z,data=data)
Call:
lm(formula = y ~ -1 + X + Z, data = data)
Coefficients:
Xx.1d0 Xx.2d0 Xx.3d0 Xx.4d0 Xx.5d0 Xx.6d0 Xx.7d0 Xx.1d1 Xx.2d1 Xx.3d1 Xx.4d1 Xx.5d1
23.510 19.912 18.860 22.177 23.080 19.794 18.727 68.572 69.185 67.693 67.082 68.642
Xx.6d1 Xx.7d1 Zz.1d0 Zz.2d0 Zz.3d0 Zz.4d0 Zz.5d0 Zz.6d0 Zz.7d0 Zz.1d1 Zz.2d1 Zz.3d1
69.159 67.496 1.381 -11.872 -19.361 -21.835 -19.698 -11.244 NA -1.329 -38.449 -62.254
Zz.4d1 Zz.5d1 Zz.6d1 Zz.7d1
-69.993 -61.438 -39.754 NA
Untuk mengatasi masalah ini, Wood, Generalized Additive Models: An Introduction with R , halaman 163-164 mengusulkan jumlah (atau rata-rata) batasan keterpusatan:
Ini dapat dilakukan dengan reparametrization jika matriks ditemukan sedemikian rupa
-matrix dapat ditemukan dengan dekomposisi-QR dari matriks kendala .
Perhatikan bahwa adalah oleh partisi unity-property.
Versi terpusat / terbatas dari B-Spline-Matrix saya adalah:
X <- data[,spline(x,min(x),max(x),5,2,by=d,intercept=TRUE)]
head(X)
x.1d0 x.2d0 x.3d0 x.4d0 x.5d0 x.6d0 x.1d1 x.2d1 x.3d1 x.4d1
[1,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000000 0.2271923 -0.3225655 -0.3225655 -0.3225655
[3,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000 0.0000000 0.0000000 0.0000000
x.5d1 x.6d1
[1,] 0.0000000 0.00000000
[2,] -0.2728077 -0.05790256
[3,] 0.0000000 0.00000000
Z <- data[,spline(z,min(z),max(z),5,2,by=d,intercept=TRUE)]
head(Z)
z.1d0 z.2d0 z.3d0 z.4d0 z.5d0 z.6d0 z.1d1 z.2d1 z.3d1 z.4d1
[1,] 0.2271923 -0.3225655 -0.3225655 -0.3225655 -0.2728077 -0.05790256 0.0000000 0.0000000 0.0000000 0.0000000
[2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000000 0.2271923 -0.3225655 -0.3225655 -0.3225655
[3,] 0.2875283 -0.3066501 -0.3079255 -0.3079255 -0.2604260 -0.05527458 0.0000000 0.0000000 0.0000000 0.0000000
z.5d1 z.6d1
[1,] 0.0000000 0.00000000
[2,] -0.2728077 -0.05790256
[3,] 0.0000000 0.00000000
Pertanyaan saya adalah: Meskipun kecocokannya sangat mirip, mengapa kolom-kolom B-Spline saya berbeda dari apa yang disediakan oleh gam? Apa yang saya lewatkan?
# comparing with gam from mgcv
mod.gam <- gam(y~d+s(x,bs="ps",by=d,k=7)+s(z,bs="ps",by=d,k=7),data=data)
X.gam <- model.matrix(mod.gam)
head(X.gam)
(Intercept) d1 s(x):d0.1 s(x):d0.2 s(x):d0.3 s(x):d0.4 s(x):d0.5 s(x):d0.6 s(x):d1.1 s(x):d1.2
1 1 0 0.5465301 -0.05732768 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000 0.00000000
2 1 1 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000 0.00000000 0.5465301 -0.05732768
3 1 0 0.5465301 -0.05732768 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000 0.00000000
s(x):d1.3 s(x):d1.4 s(x):d1.5 s(x):d1.6 s(z):d0.1 s(z):d0.2 s(z):d0.3 s(z):d0.4 s(z):d0.5
1 0.0000000 0.0000000 0.0000000 0.00000000 0.5465301 -0.057327680 -0.2351708 -0.2259983 -0.1201207
2 -0.2351708 -0.2259983 -0.1201207 -0.01043987 0.0000000 0.000000000 0.0000000 0.0000000 0.0000000
3 0.0000000 0.0000000 0.0000000 0.00000000 0.5471108 -0.031559945 -0.2302910 -0.2213227 -0.1176356
s(z):d0.6 s(z):d1.1 s(z):d1.2 s(z):d1.3 s(z):d1.4 s(z):d1.5 s(z):d1.6
1 -0.01043987 0.0000000 0.000000000 0.0000000 0.0000000 0.0000000 0.00000000
2 0.00000000 0.5465301 -0.057327680 -0.2351708 -0.2259983 -0.1201207 -0.01043987
3 -0.01022388 0.0000000 0.000000000 0.0000000 0.0000000 0.0000000 0.00000000
Garis putus-putus sesuai dengan pas saya, garis lurus ke versi gam
Jawaban:
Berikut adalah contoh sederhana menggunakan tautan dari Nemo. Pertanyaan yang saya jawab adalah
Saya menjawab ini karena ini adalah judul dan judul
agak tidak jelas karena alasan yang saya berikan pada akhirnya. Inilah jawaban untuk pertanyaan di atas
Anda dapat melakukan lebih baik dalam hal kecepatan komputasi daripada komputasi secara eksplisit
seperti yang dijelaskan pada halaman 211 dari
Ada beberapa masalah dalam kode OP
Untuk
maka saya tidak mengerti bagaimana Anda berharap untuk mendapatkan yang sama. Anda mungkin telah menggunakan simpul yang berbeda dan saya tidak melihat bagaimana
spline
fungsi akan menghasilkan hasil yang benar di sini.Jika yang terakhir dipasang
lm
maka tidak dikenakan sanksi sehingga hasilnya harus berbeda?sumber
spline
-fungsi adalah kebiasaan