R Classical Statistics: Estimation, Tests, Linear Models, and GLMs

This English version was machine-translated from the Chinese original. Technical terms may require verification.

EDA and descriptive statistics

And here we'll introduce the parts of the EDA technology that are mainly descriptive, and there are other parts that are also relevant.

Probability Map (PDF)

Although we keep saying that PDF and the CDF have a similar role to play, we use PDF as the main thing, and we do it when we make the drawings, and we don't really study the CDF because it's not intuitive enough.

Density function graphs of the usual distribution

Understanding the pattern of the overall distribution helps to capture the basic characteristics of the sample. We first look at some of the probabilistic functions of the commonly used distribution mentioned in chapter III through specific examples (see below). For the discrete distribution point method, for the continuous distribution means its density function. We use the functions of the PDF to make our graphics work.

Here we have some R code to help us understand the usual PDF maps.

Two distributions

n<-20
p<-0.2
k<-seq(0,n)
plot(k,dbinom(k,n,p),type='h', main='Binomial distribution, n=20, p=0.2',xlab='k')

Porcelain distribution

lambda<-4.0
k<-seq(0,20)
plot(k,dpois(k,lambda),type='h', main='Poisson distribution, lambda=5.5',xlab='k')

Geometric distribution

p<-0.5
k<-seq(0,10)
plot(k,dgeom(k,p),type='h', main='Geometric distribution, p=0.5',xlab='k')

Supergeometric distribution

