# LM2 Mallin kaava ja mallimatriisi # # Synteettisten kuitujen lujuus # Selittäjän eli faktorin a tasoina puuvillan %-osuudet kuidussa: # 15,20,25, 30 ja 35 prosenttia # y on vetolujuus y <- scan() 7 7 15 11 9 12 17 12 18 18 14 18 18 19 19 19 25 22 19 23 7 10 11 15 11 cotton <- data.frame(a = factor(rep(1:5,rep(5,5))),y) #rep(1:5,rep(5,5)) # [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 obj <- lm(y ~ a, cotton) alf.star <- coef(obj) alf.star a <- factor(rep(1:5,rep(5,5))) tapply(y,list(a),mean) # 1 2 3 4 5 # 9.8 15.4 17.6 21.6 10.8 # ryhmäkeskiarvot alf.star # alfa-tähti parametrien estimaatit #(Intercept) a2 a3 a4 a5 # 9.8 5.6 7.8 11.8 1.0 # Intercept=9.8 1.rka 2.rka=9.8+5.6 jne Ca <- contrasts(cotton$a) # contrast matrix for `a' #Ca # 2 3 4 5 1 0 0 0 0 2 1 0 0 0 3 0 1 0 0 4 0 0 1 0 5 0 0 0 1 alf <- drop(Ca %*% alf.star[-1]) # alf # 1 2 3 4 5 # 0.0 5.6 7.8 11.8 1.0 dummy.coef(obj) # Full coefficients are # # (Intercept): 9.8 # a: 1 2 3 4 5 # 0.0 5.6 7.8 11.8 1.0 ###################### # Estimoidaan rajoitteen c^T*C_a=0 vallitessa # c^T=(1 0 0 0 0) # Muodostetaan mallimatriisi # ja estimoidaan parametrit rajoitteen aplha_1=0 # suoraan laskemalla pns-laskukaavalla ################### # TEHTÄVÄ 1.14 # ################### X <- matrix(0,25,6) # Mallimatriisin pohja X[,1] <- rep(1,25) rep(2:6,rep(5,5)) i <- cbind(1:25,rep(2:6,rep(5,5))) # indeksimatriisi X[i] <- 1 # contrast.treatment cx <- c(0,1,0,0,0,0) # rajoitevektori #cx #[1] 0 1 0 0 0 0 Xc <- rbind(X,cx) dim(y) <- c(25,1) yc <- rbind(y,0) XXc <- t(Xc)%*%Xc b <- t(Xc)%*%yc beta <- solve(XXc)%*%b # solve(XXc) kääntää matriisin XXc t(beta) # contrast.sum cx <- c(0,1,1,1,1,1) Xc <- rbind(X,cx) yc <- rbind(y,0) XXc <- t(Xc)%*%Xc b <- t(Xc)%*%yc beta <- solve(XXc)%*%b # solve(XXc) kääntää matriisin XXc t(beta) C <- c(-5.24, 0.36, 2.56, 6.56, -4.24) alfa+15.04 # contrast.helmert dim(alfa) <- c(5,1) library(mass) t(fractions(ginv(contr.helmert(n=5)))%*%alfa) # 2.8 1.666667 1.833333 -1.06 ############## # Vastaa alpha.star-parametrisointia. Parametrit # estimoitu rajoitteen alpha_1=0 vallitessa ############## beta <- solve(XXc,b) # Ratkaisee yhtälön XXc%*%beta = b > t(beta) [,1] [,2] [,3] [,4] [,5] [,6] [1,] 9.8 2.517839e-14 5.6 7.8 11.8 1 # Simuloitu data dat <- data.frame(a = factor(rep(1:3, 3)), y = rnorm(9, rep(2:4, 3), 0.1)) rep(1:3, 3) #[1] 1 2 3 1 2 3 1 2 3 rep(2:4, 3) #[1] 2 3 4 2 3 4 2 3 4 dat a y 1 1 1.926767 2 2 3.234859 3 3 3.927577 4 1 1.955336 5 2 3.125241 6 3 3.999333 7 1 1.895134 8 2 3.140310 9 3 4.078987 obj <- lm(y ~ a, dat) alf.star <- coef(obj) alf.star #(Intercept) a2 a3 # 1.925745 1.241058 2.076220 Ca <- contrasts(dat$a) # contrast matrix for `a' alf <- drop(Ca %*% alf.star[-1]) > alf 1 2 3 0.000000 1.241058 2.076220 > dummy.coef(obj) Full coefficients are (Intercept): 1.925745 a: 1 2 3 0.000000 1.241058 2.076220 ################################ # Kontrastimatriisit ################################ ################################ N <- factor(Nlevs <- c(0,1,2,4)) contrasts(N) contrasts(ordered(N)) contrasts(N) 1 2 4 0 0 0 0 1 1 0 0 2 0 1 0 4 0 0 1 op <- options(contrasts=c("contr.helmert", "contr.poly")) op $contrasts unordered ordered "contr.treatment" "contr.poly" ################ # TEHTÄVÄ 1.13 # ################ obj <- lm(y ~ a, cotton) obj #Coefficients: #(Intercept) a2 a3 a4 a5 # 9.8 5.6 7.8 11.8 1.0 options(contrasts=c("contr.sum", "contr.poly")) obj <- lm(y ~ a, cotton) obj #Coefficients: #(Intercept) a1 a2 a3 a4 # 15.04 -5.24 0.36 2.56 6.56 # a5 -4.24 options(contrasts=c("contr.helmert", "contr.poly")) obj <- lm(y ~ a, cotton) obj #Coefficients: #(Intercept) a1 a2 a3 a4 # 15.040 2.800 1.667 1.833 -1.060 ################################# N <- factor(Nlevs <- c(0,1,2,4)) contrasts(N) # Helmertin kontrastit, C-matriisi [,1] [,2] [,3] 0 -1 -1 -1 1 1 -1 -1 2 0 2 -1 4 0 0 3 contr.helmert(n=4) contr.sum(n=4) library(mass) fractions(ginv(contr.helmert(n=4))) fractions(ginv(contr.sum(n=4))) fractions(ginv(contr.sdif(n=4))) contr.sum(n=5) contr.helmert(n=5) # Estimoidaan nyt cotton-aineistosta va-mallin param alf.star (Intercept) a1 a2 a3 a4 15.040000 2.800000 1.666667 1.833333 -1.060000 > fractions(ginv(contrasts(a))) # C^+, C:n g-inverssi [,1] [,2] [,3] [,4] [,5] [1,] -1/2 1/2 0 0 0 [2,] -1/6 -1/6 1/3 0 0 [3,] -1/12 -1/12 -1/12 1/4 0 [4,] -1/20 -1/20 -1/20 -1/20 1/5 mu <- tapply(y,list(a),mean) > alf.star (Intercept) a1 a2 a3 a4 15.040000 2.800000 1.666667 1.833333 -1.060000 > mu 1 2 3 4 5 9.8 15.4 17.6 21.6 10.8 > 0.5*(mu[2]-mu[1]) # alf.star[1] =0.5*(alf_2-alf_1) 2 2.8 > dim(mu) <- c(5,1) > ginv(contrasts(a))%*%mu [,1] [1,] 2.800000 [2,] 1.666667 [3,] 1.833333 [4,] -1.060000 > contrasts(ordered(N)) .L .Q .C 0 -0.6708204 0.5 -0.2236068 1 -0.2236068 -0.5 0.6708204 2 0.2236068 -0.5 -0.6708204 4 0.6708204 0.5 0.2236068 N2 <- N contrasts(N2, 2) <- poly(Nlevs, 2) N2 <- C(N, poly(Nlevs, 2), 2) # alternative contrasts(N2) fractions(ginv(contr.helmert(n = 4))) ############# # Halutaan estimoida erotukset # alf_2-alf_1, alf_3-alf_2, jne ############# Cp <- diag(-1, 4, 5) Cp[row(Cp) == col(Cp) - 1] <- 1 Cp > Cp <- diag(-1, 4, 5) > Cp [,1] [,2] [,3] [,4] [,5] [1,] -1 0 0 0 0 [2,] 0 -1 0 0 0 [3,] 0 0 -1 0 0 [4,] 0 0 0 -1 0 > Cp[row(Cp) == col(Cp) - 1] <- 1 > Cp # Haluttu C-matriisi [,1] [,2] [,3] [,4] [,5] [1,] -1 1 0 0 0 [2,] 0 -1 1 0 0 [3,] 0 0 -1 1 0 [4,] 0 0 0 -1 1 > fractions(ginv(Cp)) [,1] [,2] [,3] [,4] [1,] -4/5 -3/5 -2/5 -1/5 [2,] 1/5 -3/5 -2/5 -1/5 [3,] 1/5 2/5 -2/5 -1/5 [4,] 1/5 2/5 3/5 -1/5 [5,] 1/5 2/5 3/5 4/5 > options(contrasts=c("contr.sdif", "contr.poly")) > contrasts(a) 2-1 3-2 4-3 5-4 1 -0.8 -0.6 -0.4 -0.2 2 0.2 -0.6 -0.4 -0.2 3 0.2 0.4 -0.4 -0.2 4 0.2 0.4 0.6 -0.2 5 0.2 0.4 0.6 0.8 > fractions(contrasts(a)) 2-1 3-2 4-3 5-4 1 -4/5 -3/5 -2/5 -1/5 2 1/5 -3/5 -2/5 -1/5 3 1/5 2/5 -2/5 -1/5 4 1/5 2/5 3/5 -1/5 5 1/5 2/5 3/5 4/5 > fractions(ginv(contrasts(a))) [,1] [,2] [,3] [,4] [,5] [1,] -1 1 0 0 0 [2,] 0 -1 1 0 0 [3,] 0 0 -1 1 0 [4,] 0 0 0 -1 1 obj <- lm(y ~ a, cotton) alf.star <- coef(obj) > alf.star (Intercept) a2-1 a3-2 a4-3 a5-4 15.04 5.60 2.20 4.00 -10.80 > t(Cp%*%mu) [,1] [,2] [,3] [,4] [1,] 5.6 2.2 4 -10.8 ############# # Halutaan estimoida erotukset # alf_i-mean(alf), i=1,...,5 eli poikkeamat kok.keskiarvosta ############# > options(contrasts=c("contr.sum", "contr.poly")) > contrasts(a) [,1] [,2] [,3] [,4] 1 1 0 0 0 2 0 1 0 0 3 0 0 1 0 4 0 0 0 1 5 -1 -1 -1 -1 > fractions(ginv(contrasts(a))) [,1] [,2] [,3] [,4] [,5] [1,] 4/5 -1/5 -1/5 -1/5 -1/5 [2,] -1/5 4/5 -1/5 -1/5 -1/5 [3,] -1/5 -1/5 4/5 -1/5 -1/5 [4,] -1/5 -1/5 -1/5 4/5 -1/5 > obj <- lm(y ~ a, cotton) > alf.star <- coef(obj) > alf.star (Intercept) a1 a2 a3 a4 15.04 -5.24 0.36 2.56 6.56 [1] 4.24 # a5 = 4.24 > t(mu) [,1] [,2] [,3] [,4] [,5] [1,] 9.8 15.4 17.6 21.6 10.8 > mean(mu) [1] 15.04 > obj <- lm(y ~ a-1, cotton) > coef(obj) a1 a2 a3 a4 a5 9.8 15.4 17.6 21.6 10.8 > obj <- lm(y ~ a, cotton) > gm <- rep(1,5)*mean(mu) > gm [1] 15.04 15.04 15.04 15.04 15.04 > alf.star[-1] a1 a2 a3 a4 -5.24 0.36 2.56 6.56 alf <- c(alf.star[-1],-sum(alf.star[-1])) > alf a1 a2 a3 a4 -5.24 0.36 2.56 6.56 -4.24 > sum(alf) [1] -1.665335e-16 > gm+alf # kok.ka + alpha a1 a2 a3 a4 9.8 15.4 17.6 21.6 10.8