N<-30
M<-10
n<-10
k<-seq(0,10)
plot(k,dhyper(k,N,M,n),type='h', main='Hypergeometric distribution,
     N=30, M=10, n=10',xlab='k')

Negative binary distribution

n<-10
p<-0.5
k<-seq(0,40)
plot(k, dnbinom(k,n,p), type='h',
     main='Negative Binomial distribution,
     n=10, p=0.5',xlab='k')

Normal distribution

curve(dnorm(x,0,1), xlim=c(-5,5), ylim=c(0,.8),col='red', lwd=2, lty=3)
curve(dnorm(x,0,2), add=T, col='blue', lwd=2, lty=2)
curve(dnorm(x,0,1/2), add=T, lwd=2, lty=1)
title(main="Gaussian distributions")
legend(par('usr')[2], par('usr')[4], xjust=1, c('sigma=1', 'sigma=2', 'sigma=1/2'),
       lwd=c(2,2,2), lty=c(3,2,1),col=c('red', 'blue', par("fg")))

t Distribution

curve(dt(x,1), xlim=c(-3,3), ylim=c(0,.4),col='red', lwd=2, lty=1)
curve(dt(x,2), add=T, col='green', lwd=2, lty=2)
curve(dt(x,10), add=T, col='orange', lwd=2, lty=3)
title(main="Student T distributions")
legend(par('usr')[2], par('usr')[4], xjust=1, c('df=1', 'df=2', 'df=10', 'Gaussian distribution'),
lwd=c(2,2,2,2), lty=c(1,2,3,4),
col=c('red', 'blue', 'green', par("fg")))

Carside distribution

curve(dchisq(x,1), xlim=c(0,10), ylim=c(0,.6), col='red', lwd=2)
curve(dchisq(x,2), add=T, col='green', lwd=2)
curve(dchisq(x,3), add=T, col='blue', lwd=2)
title(main='Chi square Distributions')

F distribution

curve(df(x,1,1), xlim=c(0,2), ylim=c(0,.8), lty=1)
curve(df(x,3,1), add=T, lwd=2,lty=2)
curve(df(x,6,1), add=T, lwd=2, lty=3)
title(main="Fisher's F")

Histogram and density function estimation

Histogram

Histogram is the basic tool for exploratory data analysis, giving a frequency distribution diagram of data, with a long rectangle of equal widths commonly used in group distance situations, where the rectangle represents the size of the frequency; on the graphic, the cross-references represent the range of values to be taken from the variable of interest, and the coordinates indicate the frequency (or frequency) size, so that the frequency (or frequency) is the right direction. Figure

hist(x, breaks = "Sturges", freq = NULL, probability = !freq,
     col = NULL, main = paste("Histogram of" , xname),
     xlim = range(breaks), ylim = NULL,
     xlab = xname, ylab, axes = TRUE, nclass = NULL)

Where the breaks are used to specify partitions (integer numbers are number of spaces) col indicates colours freq indicates whether to use frequency-numbers in the square Figure

Nuclear density estimates

density(x, bw = "nrd0",
        kernel = c("gaussian", "epanechnikov", "rectangular",
                   "triangular", "biweight", "cosine", "optcosine"),
        n = 512, from, to)

kernel decides that the smooth function is used by default is normal; n gives the number of nuclear density estimates at equal intervals from the left and right ends of the nuclear density estimates to be calculated separately from to Points

Example

N <- 100000
n <- 100
p <- .9
x <- rbinom(N,n,p)
hist(x, xlim=c(min(x),max(x)), probability=T,
     nclass=max(x)-min(x)+1, col='lightblue',
     main='Binomial distribution, n=100, p=.5')
lines(density(x,bw=1), col='red', lwd=3)

Descriptive statistical analysis

Descriptive statistical analysis is an important part of the EDA, and we're here to present it, in terms of the data-type gap.

Knowledge on graphics is available for referenceR Visualization

Descriptive statistical analysis of individual data sets

Graphical missions

The distribution of the single group data can be done by the histograms described above, as well as by the nucleodensity curve and box line. Normality tests are generally done using QQ graphs.

library(DAAG)
data(possum)
fpossum <- possum[possum$sex=="f",]
par(mfrow=c(1,2))
attach(fpossum)
hist(totlngth,breaks=72.5+(0:5)*5, ylim=c(0,22),
       xlab="total length", main="A:Breaks at 72.5,77.5…")

stem(fpossum$totlngth)
boxplot(fpossum$totlngth)
qqnorm(fpossum$totlngth, main="Normality Check via QQ Plot")
qqline(fpossum$totlngth, col='red')

Data missions

We've been introduced to the Obsidian, and it's enough to give a simple description of some code.

library(DAAG)
data(possum)
fpossum <- possum[possum$sex=="f",]

趋势部分
summary(fpossum$totlngth) #汇总分析
fivenum(fpossum$totlngth) #五数分布
quantile(fpossum$totlngth) #分位数
median(fpossum$totlngth) #中位数
max(fpossum$totlngth)
min(fpossum$totlngth)
mean(fpossum$totlngth) #极大 极小 均值
#离散部分
max(fpossum$totlngth)-min(fpossum$totlngth)
IQR(fpossum$totlngth)
sd(fpossum$totlngth)
var(fpossum$totlngth)
mad(fpossum$totlngth)
#偏度峰度
library(fBasics)
skewness(fpossum$totlngth)
kurtosis(fpossum$totlngth)

Special statistical kitfBasicsFunctions inbasicStats( )It provides almost all descriptive statistics,pastecsThere's a name in the bag.stat.desc()And the function of the function has the same effect.

Descriptive statistical analysis of multiple sets of data

Graphical missions

We're more likely to have a few ways of looking at multiple data.

n<-10
d<-data.frame(y1 = abs(rnorm(n)),
y2 = abs(rnorm(n)),
y3 = abs(rnorm(n)),
y4 = abs(rnorm(n)),
y5 = abs(rnorm(n)) )
plot(d)
matplot(d, type = 'l', ylab = "", main = "Matplot")
boxplot(d)

Data missions

The data we need to look at is not much more than the single set of data, but the original function changes some of the methods used to promote it, and then adds a little bit of content.

In particular, we may need to revert to descriptive statistical analysis functions for a column for the data box, and we prefer to use functions. apply() and sapply() aggregate() They can increase programming efficiency. Apply function to matrix and data box

## 汇总分析和一些针对数据框的操作函数
summary(state.x77)
aggregate(state.x77, list(Region = state.region), mean)
aggregate(state.x77, list(Region = state.region, Cold = state.x77[,"Frost"] > 130),mean)

sd(state.x77)
#var函数不可以继续计算方差了
#相关分析
x<-c(44.4, 45.9, 46.0, 46.5, 46.7, 47, 48.7, 49.2, 60.1)
y<-c(2.6, 10.1, 11.5, 30.0, 32.6, 50.0, 55.2, 85.8, 86.8)
cor(x,y)
cor(x,y,method="spearman")
cor(x,y,method="kendall")
cor.test(x,y, method="spearman")

cor（）The function calculates three correlations at the same time, as Pearson Spearman Kendall, we simply do the description here, and then we study the correlation separately, and we also use the correlation as a hypothetical test.

Descriptive statistics for disaggregated data

List of rows

If the variables to which the data are focused are qualitative, the data are referred to as disaggregated data, which are often described in tables and serve further statistical analysis, and we consider mainly the data in the D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D-D

A very simple code allows for a column table to be created.

## 本质上就是构造了一个有行名和列名的矩阵，不能能用数据框替代
Eye.Hair <- matrix(c(68,20,15,5, 119,84,54,29, 26,17,14,14, 7,94,10,16),  nrow=4,byrow=T)
colnames(Eye.Hair) <- c("Brown", "Blue", "Hazel", "Green")
rownames(Eye.Hair) <- c("Black","Brown","Red", "Blond")
Eye.Hair

You can also construct a column combination from the original data using the table function

## table()函数从因子factor中获取频数的函数
## 当接受两个factor的时候，table函数创建一个二维的列联表，高维也可
table(menarche,tanner)

A joint list for the Govi ftable()function can output multi-dimensional arrays in a compact and attractive way, reference functions

ftable(table(factorA,factorB,factorC))

The margins of the list are very important, and we have a function to study it; besides that, it's also important to study the frequency list.

## 创建列联表
Eye.Hair <- matrix(c(68,20,15,5, 119,84,54,29, 26,17,14,14, 7,94,10,16), nrow=4,byrow=T)
colnames(Eye.Hair) <- c("Brown", "Blue", "Hazel", "Green")
rownames(Eye.Hair) <- c("Black","Brown","Red", "Blond")

1 按照行 2 按照列 计算边缘总和
margin.table(Eye.Hair,1)
margin.table(Eye.Hair,2)
1 按照行 2 按照列 计算边缘概率
prop.table(Eye.Hair,1)
prop.table(Eye.Hair,2)
#直接获取概率矩阵
Eye.Hair/sum(Eye.Hair)

Simple graphical description

The bar chart is the most basic graphic depiction of the matrix, and it's very common to have Marseic charts, which we'll present in a separate graphic section. Marcello.

data(HairEyeColor)
a <- as.table(apply(HairEyeColor,c(1,2),sum))
barplot(a, legend.text = attr(a, "dimnames")$Hair)
barplot(a, beside = TRUE, legend.text = attr(a, "dimnames")$Hair)

Centralization and standardization of data

scale(x, center = TRUE, scale = TRUE)
#中心化or标准化函数 x接受矩阵和数据框 返回中心化or标准化后的结果
#第一个参数控制中心化 第二个参数控制缩放

Statistical extrapolations

Here we're going to introduce the part of the classic statistics that we're talking about statistical inferences, which is to include parameter estimates, parameter hypothesis tests, and a few more simple non-parametric hypothesis tests, which, while we spent a whole semester on basic knowledge, we're basically good enough to use a chapter on R to achieve it.

Parameter estimation

Rectangular and very similar estimates

Rectangular estimate

By the law of Sinchin and the strong Como Gorov theorem, if the overall X-K-Rect exists, the k-Rect of the sample is reduced by probability to the overall k-Rect, and the continuous function of the sample-Rect is reduced to the continuous function of the general rect.

We don't have to study any theory here, and one example is enough.

An index distribution with a total parameter of λ, the density function is

$$ p(x|\lambda)=\lambda\exp^{-\lambda x},\quad x>0 $$

then$\lambda$The rule is estimated.

$$ \hat{\lambda}=\frac{1}{\overline{X}} $$

X<-c(0.59132754,0.12854935,0.46900228,0.29835980,0.24341462, 0.06566637,0.40085536,2.99687123,0.05278912,0.09898594)
lambda<- 1/mean(X)
lambda

So this is what the rectangular estimate does, is the original mathematical calculations are part of R's calculations, or is it a simple part of R's calculations?

Very seemingly estimated (presentation of R optimized function)

So we need to calculate the apparent function and then use software to do a hugely enhanced processing, which means that what we're actually doing in R is an optimisation problem.

#optimize( )的调用格式
optimize(f = , interval = , lower = min(interval),
        upper = max(interval), maximum = TRUE,
        tol = .Machine$double.eps^0.25, ...)

of which

• f is an apparent function that requires us to define a form of basic knowledge. Pass.
• Interval is the value of the parametric gill; the lower is the bottom of the gill, upper is the upper of the gill.
• maxim = TRUE is for a large value, otherwise (maximm = FALSE) indicates the very small value of the function.
• Tol is the exact value requested, and it's normal to default.

Optimize only applies to the optimization of single parameters, but it applies to both the maximum and the minimum.

#nlm( )的调用格式
nlm(f, p, hessian = FALSE, typsize=rep(1, length(p)),     fscale=1,print.level = 0, ndigit=12, gradtol = 1e-6,
stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000),
steptol = 1e-6, iterlim = 100, check.analyticals = TRUE, ...)

It uses the Newton-Lafson algorithm to find the minimum point of the function.

#optim( )的调用格式
optim(par, fn, gr = NULL,
method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"),
lower = -Inf, upper = Inf, control = list( ), hessian = FALSE, ...)

Optimize one of the five methods given by the method option

The last two can be used for multidimensional issues.

And here we're going to give an example of a very similar estimate, which is a one-dimensional function that is given directly.

f <- function(P){(P^517)*(1-P)^483}
optimize(f,c(0,1),maximum = TRUE)

It's operating in two parts.

• maximm is a very similar estimate.
• objective is the function of the function at this time

Inter-sectional estimates of mono-normal overall parameters

Inter-sectional estimates are a relatively special kind of problem, and he has a very close connection to the hypothetical tests, because the core and the statistical data are in a relatively close form.

In fact, many of the problem with inter-sectional estimates is that the results of the inter-sectional estimates, which are performed by the hypothetical function, are part of the results of the hypothetical test function.

But the function of some questions is not provided by R, which in fact does not have to design functions to assist our solution, and we'll just go back to the problem of R providing an estimate of the function.

The difference is unknown, and the most common estimate of the average is the average value of the study.

$$ \begin{pmatrix}\overline{X}-\frac{S}{\sqrt{n}}t_{1-\frac{\alpha}{2}}(n-1),\overline{X}+\frac{S}{\sqrt{n}}t_{1-\frac{\alpha}{2}}(n-1)\end{pmatrix} $$

It's easy to know that we're doing a t-test.

t.test(x, y = NULL,alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, ...)

• X, y is the data used for the tests, and all give is a double samplet.
• Alternative decision-type
• mu is average, only works with hypothetical tests.

Here's a simple example.

x<-c(175 , 176 , 173 , 175 ,174 ,173 , 173, 176 , 173,179 )
t.test(x)
t.test(x)$conf.int
#可以用conf.int选择只访问置信区间 本质上就是在列表上选择了一部分分量

Inter-sectional estimates of two normal overall parameters

The mathematical form of the range estimate for the difference between the two aggregates, which is unknown but equal, is the average difference.

$$ \left((\overline{X}-\overline{Y})\pm t_{1+\frac\alpha2}\sqrt{\frac1n+\frac1m}\sqrt{\frac{(n-1)S_1^2+(m-1)S_2^2}{n+m-2}}\right) $$

And at this point, we know that what is needed is a test t. Using the t test function described above, we can use the following examples:

x<-c(628,583,510,554,612,523,530,615)
y<-c(535,433,398,470,567,480,498,560,503,426)
t.test(x,y,var.equal=TRUE)
#var.equal需要设定TRUE 此时认为两个总体方差相等

The study equation is also a common feature of the two normals in general.

var.test(x, y, ratio = 1,
         alternative = c("two.sided", "less", "greater"),conf.level = 0.95, ...)

One example is as follows:

x<-c(20.5,19.8,19.7,20.4,20.1,20.0,19.0,19.9)
y<-c(20.7,19.8,19.5,20.8,20.4,19.6,20.2)
var.test(x,y)

Inter-area estimates for single-total ratio p

In many practical questions, we often have to estimate the proportion of individuals with certain characteristics in the overall population, which is a category that deserves separate research, which is very important, and the mathematical theory generally uses large samples to get a normal distribution in the form of the following:

$$ \hat{p}\pm z_{1-\frac{\alpha}{2}}\sqrt{\hat{p}(1-\hat{p})/n}-\frac{1}{2n}. $$

We have a special R function to do this. The sample properties are almost subject to hypergeometric distributions, and we can choose to use either normal or two distributions to match them.The latter is a different mathematical form.

#调用格式
prop.test(x, n, p = NULL,
          alternative = c("two.sided", "less", "greater"),
          conf.level = 0.95, correct = TRUE)
#正态检验是一种近似检验 需要大样本
binom.test(x, n, p = NULL,
          alternative = c("two.sided", "less", "greater"),
          conf.level = 0.95, correct = TRUE)
#二项分布是一种精确检验 不需要大样本 这里我们本质上是调用了二项分布的估计和检验

• x is the number of samples n is the total number
• Correct is whether to use a continuous distribution approximation
• P is the probability of the original hypothesis, which is useful in the hypothetical test.

Estimated inter-area differences in the two overall ratios

In the case of large samples, they are almost subject to normal distribution, so they can give the formula that has been used to make a case for the distribution of the normal.

$$ (\hat{p}_1-\hat{p}2)\pm z{1-\frac{\alpha}{2}}\sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2}.} $$

We have only one form to deal with this type of problem.

like<-c(478, 246)
people<-c(1000, 750)
prop.test(like, people)

Sample capacity determination

This is a counter-question of a range of estimates, and we give the maximum permissible error of the parameter estimates, and we calculate the sample capacity required; mathematically, there's no change in the calculation, but there are some new functions that are needed to deal with the problem.

R base does not provide any function that would help us deal with this type of problem, suggesting manual extrapolation, using R-based functions to do some ancillary computation; of course, it can also ask if there are any packages that help us with this type of problem when needed.

Parameter hypothetical tests

Another important element of statistical extrapolation is the hypothetical tests, the information provided by the sample, the construction of the appropriate statistical quantity, and the testing of the assumptions provided;

Assuming that the tests are divided into two broad categories of parameters, the overall distribution at this time (classical mathematical studies) is known as the case of the parameters we study;

We're testing the distribution type for a detailed study in non-parametric statistics, and we need more theoretical support to understand it better.

Some important theoretical knowledge in the hypothetical test

Our core step is

1. Give the original assumptions, get the corresponding alternative assumptions.
2. Determination of the level of visibility$\alpha$
3. Study the numbers, determine his distribution.
4. Gives a field of rejection of the original rejection (statistical down to the original rejection)
5. Calculates the value of the test statistics to which the sample points correspond
6. Disclaimer of the original hypothesis

- Assuming it's a test.$p$Value is a different kind of hypothetical test from our previous studies; at this point we do not determine the level of prominence in advance, but we calculate it.$p$ Value: To decide whether to reject the original hypothesis by comparing it with the common situation of several levels of prominence

$p$The smaller the value, the more it should be rejected, the more common it is to have a few significant levels of 0.1 0.05 0.01, and every less, the more significant trend is to reject the original assumptions, and the more important that we are in the later hypothesis test is that we are looking at the following:$p$Value

We're just going to go back and introduce some of the more important hypothetical questions.

Test for single-normal overall parameters

The hypothetical test of the single normal overall average is easy to verify.

$$ T=\frac{\overline{X}-\mu_{0}}{S^{*}/\sqrt{n}} \sim t(n-1) $$

And we're still calling the functions we've described.

salt<-c(490 , 506, 508, 502, 498, 511, 510, 515 , 512)
t.test(salt, mu=500)

And the mu is the average of the original assumption, which is the same function as the estimated spatially of the average, and in fact, they do the same thing, and if they need to, we can do a single-sided examination, and adjust the parameters inside.

Test of two normal overall parameters

The test of whether the two normal general averages are zero can also be conducted directly as a hypothetical test.

x<-c(628,583,510,554,612,523,530,615)
y<-c(535,433,398,470,567,480,498,560,503,426)
t.test(x,y,var.equal=TRUE)
#var.equal需要设定TRUE 此时两个总体方差相等 否则我们需要修改参数为FALSE

The problem with the difference can also be calculated directly from the function in front, the assumption is still one.

x<-c(20.5,19.8,19.7,20.4,20.1,20.0,19.0,19.9)
y<-c(20.7,19.8,19.5,20.8,20.4,19.6,20.2)
var.test(x,y)

Test of data in pairs t

We can use the t.test function that we've already described, but we just need to modify some parameters, and the difference doesn't matter at this point, but the pairs need to be declared individually.

x<-c(20.5, 18.8, 19.8, 20.9, 21.5, 19.5, 21.0, 21.2)
y<-c(17.7, 20.3, 20.0, 18.8, 19.0, 20.1, 20.0, 19.1)
t.test(x, y, paired=TRUE)

Test of single sample ratio

We're still ahead of us, with two kinds of precision tests and a proximate test.

binom.test(c(7, 5), p=0.4)
#这是另一种调用格式 输入向量 前面是成功次数 后面是失败次数
prop.test(7, 12, p=0.4, correct=TRUE)
#近似检验可能不准 R会在任何需要的时候给出警告

Test of two sample ratios

In fact, it's very simple, exactly the same pattern as the previous parameter estimates.

like<-c(478, 246)
people<-c(1000, 750)
prop.test(like, people)

Relevant analysis

Relevance of classification variables

First we need to test the independence of the classification variables, that is, the independence test of the matrix.

R offers a variety of things.Test for type variable independenceThe way we're going to be here to introduce you.

The Carabineros Independence Test

Availablechisq.test()function to test the column and column variables of the 2-dimensional table for the card independence of the column, as shown below

compare<-matrix(c(60,32,3,11), nr = 2, dimnames = list(c("cancer", "normal"),c("smoke", "Not smoke")))

chisq.test(compare, correct=TRUE) #检验函数接受一个二维表

Fisher's exact inspection.

The previous C.I.S. test allows only 20% of the two-dimensional list to have an expected frequency of less than 5 or a warning that we can use the Fisher test, instead of a similar calibration, to be used on a two-dimensional list of any number of columns greater than 2 but not for a warning.$2\times2$

compare<-matrix(c(60,32,3,11), nr = 2, dimnames = list(c("cancer", "normal"),c("smoke", "Not smoke")))

fisher.test(compare, correct=TRUE)

Cochran-Mantel-Haenszel Test

mantelhaen.test()The function is used to perform the Cochran-Mantel-Haenszel card-check, the assumption being,Two nominal variables are independent in each of the third variables.I'm sorry. The following codes test the independence of treatment and improvement at each level of gender.

## 构建列联表
mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis)
## 检验
mantelhaen.test(mytable)

Calculation of relevant coefficients

The distinguishing test in the previous section assessed whether there was sufficient evidence to reject the original assumption that the variables were independent of each other. If you can reject the assumption, then your interest will naturally shift to the measure of relevance that measures the strength and weakness of relevance.vcdPackageassocstats()function to calculate the phi coefficient, the column link and the Cramer's V coefficient of the 2-dimensional array table. Example code

library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
assocstats(mytable)

Pearson, Spearman and Kendall are related.

Pearson ' s coefficients measure linear correlation between two quantitative variables. The Specarman grade-related coefficient measures the degree of correlation between the hierarchical sequence variables. Kendall's Tau-related coefficient is also a grade-related measure of non-parameters. He's also our most basic three related coefficients. Availablecorfunction calculates,methonParameters are used to select the method of calculation, as follows:

states<- state.x77[,1:6]
#cor函数可以对一个数据框进行操作
cor(states)

After calculating the coefficient, we should also consider the issue of visibility, which can be tested using the Cor.test() function for individual Pearson, Spearman and Kendall-related coefficients. The simplified format for use is:

cor.test(x, y, alternative = , method = )

Offset factor

Relevance refers to the interrelationship between the other two quantitative variables when controlling one or more quantitative variables. You can use it.ggmPackagepcor()function calculates the relative coefficient The function calls in the format:

pcor(u, S)

Of whichuis a numerical vector, the first two values represent the subscript of the variable to which the relevant coefficient is to be calculated, and the remaining values are the subscript of the conditional variable (i.e. the variable to which the impact is to be excluded).Sis the synapse array of variables. Example:

library(ggm)
colnames(states)
pcor(c(1,5,2,3,6), cov(states))

Corresponds, we have a hypothetical test function.

pcor.test(r, q, n)

Of whichrBypcor()function calculates the coefficient of bias,qFor the number of variables (in place of value),nIt's the size of the sample.

Difference Analysis

The equation analysis (analysis of analysis, short of ANOVA) is an effective statistical method for analysing experimental data in industrial and agricultural production and scientific research.

There are two main causes of differences in observations (variability): uncontrolled fluctuations caused by random factors or observations errors during the testing, and manageable fluctuations caused by different treatments in the testing or by different conditions.

The main work of the variance analysis is to decompose the total variation (variant) of the observation data into factor effects and test errors according to the causes of the variation and to provide quantitative analysis of the factors, comparing the importance of the various causes in the total variation as a basis for further statistical inferences.

ANOVA is widely used in various experimental and quasi-experimental designs, and is another regression model in the case of a classification variable.

Single factor variance analysis

Basic methodology

We've already introduced specific mathematical models in our pilot design, and here's what R should do.

aov(formula, data=NULL, projections=FALSE,qr=TRUE, contrasts=NULL, ...)

And we'll be back to tell you how this formula should be used.

Gives a single factor, five-level example of a differential analysis, and this code should be able to understand the meaning.

X<-c(25.6, 22.2, 28.0, 29.8, 24.4, 30.0, 29.0, 27.5, 25.0, 27.7,
23.0, 32.2, 28.8, 28.0, 31.5, 25.9, 20.6, 21.2, 22.0, 21.2)
A<-factor(rep(1:5, each=4))
#rep函数 times参数是控制整体重复的次数 each是控制每个元素重复的次数
A
miscellany<-data.frame(X, A)
aov.mis<-aov(X~A, data=miscellany)
#X列标识了数据 A列标识了数据对应的因子水平
summary(aov.mis)

plot(miscellany$X ~ miscellany$A)
#绘制分组图形的箱线图，直观的比较差异
plot(miscellany$A, miscellany$X)
#这个和上面是一样的 自变量为因子的plot函数也会绘制boxplot

Multiple comparisons of average values

When we do the differential analysis, we find that there are significant differences in the mean values of the effects, and we can only know that there are certain values that differ, but we can't tell which is different, and the following approach helps us.Find out which values are different when the variance is analysed.

It's actually a comparison between two levels of a factor A, using a test test that's not fundamentally different from the test we're using in front of us.

However, in the course of multiple t tests, the p-values need to be adjusted to ensure normal judgment, as follows:

p.adjust.methods

He'll tell us what the p-value adjustment is, and there's more than a Bonferroni, and now we're doing multiple t tests, and here's the way we do it.

pairwise.t.test(x, g, p.adjust.method=p.adjust.methods,pool.sd=TRUE, ...)

xIt's the vector that responds to the variable. gis the group vector (factor) p.adjust.methodIt's the way to adjust the p value mentioned above.

The amount of the matrix returned is allpValue

The Tukey Act

We rejected the assumption that there was no difference between levels and determined who was different. Now we want to make a trust gap in the effects. It's actually another way to determine what effects are different.

function in the

TukeyHSD(x, which, ordered=FALSE, conf.level=0.95)

x It's the result of the variance analysis. which is the factor vector that needs to be calculated between the comparative zones ordered is the logical value, if"true", the factor level is a confidence level

We'll use a simple example to explain the use.

sales<-data.frame( X=c(23, 19, 21, 13, 24, 25, 28, 27, 20, 18, 19, 15, 22, 25, 26, 23, 24, 23, 26, 27), A=factor(rep(1:5, c(4, 4, 4, 4, 4))) )
#数据集生成
summary(aov(X~A, sales))
#方差分析
pairwise.t.test(sales$X, sales$A, p.adjust.method="bonferroni")
#单组比较
TukeyHSD(aov(X~A, sales))
#计算所有均值差的计算区间 因为我们没改which
plot(TukeyHSD(fit))
#绘图

Sqincy test

To do the differential analysis, we need to make sure that the following three conditions are met.

• Addability (variability plus)
• Independent normality (level of independence, internal normality)
• Sq. Alignness (equals of different horizontal differences)

Now we're studying how to test each other for compatibility.

In fact, regression analysis is a relaxed equation, but it's not the same. So we usually only study the differences in the equation analysis.

The normality of the disability can be analysed directly by reference to the issue of return

Bartlett Test

Function format is that the meaning of the parameters need not be explained, as many places in the front line.

bartlett.test(x, g, ...)
bartlett.test(formula, data, subset, no.action, ...)

Levene's testing.

Function format is the same as the argument

leveneTest(x, group)

Example

bartlett.test(X~A, data=sales)
library(car)
leveneTest(sales$X, sales$A)

Both p-values are very large, which means they didn't reject the original hypothesis, which is the equation.

Remarks

The variance analysis model can be considered a special linear model, so the differential analysis can also use linear model functions.lm( ), and also functionanova( )Extract the variance analysis table, thereforeaov(formula)Equivalent toanova(lm(formula)) A single factor differential analysis can also be used for functionsoneway.test( ), if the difference between the data below each level is equal (use options)var.equal=TRUE) , it is equivalent to using a functionaov( )perform normal variance analysis; if the differences in the data below each level are not equal (use options)var.equal=FALSE(iii) Welch (1951), which uses the approximation method;

When the distribution below each level is unknown (generally default normal, the above method is used), the variance analysis is performed using Kruskal-Wallis, etc. They're non-parametric statistics.

Double Factor Difference Analysis

No interaction

Theoretically, we've omitted the following code lines.

juice<-data.frame(
X = c(0.05, 0.46, 0.12, 0.16, 0.84, 1.30, 0.08, 0.38, 0.4, 0.10, 0.92, 1.57, 0.11, 0.43, 0.05, 0.10, 0.94, 1.10, 0.11, 0.44, 0.08, 0.03, 0.93, 1.15), A = gl(4, 6),B = gl(6, 1, 24)
)
#数据建立
juice.aov<-aov(X~A+B, data=juice)
summary(juice.aov)
#方差分析 把公式部分改了
bartlett.test(X~A, data=juice)
bartlett.test(X~B, data=juice)
#两个方差齐性检验

Situations with interactive effects Down

Theoretically, we've omitted the following code lines.

rats<-data.frame(
Time=c(0.31, 0.45, 0.46, 0.43, 0.82, 1.10, 0.88, 0.72, 0.43, 0.45, 0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.38, 0.29, 0.40, 0.23, 0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.40, 0.56, 1.02, 0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.30, 0.37, 0.38, 0.29, 0.23, 0.25, 0.24, 0.22, 0.30, 0.36, 0.31, 0.33),
Toxicant=gl(3, 16, 48, labels = c("I", "II", "III")),
Cure=gl(4, 4, 48, labels = c("A", "B", "C", "D")) )
#数据集构建
op<-par(mfrow=c(1, 2))
#设置绘图参数并存储
plot(Time~Toxicant+Cure, data=rats)
#一类特殊的boxplot绘制方法
with(rats,interaction.plot(Toxicant, Cure, Time, trace.label="Cure"))
with(rats, interaction.plot(Cure, Toxicant, Time, trace.label="Toxicant"))
#绘制交互效应图 如果不出现明显交叉就基本认为没有交互作用
rats.aov<-aov(Time~Toxicant*Cure, data=rats)
summary(rats.aov)
#虽然认为没有交互 但是还是进行了带有交互效应的方法分析 上面的函数等价于
## rats.aov<-aov(Time~Toxicant+Cure+Toxicant:Cure, data=rats)
## 只考虑交互效应则为
## rats.aov<-aov(Time~Toxicant:Cure, data=rats)

The results of the variance analysis suggest that two factors are significant, but the interaction is not significant.

Coordinated differential analysis

The hypothetical tests of comparison of two or more groups of average values in the range analysis methods described in the preceding two sections, which are generally manageable, are sometimes handled by factors that are not controlled in practice, how to deduct or balance the effects of these uncontrollable factors while comparing the differences between two or more groups of equal values, and the method of co-ordinated differential analysis may be considered.

Analysis of Covariance, ancovaA statistical analysis methodology combining linear regression analysis with differential analysisThe underlying idea is to consider some variables (i.e. unknown or uncontrollable) that have an impact on the response variable Y as covariate, to create linear regression relationships that respond to Y and the change in X, and to use this regression to equalize X values before hypothetically testing the difference between the values of the revised Y for each processing group, which is the substance of the reduction of X-Y squared from the overall equation of Y, and to analyse the difference squared and further decomposed before the squared analysis of the equation is made to better evaluate the effects of this treatment.

We use functions to explain the problem. Use the function in the package HH ancova as

ancova(formula, data.in = sys.parent(),x, groups)

Formula is the formula for the coordinated differential analysis, Data.in is the data, x is the variance analysis, groups is the factor.

Here are some examples.

feed<-as.factor(rep(c("A","B","C"),each=8) )
Weight_Initial <- c(15,13,11,12,12,16,14,17,17,16,
                    18,18,21,22,19,18,22,24,20,23, 25,27,30,32)
Weight_Increment <-c(85,83,65,76,80,91,84,90,97,90,
                     100,95,103,106,99,94,89,91,83,
                     95,100,102,105,110)
data_feed<-data.frame(feed,Weight_Initial,Weight_Increment)
#数据集构建 其中Weight_Initial是我们想要考虑的协变量
ancova(Weight_Increment ~ Weight_Initial+feed , data=data_feed)
#不考虑交互
ancova(Weight_Increment ~ Weight_Initial*feed , data=data_feed)
#考虑交互

Test design and variance analysis in progress

Checklist Test

Let's give an example of how data from the active test was read into the data box.

rate<-data.frame(
A=gl(3,3),
B=gl(3,1,9), C=factor(c(1,2,3,2,3,1,3,1,2)),
Y=c(31, 54, 38, 53, 49, 42, 57, 62, 64)
)
#正交试验数据的建立 存储了每个数据对应的各个因子水平
K<-matrix(0, nrow=3, ncol=3, dimnames=list(1:3, c("A","B","C")))
for (j in 1:3)
  for (i in 1:3)
    K[i,j]<-mean(rate$Y[rate[j]==i])
#计算了每个水平的均值（实际上可以使用tapply函数简化 如下）
#K <- tapply(rate$Y, rate$A, mean)
plot(as.vector(K), axes=F, xlab="Level", ylab="Rate")
xmark<-c(NA,"A1","A2","A3","B1","B2","B3","C1","C2","C3",NA)
axis(1,0:10,labels=xmark)
axis(2,4*10:16)
axis(3,0:10,labels=xmark)
axis(4,4*10:16)
lines(K[,"A"]); lines(4:6, K[,"B"]); lines(7:9,K[,"C"])
#因子各个水平 指标均值情况

This is the visual method of analyzing the direct test, which is actually the best indicator available, but not as rigorous as the differential analysis.

I'm just trying to get a difference.

We didn't change the R function.

rate.aov<-aov(Y~A+B+C, data=rate)
summary(rate.aov)

The test is also an analysis of the interaction.

Repeat the experiment.

The so-called duplicate measurement differential analysis, i.e. the testee was measured more than once. This section focuses on the analysis of the difference in measurements (a common design) with a group and an inter-group factor.

Because the variable is the carbon dioxide absorption (uptake) variable is the plant type Type and the CO2 concentration (conc) at seven levels, Type is the inter-group factor, conc is the intra-group factor, Plant is the individual symbol

We have a separate approach to this problem involving repeated tests, and we first ask for data structure.Still one line per observation, which requires simultaneous intra-group factors, inter-cluster factors and individual symbols And the variance analysis code becomes

## A组内因子 W组内因子 B组间因子
y ~ A + Error (Subject/A) #单因素组内 ANOVA
y ~ B * W + Error (Subject/W) #含单个组内因子（w）和单个组间因子（B）的重复测量ANOVA

Regressive analysis

The analysis only leads to a correlation between two variables, but it does not answer how they relate to each other, i.e., they cannot identify the function of a causal relationship between them.

OLS and derivatives

Symbol used in R expression

Here is the more complex system of functions that we first came into contact with R, and there are symbols that we need to understand.

We're here to create a control sheet to show the more common symbols.

Symbol	Expression Meaning
~	Separator, left to respond to variable, right to interpret variable
+	Separating projection variables
:	Intersections that represent the projection variables
*	It's not useful to suggest a simple way to all possible interactive items.
^	Means that the interactive item reaches a certain number, for example $\text{代码} y \sim(x + z + w)^2 \text{可展开为} y\sim x+ z + w + x:z + x:w + z:w$
-1	Remove Intersection
I()	Explain the bracketed elements from the arithmetical point of view, avoiding symbol conflicts

One-linear regression

Math forms need not be wasted on the presentation of the basic but important function of regression, which is the basis for our continuing learning of more functions.

#回归函数
lm(formula, data, subset, weights, na.action,method="qr",
   model=TRUE, x=FALSE, y=FALSE, qr=TRUE, singular.OK=TRUE, contrasts=NULL, offset)

#返回模型的参数
coefficients(object)
#模型参数的置信区间
confint(object, level=0.95, …)
#汇总分析函数
summary(object)
#返回预测残差
residuals()
rstandard(model, infl=lm.influence(model, do.coef=FALSE),
          sd=sqrt(deviance(model)/df.residual(model)), …)
rstudent(model, infl=lm.influence(model), do.coef=FALSE)
#列出拟合模型的预测值
fitted()
#预测
predict(object, newdata, interval = "confidence", level = 0.95)
#绘制回归曲线图 一般和plot联合使用
abline(object)
#手动计算p-value
f_statistic <- summary(X.lm)$fstatistic
f_value <- f_statistic[1]
p_value <- pf(f_value, f_statistic[2], f_statistic[3], lower.tail = FALSE)

lm()Function is the core function of the equation to which the regression equation is based

• Formula is the choice of regression models.
• Data is a data frame
• Subset is a subset of sample observations.
• Weights are weighted vectors for the assembly
• na.action shows whether the data contains missing values
• Method is pointing out the method used to make the match.
• The logical value behind is whether to return the value

summary（）The information that is used to answer the entire model is the classic summy function of the system.

• Model parameter estimates
• Model hypothetical test (without including the equation p value that can be directly quoted, but with the f value that can be used to design a function to address this)

Multiple Returns

We can easily construct a multi-form return form, as follows:

fit2 <- lm(weight ~ height + I(height^2), data=women)
## I的含义是里面增加了一个算术项 我们构建是一个多元回归

It's still a linear regression, a multi-form regression.

In fact,Whatever we construct on the right side of the equation, as long as the parameter item is linear, it does not affect the properties of linear regression.

Multi-linear regression and variable selection

The function has no formal change here.

lm.reg<-lm(y~x1+x2+x3+x4, data=blood)

If you want to study interaction, the code should be changed to

fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)

And here we have a little extra question about the selection of variables.

step(object, scope, scale=0,direction=c("both", "backward", "forward"),trace=1, keep=NULL, steps=1000, k=2)

Parameter Interpretation

• object is the result of a linear model or a broad linear model analysis
• The scope means whether or not to limit the scope of model selection
• Direction is the choice of the method, forward, backward or back.

This function changes the model directly.

Re-entry diagnosis

The main elements are: disability analysis, impact analysis, colinear diagnosis; we're still here to study regression diagnosis, because at this point in time, many of our research is really focused on the previous presentation of the OLS, and some of the elements can certainly be extrapolated seamlessly, and can be used well in other models, and diagnosis, especially in the case of disability analysis, is a very important link in the overall return, and it deserves very deep study.

Standard analytical methodology

The disability analysis is a very large module, and we're here to present only some of the more common parts.

The most detailed disability analysis should be based on residuals（）function to perform a custom analysis, we do not present here.

The most basic disability analysis is based on plot He's been providing us with four of the most common residual analysis. Figure

fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)

It contains the most popular graphical tool for model diagnosis.

• Function of the difference pair of y
• Disability qq graph testing
• Distribution of square root of standardized residuals
• Cook Distance

So we'll study it separately.

• Disability and the independence of variables, i.e. whether important regression items are missing
• Normality of the disability
• Whether the disability is equal to the difference
• Impact analysis issues

Impact analysis

We have a lot of ways to study impacts. But none of these can be directly concluded by the fact that we're not doing anything to make a specific analysis of the problem. The following is the text:

lm.influence(model, do.coef=TRUE)

It gives a model regression factor after a certain point of observation, which can be used to determine the impact.

cooks.distance(model, infl=im.influence(model, do.coef=FALSE),
               res=weighted.residuals(model), sd=sqrt(deviance(model)/df.residual(model)),
               hat=infl$hat, ...)

Cook statistics are also very common statistics on the impact of judgement

dffits(model, infl=..., res=...)

This is the DFFITS Code.

covratio(model, infl=lm.influence(model, do.coef=FALSE),res=weighted.residuals(model))

COVRATIO Guidelines

influence.measures(model)

A summary of the statistics on impact, including the above.

library(car)
influencePlot(fit, id.method="identify", main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")

carA function provided by the package that integrates information on the discrete points, leverage values and powerful impact points in a visualized chart Medium

Symmetrical linear diagnosis

The more common combination of linear diagnostics is the characteristic value kappa, the differential expansion factor VIF, which functions as follows:

eigen(x, symmetric, only.values=FALSE, EISPACK=FALSE)
#计算矩阵的特征值 辅助判断复共线性
kappa(x, exact=FALSE, ...)
#计算矩阵的kappa值 也就是条件数 100以上就是强相关 30以上中度相关
vif(lmobj, digits=5)
#计算方差膨胀因子VIF 这个函数来自于DAGG包 10意味着强相关

Comprehensive test for regression

gvlma()The functions, which were prepared by Pena and Slate (2006), allow for a comprehensive validation of linear model assumptions, together with an evaluation of slopes, peaks and heterogeneity. In other words, it provides a separate comprehensive test (passed/not passed) for model assumptions from the package gvlma

It's very convenient to use.

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

GLM Return

Here's where we're going to discuss the return of the GLM.

On broad linear models

One of the broad linear models (Generalized Linear Model) is the promotion of normal linear models, which require that the response variable is only dependent on interpretation of the variable in linear form. It is understandable here that GLM retains the structure of the linear projection sub-program, but links the response variable expectations to the explanation variable linear combination with a connecting function, while allowing the response variable to come from the index distribution group.

R directly provides the function of matching and calculating broad linear modelsglm( ), it is called in the format

log<-glm(formula, family=family.generator,data=data.frame)

• formulaFor the purpose of formulating formulae, the meaning is the same as that of linear models;
• familyFor the distribution group, including normal distribution, bi-biNMial, porpoise distribution and gamma distribution, the distribution group can also specify the connecting functions to be used by using option link=, refer to the table below
• Data is the data box.

The common distribution family and default connection functions in GLM can be understood by following the following table:

Distribution Type	Default Connect Functions	Common model
binomial	logit	Logit Return
gaussian	identity	Normal linear regression
gamma	inverse	Gamma GLM
inverse.gaussian	1/mu^2	Anti-Gorgos GLM
poisson	log	Porcelain returns.
quasi	identity Difference with constant	Aquarified GLM
quasibinomial	logit	Paradispersion Logit
quasipoisson	log	Accurate Porcelain

A few examples of the way in which the distributional community is called are given below.

#正态分布 恒等连接
fm <- glm(formula, family = gaussian(link = identity), data = data.frame)

#二项分布 logit连接 是logistics回归的形式
log<-glm(formula, family = binominal(link = logit),data = data.frame)
#Possion分布
log<-glm(formula, family = poisson(link = log),data = data.frame)
#Gamma分布
log<-glm(formula, family = gamma(link = inverse),data = data.frame)

More Reference Functions

The extended Logistic returns and variants in R are as follows: Icon

• The glmRob() function in the robust Logistic returns robust package can be used to develop broadly linear models that are robust, including robust Logistic returns. When the proposed Logistic regression model data are isolated and strongly influenced, a robust Logistic return can be useful.
• Multiple distribution regressions can be combined with multiple Logistic returns using the mlogit() function in the mlogit package if the response variable contains more than two disorderly categories (e.g. married/widow/divorce).
• The Logistic returns if the response variable is an orderly group (e.g., credit risk is differential/good/good), the lrm() function in the rms package is used to combine the Logistic returns.

R provides some useful extensions to the basic porcelain regression model

• We have a processing habit of converting to the following formulation model.$\log_\mathrm{e}\left(\frac{\lambda}{time}\right)=\beta_0+\sum_{j=1}^p\beta_jX_j$
• Zeroinfol() function in pscl package allows zero-inflating borose return
• The glm Rob() function in the robust package can be designed to match a robust broad linear model with a robust permersal pine regression

Non-linear regression model

Internal online sex return

The most classic internal online regression is multi-return; the normal multi-returnal deformation that can be solved by the multi-linear regression method described above, which we will now present in the next section is the multi-dimensional calculation method.

Function form is

poly(x, ..., degree = 1, coefs = NULL)
#计算正交多项式 degree是阶数

Inlinear nonlinear regression

This is about the most optimized questions. The functions are not unique, one introduction, we create what we think is a regression function, then optimize the acquisition of parameters.

nls(formula, data = parent.frame(), start, control = nls.control(), algorithm = "default", trace = FALSE, subset, weights, na.action, model = FALSE)
#nls函数对于实现内在非线性回归非常的实用

nlm(f, p, hessian = FALSE,
typsize=rep(1, length(p)), fscale=1, print.level = 0, ndigit=12, gradtol = 1e-6, stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000), steptol = 1e-6, iterlim = 100, check.analyticals = TRUE, …)
#nlm函数也可以处理这个问题 当然它本身就是用来处理最优化问题的 所以需要转化问题形式

There are examples.

cl<-data.frame(
X=c(rep(2*4:21, c(2, 4, 4, 3, 3, 2, 3, 3, 3, 3, 2, 3, 2, 1, 2, 2, 1, 1))), Y=c(0.49, 0.49, 0.48, 0.47, 0.48, 0.47, 0.46, 0.46, 0.45, 0.43, 0.45, 0.43, 0.43, 0.44, 0.43, 0.43, 0.46, 0.45, 0.42, 0.42, 0.43, 0.41, 0.41, 0.40, 0.42, 0.40, 0.40, 0.41, 0.40, 0.41, 0.41, 0.40, 0.40, 0.40, 0.38, 0.41, 0.40, 0.40, 0.41, 0.38, 0.40, 0.40, 0.39, 0.39))
nls.sol<-nls(Y~a+(0.49-a)*exp(-b*(X-8)), data=cl, start = list( a= 0.1, b = 0.01 ))
summary(nls.sol)

fn<-function(p, X, Y){
f <- Y-p[1]-(0.49-p[1])exp(-p[2](X-8))
res<-sum(f^2)
f1<- -1+exp(-p[2](X-8))
f2<- (0.49-p[1])exp(-p[2](X-8))(X-8)
J<-cbind(f1,f2)
attr(res, "gradient") <- 2t(J)%%f
res
}
#建立最优化函数
out<-nlm(fn, p=c(0.1, 0.01), X=cl$X, Y=cl$Y, hessian=TRUE); out

Multi-statistical analysis

We're here to show how common methods in multiple statistical analysis can be achieved in R.

Main ingredient distribution and factor analysis

As two classic techniques of relief, we are here to present their R-realization, actually, in a very common way. So we could have put it together.

Base R-provided algorithm

Main ingredient analysis

#PCA的计算
princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,subset = rep(TRUE, nrow(as.matrix(x))), ...)

#提取主成分信息
summary(object, loadings = FALSE, cutoff = 0.1, …)
#分析载荷矩阵
loadings(x)
#预测新数据主成分的值
predict(object, newdata, …)
#绘制主成分的碎石图
screeplot(x, npcs = min(10, length(x$sdev)),type = c("barplot", "lines"), main = deparse(substitute(x)), …)
#绘制数据关于主成分的散点图
biplot(x, choices = 1:2, scale = 1, pc.biplot = FALSE, …)

• xIt's data for the main ingredient analysis.
• cor T and F determine whether to use the sample for the main component analysis or the matrix for the main ingredient analysis.

Factor analysis

The function of the factor analysis is

factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
         subset, na.action, start = NULL, scores = c("none", "regression", "Bartlett"),
         rotation = "varimax", control = NULL, ...)

• xis the data, as expressed in the data box
• factorsMeaning factor
• scoresThis means that you have to use the factor score.
• rotation = "varimax"This means rotate with the maximum variance

The analysis of the factor analysis is basically the same as that of the main ingredient, with differential contribution rates, and the load matrix, which is available for analysis, has all been stored in a list form for subsequent analysis.

Functions provided by psych packages

The function here has a slightly higher degree of freedom and a slightly more detailed information than the underlying R, and they combine the function system.

#含多种可选的方差旋转方法的主成分分析
principal()

#可用主轴、最小残差、加权最小平方或最大似然法估计的因子分析
fa()
#含平行分析的碎石图（做随机数据矩阵相应的平均特征值，辅助选择主成分个数或辅助因子分析的进行，毕竟两者本质接近）
fa.parallel()
#绘制因子分析或主成分分析的结果
factor.plot()
#绘制因子分析或主成分的载荷矩阵
fa.diagram()
#因子分析和主成分分析的碎石图
scree()

Disaggregation

The most common method of diagnosing is distance and Fisher.

Fisher says goodbye.

It's basically a linear LDA method of classification, which is generally called Fisher's method in multiple statistics, and it's a function of LDA.

The function is as follows:

lda(formula, data, ... , subset, na.action)

• Formula is the formula for the sorting of the variable that is described as the return of the source to the classification.
• Subset indicates training samples

Use iris data sets as an example

data(iris)
attach(iris)
names(iris)
library(MASS)
iris.lda <- lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width)
iris.lda
iris.pred=predict(iris.lda)$class
#用于预测 此时预测的是训练用数据集
table(iris.pred, Species)
detach(iris)

The final predictive matrix shows the difference and the original situation.

Distance-segment

We know that the core of distance determination is the calculation of distance, and so R does not design a function for distance separation, but there are functions that help us calculate the distance of the horse.

mahalanobis(x, center, cov, inverted=FALSE, ...)

It accepts a data box as input, and center is the data centre cov is the matrix of the co-ordinated array, which is output in the form of a matrix that reflects the distance between the two elements.

Cluster analysis

System Cluster

There are very simple functions for system cluster analysis that can help us do this.

#计算距离矩阵使用
dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)

#计算聚类结果用
hclust(d, method = "complete", members=NULL)
#绘制聚类图 它有着聚类专用形式为
plot(object, hang=-1)
#用来对聚类结果进行切割 给出我们需要的类个数或者高度就可以了
plclust(object, hang=-1)
rect.hclust(tree, k = NULL, which = NULL, x = NULL, h = NULL,border = 2, cluster = NULL)
#聚类结果转化为树状的谱系图 使用plot绘制
as.dendrogram(object, hang = -1, …)

• The D is the distance structure, the method is the choice of the system cluster method, the default maximum distance.
• x is the data box
• Method is the method of calculating distance.
• The diag andupper logical variables control whether the output is only diagonal or the output is up triangle.

Dynamic Cluster

The classic dynamic cluster method is k-means.

kmeans(x, centers, iter.max = 10, nstart = 1,
algorithm = c("Hartigan-Wong", "Lloyd", "Forgy", "MacQueen"))

• x is the data box
• Canters are a group number or an initial cluster centre.
• It's the maximum number of words.
• Algorithm is an algorithm for dynamic clustering.

Number of clusters

One of the most desirable types of clustering is the number of types we need to group samples, which, if too few, may place data of a serious heterogeneity in the same group, or if too many, may not be fully classified, and the real cluster scenario may well require some knowledge of the field. We're here to present some stable methods of selecting the number of categories.

NbClustThe package provides a large number of indicators to determine the optimal number of categories in a cluster analysis. There is no guarantee that the results of these indicators will be consistent. In fact, they may be different. However, the results can be used as a reference for selecting the group of K-digit values.NbClust()The function input includes a matrix or data frame that needs to be used for a cluster, the distance measure and the group method used, and the number of the smallest and largest grouping to be used for a cluster. It returns each cluster index, and it also outputs the optimal number of recommended clusters.

An example of a code is

library(NbClust)
nc <- NbClust(nutrient.scaled, distance="euclidean",
                  min.nc=2, max.nc=15, method="average")
#返回结果包含了各种指标决定的聚类类别个数，以及他们的投票结果

Typical relevant analysis

The typical relevant analysis is a multi-dimensional statistical approach to the relationship between the two sets of variables, and the only way we can learn to express the correlation between the two groups in as simple a form as possible.

function in the

cancor(x, y, xcenter = TRUE, ycenter = TRUE)

• Where x y is a data matrix of two variables

Results include:

• Typical correlation factor
• The payload factor used to construct typical correlations

The typical relevance of the procedure is also the hypothetical test.

Corresponding analysis

Function is in the form of a MASS package

corresp(x, nf = 1, ...)

• x is the form of the data matrix.
• The nf is the factor.

One simple example is

x.df=data.frame(HighlyFor=c(2, 6, 41, 72, 24), For =c(17, 65, 220, 224, 61), Against=c(17, 79, 327, 503, 300), HighlyAgainst=c(5, 6, 48, 47, 41))
rownames(x.df)<-c("BelowPrimary", "Primary", "Secondary", "HighSchool","College")
biplot(corresp(x.df, nf=2))
#最后这是绘制了对应分析图 怎么分析我们在理论研究的时候证明过了